We present an efficient time-stepping scheme for simulations of the coupled Navier–Stokes Cahn–Hilliard equations for the phase field approach. The scheme has several attractive characteristics: (i) it is suitable for large density ratios, and numerical experiments with density ratios up to 1000 have been presented; (ii) it involves only constant (time-independent) coefficient matrices for all flow variables, which can be pre-computed during pre-processing, so it effectively overcomes the performance bottleneck induced by variable coefficient matrices associated with the variable density and variable viscosity; (iii) it completely de-couples the computations of the velocity, pressure, and the phase field function. Strategy for spectral-element type spatial discretizations to overcome the difficulty associated with the large spatial order of the Cahn–Hilliard equation is also discussed. Ample numerical simulations demonstrate that the current algorithm, together with the Navier–Stokes Cahn–Hilliard phase field approach, is an efficient and effective method for studying two-phase flows involving large density ratios, moving contact lines, and interfacial topology changes.