We present radiation hydrodynamics simulations of the collapse of massive pre-stellar cores .
We treat frequency dependent radiative feedback from stellar evolution and accretion luminosity at a numerical resolution down to 1.27 AU .
In the 2D approximation of axially symmetric simulations , it is possible for the first time to simulate the whole accretion phase of several 10 ^ { 5 } yr for the forming massive star and to perform a comprehensive scan of the parameter space .
Our simulation series show evidently the necessity to incorporate the dust sublimation front to preserve the high shielding property of massive accretion disks .
Our disk accretion models show a persistent high anisotropy of the corresponding thermal radiation field , yielding to the growth of the highest-mass stars ever formed in multi-dimensional radiation hydrodynamics simulations .
Non-axially symmetric effects are not necessary to sustain accretion .
The radiation pressure launches a stable bipolar outflow , which grows in angle with time as presumed from observations .
For an initial mass of the pre-stellar host core of 60 , 120 , 240 , and 480 \mbox { M } _ { \odot } the masses of the final stars formed in our simulations add up to 28.2 , 56.5 , 92.6 , and at least 137.2 \mbox { M } _ { \odot } respectively .