In this paper, a fractional step lattice Boltzmann method is proposed to model two-phase flow with large density differences by solving a modified Cahn-Hilliard phase-field equation and the incompressible Navier-Stokes equations. In order to maintain the interface profile and conserve the volume, the normal flux in the original Cahn-Hilliard equation is removed and a profile-corrected term to enforce the interface to be a hyperbolic tangent profile is added. The modified Cahn-Hilliard equation is split to two sub-equations. One is the original CH equation and the other is the correction equation. Both the original CH and the hydrodynamic equations are solved by the lattice Boltzmann method. The correction equation is solved by the finite difference method. With the multi-relaxation-time collision model and a high-order compact selective filter operation employed in the lattice Boltzmann method, numerical stability is much improved. Compared with the previous lattice Boltzmann methods, the proposed method is able to maintain the order parameter within a physically meaningful range, which is capable of tracking the interface accurately. The proposed method can also simulate two-phase fluid flows with the density ratio up to 1000. Several benchmark problems, including single vortex deform of a circle, translation of a drop, stationary droplet, capillary wave and rising bubble with large density ratios, are presented to validate the accuracy and efficiency of the present method. (C) 2019 Elsevier Ltd. All rights reserved.