Direct numerical simulations of two-dimensional (2D) and three-dimensional (3D), single-mode and multi-mode, incompressible immiscible Rayleigh–Taylor (RT) instabilities are performed using a phase-field approach and high-order finite-difference schemes. Various combinations of Atwood number, Reynolds number, surface tension, and initial perturbation amplitude are investigated. It is found that at high Reynolds numbers, the surface tension, if significant, could prevent the formation of Kelvin–Helmholtz type instabilities within the bubble region. A relationship is proposed for the vertical distance of the bubble and spike vs the Atwood number. The spike and bubble reaccelerate after reaching a temporary plateau due to the reduction of the friction drag as a result of the formation of the spike vortices and also the formation of a momentum jet traveling upward within the bubble region. The interface for a 3D single-mode instability grows exponentially; however, a higher Reynolds number and/or a lower Atwood number could result in a noticeably larger surface area after the initial growth. It is also shown that a 3D multi-mode RT instability initially displays an exponential interface growth rate similar to single-mode RT instabilities. Due to the collapse and merging of individual single-mode instabilities, the interface area for a multi-mode RT instability is strongly dependent to the mesh resolution after the exponential growth rate. However, the ratio of kinetic energy over released potential energy exhibits an almost steady state after the initial exponential growth, with values around 0.4, independently of the mesh resolution.