© R. Mignon-Risse et al. 2020

1 Introduction

Massive stars shape the dynamical and chemical evolution of galaxies because of their powerful feedback in radiation, winds, explosions in supernova, and metal-enrichment. However, their formation remains a long-standing problem. Observationally, massive stars are embedded in dense clouds; they form on timescales much shorter than their low-mass counterpart (Motte et al. 2007) and are likely to be located at distances larger than 1 kpc from us, which makes their formation process challenging to observe. Two major scenarios are under active studies: the competitive accretion model and the turbulent core accretion model. In the competitive accretion model (Bonnell et al. 2004), all stars form in clusters and stars located at the center of the gravitational potential gain more mass and eventually become massive stars, via accretion and possibly merging processes. In this scenario, the initial-mass function (IMF) is built-up naturally. On the other hand, the turbulent core model (McKee & Tan 2003) is an extension of the low-mass star formation scenario. A massive and turbulent prestellar core gravitationally collapses while fragmentation is limited by turbulent, radiative, and magnetic support (e.g. Commerçon et al. 2011a). The IMF is therefore linked to the prestellar core mass function since more massive stars form in more massive cores. For reviews on theories of massive star formation, we refer the reader to Beuther et al. (2007), Zinnecker & Yorke (2007), Tan et al. (2014), and Krumholz (2016).

There is no consensus regarding the accretion process and a way to probe it is to study outflows. Indeed, magnetic outflows are often associated with accretion (Blandford & Payne 1982; Pelletier & Pudritz 1992) as they remove the angular momentum from an accretion disk and they do not require strong magnetic fields. It is then possible to study the physics of accretion via the outflow properties (Pudritz et al. 2006) and in particular the accretion rate from the outflow velocity (Pelletier & Pudritz 1992). In the case of massive stars, they can act as a channel for radiation to escape. Moreover, accretion modes differ from one scenario to another. Disk accretion is more likely to occur in the turbulent core accretion model, with high accretion rates (Ṁ ~ 10⁻⁴−10⁻³ M_⊙ yr⁻¹; McKee & Tan 2003). More chaotic accretion mechanisms are associated with the competitive accretion model or with models including accretion via filaments. Recent observations by Goddi et al. (2018) reveal signatures of outflows whose direction varies through time. If they are perpendicular to the accretion disk (Blandford & Payne 1982) this indicates that the plane of accretion changes with time, favoring accretion via filaments and competitive accretion. The study of outflows can indeed help to distinguishing the accretion modes and these questions highlight the need for realistic outflow models (magnetic and radiative) in numerical simulations. Regarding the radiative outflows, this means the use of a radiative transport method well-suited in the optically-thin regime. Current and past studies of massive star formation have mainly focused on its radiative aspects. In a 1D spherically-symmetric approach, the radiative force of a massive (proto)star is expected to counteract gravity up to the point where accretion is stopped as computed analytically (Larson & Starrfield 1971) and then numerically (Kuiper et al. 2010a). Their results showed that the highest mass reached was 40 M_⊙. Two-dimensional and 3D numerical simulations have permitted the emergence of a new accretion mode, the “flashlight effect” (Yorke & Sonnhalter 2002) which allows the radiation to escape freely in the poles while material is accreted through the disk (Krumholz et al. 2009; Kuiper et al. 2010a; Rosen et al. 2016; Harries et al. 2017).

Monte-Carlo approaches are often used for solving radiative transfer problems for their accuracy but they are particularly expensive and their computational time scales with the number of radiative sources. This justifies the use of fluid description models suchas the flux-limited diffusion (FLD) and the M1 methods for radiation-hydrodynamics (RHD). The first RHD calculations relied on the FLD closure relation for its simplicity and advantageous computational cost. However, it is more suited for the optically-thick regime while the use of a flux-limiter corrects the propagation speed in the optically-thin regime (Levermore & Pomraning 1981). In addition, the FLD method does not permit to capture shadows behind very optically-thick gas. In the context of massive star formation, the flashlight effect is due to the nonisotropic character of the radiation field because the optical thickness is very different in the disk direction and in the cavities direction. Therefore, numerical developments regarding radiative transfer have been made, especially in the optically-thin limit. Recent approaches treat stellar irradiation in a more consistent way with ray-tracing (Kuiper et al. 2010b; Kim et al. 2017), long-characteristics (Rosen et al. 2017), Monte-Carlo radiative transfer (Haworth & Harries 2012; Harries et al. 2017, who took advantage of the independency between photon packets to parallelize efficiently the Monte-Carlo step), and the M1 closure relation (Levermore 1984; González et al. 2007; Aubert & Teyssier 2008; Rosdahl et al. 2013; Kannan et al. 2019; this work).

Multi-dimensional simulations using the FLD approximation or hybrid approaches show stars with mass above the 40 M_⊙ limit obtained in 1D: with the FLD method only, Yorke & Sonnhalter (2002) and Krumholz et al. (2009) form a star of ≈ 42 M_⊙ from a 120 M_⊙ and 100 M_⊙ prestellar core, respectively. The additional treatment of direct irradiation in hybrid approaches has been shown not to impact the stellar mass significantly: Klassen et al. (2016) have obtained stars as massive as 43.7 M_⊙ from an initial mass of 100 M_⊙, the simulations of Rosen et al. (2016) show a 40 M_⊙ star and Kuiper et al. (2010a) obtain a 56.5 M_⊙ star from a 120 M_⊙ prestellar core in several free-fall times. Most of these works put lower-limit on the stellar mass because the accretion phase is not finished yet at the end of the run (except for Kuiper et al. 2010a). These investigations have also noted the formation of polar cavities dominated by the stellar radiative pressure, enhanced by the particular treatment of stellar irradiation.

In addition, Krumholz et al. (2009) and Rosen et al. (2016) observe the onset of radiative Rayleigh–Taylor instabilities in the radiation-pressure-dominated cavities that feed the star-disk system and help accreting mass onto the central star via the flashlight effect. However, these instabilities have not been observed in the work of Kuiper et al. (2010a) and Klassen et al. (2016). Kuiper et al. (2012) argue that a gray FLD model, as used in Krumholz et al. (2009), underestimates the radiation force in the cavity and can artificially lead to the apparition of these instabilities. With a frequency-dependent hybrid model and a Cartesian grid, Klassen et al. (2016) did not obtain such instabilities while Rosen et al. (2016) did. The difference in their results can be explained by the use of refinement by Rosen et al. (2016) to resolve the seeds of the instabilities, the smallest modes being the more unstable (Jacquet & Krumholz 2011). Contrarily, the spherical grid without additional refinement in the cavities used by Kuiper et al. (2012) would not permit to refine them. It is thus unclear yet whether this mechanism is at work during the formation of massive stars. Meanwhile, it is clear that disk-accretion is sufficient to reach masses consistent with the massive stars observed. Our hybrid method, implemented in the Cartesian adaptive-mesh refinement (AMR) code RAMSES (Teyssier 2002) help us to establish the importance of these accretion mechanisms by capturing the nonisotropy of the radiation field.

This paper is organized as follows. Section 2 presents the equations of the flux-limited diffusion method and the M1 method, along with their coupling and implementation in the RAMSESxd code. We present the tests and the validation of our hybrid approach in Sect. 3 and its application to the collapse of a massive prestellar core leading to the formation of a massive star in Sect. 4. We discuss our results in Sect. 5

2 Methods

In this section, we present the equations of the hybrid radiative transfer approach and its implementation: the M1 method for the stellar irradiation and the flux-limited diffusion method for the dust emission.

2.1 Coupling flux-limited diffusion and M1

The FLD method (Levermore & Pomraning 1981) and the M1 method (Levermore 1984) are fluid descriptions of the radiation field. They are based on moments of the equation of radiative transfer, that is, the equation of conservation of the radiation specific intensity I_ν (x, t;n) with the propagation, the absorption and the emission (Mihalas 1984) $$\frac{1}{c}\frac{\partial I_{\nu }}{\partial t}+n\cdot \nabla I_{\nu }=-{\kappa }_{\nu }\rho I_{\nu }+{\eta }_{\nu }.$$(1)

Here, I_ν(x, t;n) is the amount of energy of a photon beam at a given position x and time t, in direction n and per unit frequency. c is the speed of light, κ_ν is the extinction coefficient (absorption and scattering contributions), ρ is the local density and η_ν is the emission coefficient. The sum of dust and gas contributions to the medium opacity weighted with the dust-to-gas ratio are encapsuled in κ_ν. We assume local thermodynamical equilibrium (LTE) and do not consider scattering (see the justification in Appendix A), hence the emission coefficient η_ν is a source function proportional to the Planck function B_ν(T) and the extinction coefficient κ_ν is just an absorption coefficient. Equation (1) must be solved at each hydrodynamical timestep to evolve the radiation field, which still depends on six variables (x, n, ν). In addition, we want to couple it to hydrodynamics. This motivates the need for taking moments of the equation of radiative transfer. Hence we lose some of the angular information but it reduces the number of variables to four, at each timestep. The radiative energy density E_ν, the radiative flux F_ν and the radiative pressure tensor ${\mathbb{P}}_{\nu }$ are defined as the 0th, 1st and 2nd moments of the radiative intensity I_ν, respectively.

Each system of moment equations involves the ith and the (i + 1)th moments, hence we need a closure relation. The FLD scheme is based on the diffusion approximation, which is suited for high optical depths, when photons propagate in a random walk in the material (e.g., in stellar interiors). It has been commonly used as a first step to introduce radiation into hydrodynamical codes (Krumholz et al. 2007; Kuiper et al. 2008; Tomida et al. 2010). In the FLD model, the equation to be solved is theequation of conservation of the radiative energy. Once integrated over all frequencies (often called a gray approximation) it gives $$\frac{\partial E_{\text{r}}}{\partial t}-\nabla \cdot \left(\frac{c\lambda }{{\kappa }_{\text{R}}\rho }\nabla E_{\text{r}}\right)={\kappa }_{\text{P}}\rho c\left(\text{a}{T}^4-E_{\text{r}}\right),$$(2)

where E_r is the frequency-integrated radiative energy, λ is the flux-limiter and is built to recover the right propagation speed in optically-thin and -thick media (Levermore & Pomraning 1981). The opacities aregiven by κ_P and κ_R, which are are Planck’s and Rosseland’s mean opacities, respectively. Thermal radiation is modeled as the Planck function B(ν, T), therefore under the gray model Planck’s and Rosseland’s opacities are respectively defined as $${\kappa }_{\text{P}}=\frac{{\displaystyle {\int }_0^{\infty }{\kappa }_{\nu }}B(\nu )\text{d}\nu }{{\displaystyle {\int }_0^{\infty }B}(\nu )\text{d}\nu },$$(3)

and $${\kappa }_{\text{R}}=\frac{{\displaystyle {\int }_0^{\infty }\frac{\partial B(\nu ,T)}{\partial T}}\text{d}\nu }{{\displaystyle {\int }_0^{\infty }\frac{1}{{\kappa }_{\nu }}}\frac{\partial B(\nu ,T)}{\partial T}\text{d}\nu },$$(4)

where the temperature derivatives appear from the chain rule applied to ∇B. The a T⁴ term in Eq. (2) arises from the integral of the Planck function over all frequencies, with a the radiation constant. We approximate the mean opacity of the radiative energy term as the Planck mean opacity.

On the other hand, with the M1 method we take the zeroth and first moments of the equation of radiative transfer (Levermore 1984). Within the M1 method we obtain the following system for the radiative energy and flux conservation, in the gray approximation as well $$\frac{\partial E_{\text{r}}}{\partial t}+\nabla \cdot F_{\text{r}}=-\rho {\kappa }_{\text{P}}\text{c}{E}_{\text{r}}+{\dot{E}}_{\text{r}}^{\star },\text{nn \&}\frac{\partial F_{\text{r}}}{\partial t}+{\text{c}}^2\nabla \cdot {\mathbb{P}}_{\text{r}}=-\rho {\kappa }_{\text{P}}\text{c}{F}_{\text{r}},$$(5)

where ${\dot{E}}_{\text{r}}^{\star }$ is the rate of radiative energy injected from stellar sources. F_r and ${\mathbb{P}}_{\text{r}}$ are the frequency-integrated radiative flux and pressure respectively.

One main asset is that the directionality of the photons beam is well-retained. The M1 method is able to model shadows to some extent (see González et al. 2007), in an irradiated accretion disk for instance, while FLD is not. In addition, as a moment method the computing cost does not scale with the number of sources.

Our goal is to take advantage of both methods, that is, the FLD method for an optically-thick medium and M1 for irradiation. Both FLD and M1 methods described above can involve several groups of photons or only one with frequency-averaged opacities. In our study however we restrict ourselves to one group of photons treated with each method.

In massive star formation simulations the dynamical influence exerted by the radiative feedback is of main importance as well as the thermal structure of the accretion flow (e.g., for fragmentation). However, doing so requires to retain to some extent the directionality of the photons emitted by the star to compute the direct radiative force. Breaking this isotropy is consistent with probing nonisotropic modes of accretion (disk or filaments). Secondly, it requires to distinguish the opacities between stellar photons, which have a UV-like energy and relatively high opacities, and photons emitted by the dust, which have a IR-like energy and relatively low opacities.

Our method is to inject the stellar photons into the group of photons treated with the M1 scheme. The gray opacity used with the M1 corresponds to the Planck mean opacity at the stellar temperature, κ_P (T_⋆), written κ_P,⋆ for the sake of readability. Once these photons are absorbed by the medium they are depleted from the M1 group as they heat the gas. The gas reemission is treated with the FLD method. In a first approach we do not deal with ionization states and leave this to further work. The set of equations that are to be solved are $$\frac{\partial E_{\text{M1}}}{\partial t}+\nabla \cdot F_{\text{M1}}=-{\kappa }_{\text{P},\star }\rho \text{c}{E}_{\text{M1}}+{\dot{E}}_{\text{M1}}^{\star },\text{nn \&}\frac{\partial F_{\text{M1}}}{\partial t}+{\text{c}}^2\nabla \cdot {\mathbb{P}}_{\text{M1}}=-{\kappa }_{\text{P},\star }\rho \text{c}{F}_{\text{M1}},\text{nn \&}\frac{\partial E_{\text{fld}}}{\partial t}-\nabla \cdot \left(\frac{c\lambda }{{\kappa }_{\text{R},\text{fld}}}\nabla E_{\text{fld}}\right)={\kappa }_{\text{P},\text{fld}}\rho \text{c}\left(\text{a}{T}^4-E_{\text{fld}}\right),\text{nn \&}{C}_{\text{v}}\frac{\partial T}{\partial t}={\kappa }_{\text{P},\star }\rho \text{c}{E}_{\text{M1}}+{\kappa }_{\text{P},\text{fld}}\rho \text{c}\left(E_{\text{fld}}-\text{a}{T}^4\right),$$(6)

where ${\dot{E}}_{\text{M1}}^{\star }$ is the stellar radiation injection term, and κ_P,⋆ρcE_M1 couples the M1 and the FLD methods via the equation of evolution of the internal energy. We use the ideal gas relation for the internal specific energy ɛ = C_vT where C_v is the specific heat capacity at constant volume. This equation closes the system and is used to evolve the gas temperature together with the radiative quantities.

2.2 Radiative acceleration

In addition of improving the thermal coupling between stellar irradiation and gas, our implementation is meant to affect the gas dynamics via a more accurate and less isotropic approach for the radiative acceleration than the FLD approximation thanks to the equation of evolution of the stellar radiative flux. We are interested in comparing the radiative acceleration with the hybrid method and with the pure FLD method.

In the frame of RHD (for the full expression of RHD equations we refer the reader to Mihalas 1984), the radiative acceleration at a given frequency is equal to $$a_{\text{rad},\nu }={\kappa }_{\nu }{F}_{\nu }/\text{c}.$$(7)

However, gray FLD and M1 methods do not share the same expression for the radiative flux. On one hand, the M1 model includes the Planck mean opacity as the flux-averaged opacity – which means that momentum is transferred each time a photon is absorbed – so the radiative acceleration is given by $$a_{\text{rad},\text{M1}}={\kappa }_{\text{P},\star }{F}_{\text{M1}}/\text{c}.$$(8)

On the other hand, the gray radiative acceleration in the FLD approximation is given by $$a_{\text{rad},\text{fld}}=-\frac{\lambda }{\rho }\nabla E_{\text{fld}}.$$(9)

Taking the asymptotic values of λ into account (see Levermore 1984), we get $\Vert a_{\text{rad},\text{thin}}\Vert ={\kappa }_{\text{R}}{E}_{\text{fld}}$ in the limits of low optical depth and $a_{\text{rad},\text{thick}}=\frac{1}{3\rho }\nabla E_{\text{fld}}$ for high optical depth.

We recall that κ_R is a harmonic mean which favorslow absorption bands while κ_P is an arithmetic mean which favors high absorption bands. As a consequence, we expect a higher radiative acceleration with the hybrid method than with the FLD method.

2.3 Implementation

The RAMSES code (Teyssier 2002) is a 3D adaptive-mesh refinement Eulerian code. We use a version of RAMSES which has been widely used for star formation simulations (Commerçon et al. 2011a, 2012; Joos et al. 2012; Hennebelle et al. 2016; Vaytet et al. 2018).

The hydrodynamical solver of RAMSES relies on finite volume methods (variables are volume-averaged over the cell) and a second-order Godunov method is used to evolve hydrodynamical variables. This code includes theFLD (Commerçon et al. 2011b, C11 hereafter) and the M1 method within RAMSES-RT (Rosdahl et al. 2013¹, R13 hereafter), both coupled to the hydrodynamics. For the FLD method, a time-splitting approach is performed. A predictor-corrector MUSCL scheme is used, where the predictor step is made under the diffusion approximation so the radiative pressure is isotropic and nonisotropy is taken into account in the corrector step. The hyperbolic part of the FLD solver relies on the second-order Godunov scheme of RAMSES and fluxes are estimated with an approximate Riemann solver (Lax-Friedrich, HLL, HLLD, etc.). The diffusion and radiation-matter coupling are handled in the implicit part of the time-splitting scheme. The diffusion part of the FLD solver is second-order accurate in space. The M1 module estimates fluxes with a Riemann solver (HLL or GLF). In this work we use the GLF solver because it captures better the isotropy of stellar radiation than HLL, and the reduced fluxapproximation for the direct radiative force (see Appendix B of Rosdahl et al. 2015 and discussion in Hopkins & Grudić 2019).

The radiative transfer puts a heavy constraint on the timestep because the Courant-Friedrich-Lewy (CFL) condition forbids the propagation of signals through more than one cell in one timestep, for explicit schemes. For the hydrodynamics this speed is the sound speed but for radiative transfer this is the speed of light and this would impose a timestep ~ 1000 times shorter. Therefore, both C11 with the FLD method and R13 with the M1 scheme have used a workaround.

On one hand, the FLD solver for diffusion and radiation-matter coupling is implicit and therefore is unconditionally stable. The system composed of the internal energy and the radiative energy equations in their discretized form leads to the inversion of a matrix computed with the conjugate gradient or biconjugate gradient algorithm, in the case of multigroup radiative transfer (González et al. 2015). The number of iterations to converge scales with the number of cells. We note that the isotropic radiative pressure contributes to the total pressure in the explicit solver of RAMSES, which therefore must satisfy the CFL condition. As a consequence, the radiative pressure must be taken into account when computing the timestep allowed by the CFL condition (Commerçon et al. 2011b).

On the other hand, the M1 solver is fully explicit with a first-order Godunov scheme, thus it obeys the CFL condition. The trick used here is the reduced speed-of-light approximation (Gnedin & Abel 2001). In this approximation, the propagation of light is not restricted by the speed of light but by the speed of the fastest wave, which is the speed of ionization front in the original paper and the fastest hydrodynamical speed in our case. An additional subcycling method relaxes this constraint. This leads to a timestep set by the hydrodynamical CFL condition.

The injection of energy from the stellar source into the M1 photons is made via a sink algorithm (Bleuler & Teyssier 2014). In this work, we retrict ourselves to one stellar source for the M1 photons. The M1 module ensures the propagation and absorption of the stellar radiation while the FLD module deals with the heating by the stellar radiation and treats the reemission. We test the accuracy of the hybrid approach with respect to the FLD module alone, since it was used in previous massive star formation calculations with the RAMSES code (Commerçon et al. 2011a).

3 Numerical tests

We test the hybrid method in a pure radiative transfer case (i.e., no hydrodynamics): a static disk irradiated by a star. We compare the temperature structure obtained with results from Monte-Carlo radiative transfer codes. We explore three levels of optical thickness integrated along the disk mid-plane: τ = 0.1, τ = 100 and τ = 10³. The parameters and results are summarized in Table 1. We refer the reader to Appendix B for performance tests.

3.1 Optically-thin and moderately optically-thick regimes: Pascucci’s test

3.1.1 Physical and numerical configurations

The first test we have performed is taken from Pascucci et al. (2004) and consists of a star irradiating a static disk made of dust and gas. Weuse it to probe the behavior and accuracy of our method in the optically-thin and moderately optically-thick regimes. In particular, we compare our results once the temperature structure is converged with respect to time with those obtained with Monte-Carlo RT codes such as RADMC-3D (Dullemond et al. 2012) and MCFOST (Pinte et al. 2006).

This is a 2D test of a static flared disk of a given analytical profile for the gas density, depending on the cylindrical radius r and on the vertical height z. The disk extends from r_in = 1 AU to r_out = 1000 AU. The density ρ(r, z) in cylindrical coordinates is given as $$\rho (r,z)={\rho }_0{f}_{\text{1}}(r)f_{\text{2}}(r,z),$$(10)

where ρ₀ is the density normalization and is linked to the only free-parameter, the integrated optical-depth throughout the mid-plane of the disk, ${\tau }_{\nu }={\displaystyle {\int }_{{\text{r}}_{\text{in}}}^{{\text{r}}_{\text{out}}}{\kappa }_{\nu }}\rho (r,z=0)\text{d}r$. The two functions f₁ and f₂ are given by $$f_1(r)={\left(\frac{r}{r_{\text{d}}}\right)}^{-1},$$(11)

and $$f_2(r,z)=\text{exp}\left(-\frac{\pi }{4}{\left(\frac{z}{h(r)}\right)}^2\right),$$(12)

where the flaring function is $$h(r)=z_{\text{d}}{\left(\frac{r}{r_{\text{d}}}\right)}^{1.125}.$$(13)

In this setup, r_d = r_out∕2 = 500 AU and z_d = r_out∕8 = 125 AU are the scale-radius and the scale-height. The star is not resolved but its luminosity is based on its physical radius and temperature. In this test, it has a radius R_⋆ =1 R_⊙ and can have two possible surface temperature: T_⋆,1 = 5800 K and T_⋆,2 = 15 000 K.

The integrated optical depth (for extinction, as in the literature) is taken to be either τ = 0.1 or τ = 100 at 550 nm to probe the optically-thin and moderately optically-thick regimes, respectively. The dust-to-gas mass ratio is equal to 0.01. Dust is made of spherical astronomical silicates of radius 0.12 micron and density of 3.6  g cm⁻³. Frequency-dependent dust opacities are taken from Draine & Lee (1984) as in Pascucci et al. (2004) and are displayed in Fig. 1. In these setups we only take the absorption into account and do not consider scattering. The corresponding Planck and Rosseland mean opacities used in the gray M1 and FLD modules are displayed in Fig. 2. We recall that we take the M1 absorption coefficient as the Planck mean opacity at the stellar temperature, κ_P (T_⋆).

Boundary conditions are chosen to be a fixed temperature of 14.8 K and a density floor of 10⁻²³ g cm⁻³. The same density floor is applied between the star and the disk edge to mimic the vacuum that RADMC-3D and MCFOST strictly apply since their respective cylindrical and spherical grids begin at R_in. The 14.8 K temperatureis applied throughout the computational domain as initial condition and is at equilibrium with radiation.

We run the simulations with AMR levels between 5 and 14, which results in a finest resolution of Δx = 0.12 AU where Δx is the cell width. This makes possible to have several (≈9) cells between the star and the disk edge and the star to have a negligible size with respect to the disk thickness (≈ 0.01 AU against ≈ 0.04 AU for the disk height at r_in). Secondly, it permits to resolve several times the mean free-path at the disk inner edge: the local optical depth is κ_P ρ_max Δx ≈ 0.15 < 1, where ρ_max is the density at the disk inner edge for the case τ = 100. Refinement is performed on the density gradient so that the disk inner edge is at the highest refinement level. We consider that the temperature structure is converged when the relative change between successive outputs decreases below 10⁻⁴ (see Ramsey & Dullemond 2015).

Fig. 1 Frequency-dependent opacities and blackbody spectra for Tdisk = 300 K and T⋆ = 5800 K. Opacities are absorption (blue pluses), scattering (red crosses) and extinction (black dots) coefficients for the dust-and-gas mixture used in the Pasccuci test. The table contains 61 frequency bins and dataare taken from Draine & Lee (1984). Apart from the broad opacity features at about 10 and 20 μm, which correspond to Si-O vibrational transitions, the opacity generally increases with the photon frequency. The opacity at stellar-like radiation frequencies is generally greater than at disk-like radiation frequencies.

Fig. 2 Planck’s (blue dashed curve) and Rosseland’s (orange dot-dashed curve) mean opacities, as a function of temperaturein the Pascucci setup.

3.1.2 Temperature structure

The RAMSES grid is Cartesian while the grids of MCFOST and RADMC-3D are cylindrical and spherical, respectively. Therefore, we interpolate temperature values on their grids to compute the relative error at the location on the RAMSES grid. Figure 3 plots the gas temperature in the disk mid-plane against the x-axis for the most optically-thin case, τ = 0.1, once the temperature structure is converged with respect to time. The location is given by the distance to the disk inner edge, r − r_in. For T_⋆,1 and T_⋆,2 the FLD run produces an important error throughout the disk mid-plane, up to ≈ 62% and ≈ 65%, respectively, and always underestimates the temperature. On the opposite, the hybrid method is quite accurate with a maximal error of ≈ 2% while RT codes (MCFOST and RADMC-3D) agree within 1% in this test (in accord with Pascucci et al. 2004). This important difference between FLD and hybrid methods comes from the regime of validity of each method: the FLD is not well-suited for optically-thin media. In the hybrid method, since there is not much absorption because of the low optical-depth, the M1 module is mainly at work and is adapted to optically-thin media (as tested in Rosdahl et al. 2013), which justifies its good accuracy.

For direct irradiation, the Planck mean opacity in the FLD implementation is computed at the local temperature even though the radiation has been emitted by the star. A direct consequence is that the stellar radiation is absorbed by the disk with an opacity coefficient computed at the disk temperature, which is much lower than the stellar temperature. As the opacity increases with the temperature (see Fig. 2), the absorption opacity with the FLD is lower and hence the temperature is lower than that given by RT codes and by the hybrid approach. The situation worsens from T_⋆,1 to T_⋆,2 at the disk edge and the error increases from ≈30 to ≈60%. It also illustrates the need for a better approach for treating massive stars irradiation.

The left panel of Fig. 4 shows the radial temperature profile in the disk mid-plane for the moderately optically-thick case, τ = 100, and for T_⋆,1. The error made by the hybrid method is higher than for the τ = 0.1 case and reaches a maximal value of ≈25% whereas the FLD method alone makes a maximal error of ≈36%. Also, the error made by the hybrid method is quite uniform with respect to the error made by the FLD method alone. For T_⋆,1 and T_⋆,2 (left and right panels of Figs. 3 and 4, respectively), the FLD method underestimates the temperature between the star and the disk edge because the mediumis optically thin and because of the Planck opacity considered, as explained above. For T_⋆,2, both methods converge toward a similar temperature at large radii. Absorption is stronger here than in the most optically-thin case, and stronger than for T_⋆,1 (see Fig. 2), so the M1 photons are quickly absorbed and the FLD module of our hybrid method is at work. Therefore a significant error is expected from the gray opacity employed in the FLD module of the hybrid method. Indeed, the temperature varies significantly throughout the disk (between ≈ 20 and 350 K for T_⋆,1 and between ≈20 and 1000 K for T_⋆,2). As a consequence, a disk cell is crossed by photons of very different frequencies and the gray approach induces errors. Conversely, the frequency-dependence of RADMC-3D method permits one to distinguish the photons that are quickly absorbed (the most energetic ones) from those at lower energy that penetrate the disk more deeply and contribute to the disk heating at larger radii.

To examine the behavior of both methods with respect to the nonisotropy of the setup, we plotted the vertical temperature profile at a cylindrical radius of 20 AU (Fig. 5). This visualization is important for this type of tests, because an optically-thick disk produces self-shielding in the mid-plane and we expect the hybrid method to capture it better than the FLD method (González et al. 2007). Here we take the temperature given by MCFOST rather than RADMC-3D because its grid is cylindrical (and not spherical) and thus errors of interpolation are avoided. The left panel of Fig. 5 shows the temperature profile for the most optically-thin case, τ = 0.1. No self-shielding is expected and the temperature should decrease slowly as the vertical height increases. Such a behavior is obtained with MCFOST as well as with the hybrid method. The temperature obtained with the FLD method is uniform with z, which is likely due to the isotropic nature of the FLD method. The relative error is comparable to the one in the radial profile: up to ≈ 47% with the FLD method and less than 1% with the hybrid approach.

On the right panel of Fig. 5, τ = 100, and MCFOST gives a lower temperature in the mid-plane than for τ = 0.1, as expected. Conversely, the FLD method does not capture at all the nonisotropic nature of the radiation onto the irradiated disk: the temperature is fairly uniform. The hybrid method reproduces partly this feature, even though the error can be as large as ≈ 20%.

We conclude that the FLD method is not capable of reproducing the temperature profile in the optically-thin and moderately optically-thick regime. The hybrid method is very accurate in the optically-thin regime (less than ≈ 2%). In the moderately optically-thick regime, the hybrid method gives a non-negligible error (up to ≈ 31% for a 15 000 K star) in the transition between optically-thin and -thick media which shows its limitations but this is a major improvement with respect to the ≈57% error made with the FLD method. In addition, the hybrid approach captures partially (≈ 20% error) the self-shielding in the disk mid-plane while the FLD approximation does not.

Fig. 3 Radial gas temperature profiles in the mid-plane of the disk following the test of Pascucci et al. (2004) for τ = 0.1. We comparethe gas temperature computed using MCFOST (black dotted-line) and RADMC-3D (red dashed-line), the hybrid method (M1+FLD, blue dots) and the FLD method alone (orange dots) in RAMSES. Left: central star temperature T⋆,1 = 5800 K; right: T⋆,2 = 15 000 K.

Fig. 4 Same as Fig. 3, but for τ = 100.

Fig. 5 Vertical gas temperature profiles at a cylindrical radius of 20 AU, following the test of Pascucci et al. (2004). We compare the gas temperature computedusing  MCFOST, the hybrid method (M1+FLD) and FLD alone in RAMSES for T⋆,1, τ = 0.1 (left) and τ = 100 (right).

3.1.3 Impact on the radiative acceleration

We look at the radiative acceleration maps obtained with the FLD and the hybrid methods for the moderately optically-thick case (τ = 100). Figure 6 shows the radiative acceleration perpendicularly to the disk plane as obtained after temperature convergence with the FLD method (left) and the hybrid method (right). The left panel shows two peculiarities of the FLD solver. First, we recall that the FLD radiative acceleration has two asymptotic values depending on the optical regime (see Sect. 2.2): in the optically-thin limit it is proportional to the radiative energy and in the optically-thick limit it is equal to the radiative energy gradient divided by the density. Farther from the star, the disk structure is visible (the dark blue zones) because of the density dependence in the radiative acceleration. Second, the aspect of the FLD acceleration closer to the star is mainly due to grid effects.

The right panel of Fig. 6 shows the sum of the FLD and M1 radiative accelerations in the hybrid case. The combination of the optically-thin and -thick methods permits to capture the nonisotropy of the radiative acceleration. The hybrid radiative acceleration is ~ 100 greater than the FLD acceleration. This result holds in the four tests: τ = 0.1, τ = 100 and T_⋆,1, T_⋆,2. It is in agreement with the study of Owen et al. (2014). This is mainly due to the temperature at which the M1 opacity is taken. Stellar photons are at a frequency that is ~10 times greater than that of photons emitted by the surrounding gas, which implies an opacity ~ 100 greater (see Fig. 1). As shown previously, the radiation transport in the optically-thin limit is accurately treated with our hybrid approach and leads to a strong improvement for the radiative acceleration due to the direct irradiation, which is one of the main contributors expected in the dynamics of massive star formation.

3.2 Optically-thick regime: Pinte’s test

The second test is a similar but more challenging setup with a higher integrated optical depth and a sharper density profile than Pascucci et al. (2004) at the disk edge, as presented in Pinte et al. (2009). The disk extends from a cylindrical radius r_in = 0.1 AU to r_out = 400 AU and the integrated optical depth (for extinction) is τ_810nm = 10³. The flared disk density profile ρ(r, z) is analytically given by $$\rho (r,z)={\rho }_0{\left(\frac{r}{r_{\text{d}}}\right)}^{-2.625}\text{exp}\left(-\frac{1}{2}{\left(\frac{z}{h(r)}\right)}^2\right),$$(14)

where the flaring function h(r) is as before, r_d = r_out∕4 = 100 AU and z_d = r_out∕40 = 10 AU. The star has a radius R_⋆ = 2 R_⊙ and the stellar surface temperatureis T_⋆ = 4000 K.

We use the opacity table from Weingartner & Draine (2001) which gives the absorption opacity of dust grains with respect to the wavelength. These opacities were calculated for spherical astronomical silicates (see Draine & Lee 1984) of size 1 micron and density 3.6 g cm⁻³.

The sharp increase in density at the disk inner edge makes this test particularly challenging because the local variation of optical depth must be resolved while the discretized equations involve locally constant absorption opacities. At the same time, it is crucial to resolve the local mean free path to prevent an excess of photon absorption, which leads to overestimating the temperature. Resolving both is even more challenging for AMR-grid codes than for cylindrical-grid codes with no material inside R_in and a logarithmic scale. In RAMSES we choose to refine the grid based on a density gradient criterion so that the disk edge is at the finest resolution and the transition from optically-thin to -thick is as resolved as possible. There is a drawback: having the greatest resolution at the disk inner edge is very computationally expensive because, first, it affects many more cells than if the refinement is operated on the central cell (as usually done because the sink particle is located there). In addition, two adjacent cells cannot differ by more than one level of refinement and generally the number of cells at the same AMR level is much higher than two. Therefore, it also means a higher resolution at larger radii. For that reason, Ramsey & Dullemond (2015) choose not to use AMR but instead, a logarithmically-scaled grid particularly adapted to this setup. Figure 7 shows the AMR level needed to resolve the local mean free path as a function of the radius in the disk midplane along with the AMR level set with RAMSES. Hence, we perform our calculations with l_max = 22, which gives a finest cell width of 1.9 × 10⁻⁴ AU, so that the mean free path at the disk inner edge is resolved. The computational cost does not make possible to extend the zone over which the mean free path is resolved, nor to compare the temperature over the entire disk radius with RADMC-3D.

The left panel of Fig. 8 plots the radial temperature profile in the disk mid-plane obtained with RAMSES with the FLD and hybrid methods versus RADMC-3D. FLD underestimates the temperature at the disk inner edge, as in the test of Pascucci et al. (2004). The temperature slope given by the hybrid method is in a better agreement with RADMC-3D than the one given by the FLD method. For the hybrid method, the temperature at the inner edge of the disk is accurately computed (up to ≈7% error) but is overestimated at larger radii where it becomes fairly constant at ≈ 65% error. It can be seen that the error made by the hybrid approach is not negligible as the mean free path becomes unresolved (Fig. 7). The temperature profile obtained with our hybrid method is very similar to what has been obtained in comparable studies (see Fig. 8 of Ramsey & Dullemond 2015).

The right panel of Fig. 8 shows the vertical temperature profile. The temperature profile shape given by the hybrid method is similar to RADMC-3D but with self-shielding partially captured (up to ≈ 61% error), unlike in the FLD method. The hybrid method then recovers the correct temperature (≈ 2% error) at a larger disk height.

This setup highlights the need to resolve the mean free path of photons to obtain the correct temperature at the disk edge and it shows that the hybrid approach is more accurate than the FLD approximation to compute the temperature structure of an optically-thick disk. Moreover, this setup is challenging for our hybrid method because most of the direct irradiation is absorbed in the inner parts of the disk so the rest of the disk temperature structure is mainly obtained with the FLD method. As shown in Ramsey & Dullemond (2015), a frequency-dependent irradiation scheme is not more accurate in this test. However, a multigroup FLD method (González et al. 2015) would improve this, as mentioned by Ramsey & Dullemond (2015).

Fig. 6 1000 AU disk edge-on slices of the radiative acceleration, following the test of Pascucci et al. (2004) obtained with RAMSES after stationarity is reached. Left: FLD run; right: hybrid run. Star and disk parameters: τ = 100 and T⋆,1. The hybrid radiative acceleration is about 100 times greater than the FLD one.

Fig. 7 AMR level needed to resolve the mean free path (mfp) of photons (red dashed-line) and effective AMR level (black dots) in the disk midplane following the test of Pinte et al. (2009).

Fig. 8 Left: radial gas temperature profile in the mid-plane of the disk following the test of Pinte et al. (2009). Right: vertical gas temperature profile at a cylindrical radius of 0.2 AU in the disk. We compare the gas temperature computed using  RADMC-3D, the hybrid method (M1+FLD) and the FLD method alone in RAMSES. The integrated optical depth in the disk mid-plane is τ810 nm = 103 and the stellar temperature is T⋆ = 4000 K.

4 Collapse of an isolated massive prestellar core

We use the newly implemented and tested hybrid method in the context of a massive star formation and study the influence of using such a radiation transport method with respect to the previously used flux-limited diffusion approximation.

4.1 Included physics

Our simulations are run with the RAMSES code (Teyssier 2002) which includes a hydrodynamics solver, sink particle algorithm (Bleuler & Teyssier 2014), and radiative transfer with either the flux-limited diffusion module alone (which we call the FLD run) or coupled to RAMSES-RT (the HY run) within our hybrid approach. The opacities were originally used in the low-mass star formation calculations of Vaytet et al. (2013) which include frequency-dependent dust opacities (T < 1500 K, Semenov et al. 2003; Draine 2003). We modify the gray opacities to account for dust sublimation, because its importance for the shielding properties of massive disks has been highlighted in Kuiper et al. (2010c). We model it in the same way as Kuiper et al. (2010a, Eqs. (21) and (22) therein) with a dust-to-gas ratio that decreases with temperature, and a sublimation temperature that increases with the density. The profile of the dust-to-gas mass ratio is given by $$\frac{M_{\text{dust}}}{M_{\text{gas}}}(\rho ,T)={\left(\frac{M_{\text{dust}}}{M_{\text{gas}}}\right)}_0\left(0.5-\frac{1}{\pi }\text{arctan}\left(\frac{T-T_{\text{evap}}(\rho )}{100}\right)\right)$$(15)

where ${\left(\frac{M_{\text{dust}}}{M_{\text{gas}}}\right)}_0$ is the initial dust-to-gas mass ratio, and the evaporation temperature is given by $$T_{\text{evap}}(\rho )=g{\rho }^{\beta }$$(16)

with g = 2000 K cm³ g⁻¹, β =0.0195 (Isella & Natta 2005). At high temperature, when all dust grains are evaporated, the gas opacity is dominant and is taken equal to 0.01 cm² g⁻¹ for comparison purposes with previous studies, such as Krumholz et al. (2009); Kuiper et al. (2014); Rosen et al. (2016); Klassen et al. (2016).

4.2 Setup

We start from initial conditions similar to Rosen et al. (2016): a 150 M_⊙ spherical cloud of radius 0.1 pc in a box of size 0.4 pc to limit boundary effects. The density profile is spherically-symmetric and ρ(r)∝ r^−1.5. The free-fall time is then $${\tau }_{\text{ff}}=\sqrt{\frac{3\pi }{32G\text{ Bar}\rho }}\simeq 42.5 \text{kyr},$$(17)

where G is the gravitational constant and $\rho$ is the mean density computed for a uniform sphere. The density at the border of the cloud is 100 times the density of the ambient medium. The cloud is in solid-body rotation around the x-axis with rotational to gravitational energy of ≈4%, typical of observed cores (Goodman et al. 1993). The initial dust-to-gas mass ratio is ${\left(\frac{M_{\text{dust}}}{M_{\text{gas}}}\right)}_0=0.01$.

The base resolution is level 7 (128³) and the finest resolution is 4096³ (i.e., level 12, five levels of refinement), which gives a physical maximum resolution of 20 AU. In order to limit artificial fragmentation (Truelove et al. 1997), we impose to have at least 12 cells per Jeans length. Sink particles can only form in cells refined to the highest level. Sink creation sites are identified with the clump finder algorithm of Bleuler & Teyssier (2014). The clump finder algorithm marks cells whose density is above a given threshold (3.85 × 10⁻¹⁴ g.cm⁻³ in these calculations). The marked cells are attached to their closest density peak, which form a “peak patch”. We check connectivity between the patches, then the significance of a peak patch is given by the ratio between the peak density and the maximum saddle density lying at a boundary of the peak patch. If the peak-to-saddle ratio is lower than a given value (2 here) the patch is attached to the neighbor patch of highest saddle density. The remaining peak patches are labeled as clumps. Each clump must then meet two conditions to lead to a sink creation: it has to be bound and subvirial. The region around a sink particle is also refined to the highest level. The accretion scheme is based on a density threshold. Consider a cell located within the accretion radius: its accreted mass by the sink is Δm = max(0.25(ρ − ρ_sink) × Δx³, 0), where Δx is the maximum resolution. We merge sinks when they are located in the same accretion volume, while the accretion radius is set to 4Δx ≈ 80 AU. The radius and luminosity of the star mimicked by the sink are computed from the pre-main sequence evolution models of Kuiper & Yorke (2013) and depend on their time-averaged accretion rate and their mass.

4.3 Results

We run one simulation with FLD only (denoted as FLD) and one with FLD+M1 (denoted as HY) until t ≃ 30 kyr ≃ 0.71 τ_ff. As the initial density profile is peaked, a sink particle is expected to form in a few kyr. Both runs lead to the formation of several sink particles, one of which is much more massive than the other sink particles and that we refer to as the main sink or star. A disk and radiative outflows form around the main sink. The criteria for determining the disk and outflows are explained below. We identify a disk on a cell-basis after converting Cartesian coordinates into cylindrical coordinates centered on the mainsink and aligned with the rotational axis, according to the several criteria of Joos et al. (2012):

• The disk is a rotationally-supported structure (i.e., not thermally supported): $\rho v_{\varphi }^2/2>f_{\text{thres}}P$, where v_ϕ is the azimuthal velocity and P is the thermal pressure. The value of f_thres = 2 is chosen, as in Joos et al. (2012). The use of f_thres > 1 leads to a stronger constraint on the identification of cells belonging to the disk;

• In order to avoid large low-density spiral arms, a gas number density threshold is set: n > 1 × 10⁹ cm⁻³;

• The gas is not on the verge of collapsing radially: v_ϕ > f_thresv_r, where v_r is the radial velocity;

• The vertical structure is in hydrostatic equilibrium: v_ϕ > f_thresv_z, where v_z is the vertical velocity.

We define outflows as gas flowing away from the central star at a velocity greater than the escape velocity, which corresponds to $v_r>v_{\text{esc}}=\sqrt{2G{M}_{\star }/r}$, where r is the distance between the sink and thecell and M_⋆ is the sink mass.

The time evolution of the main sink, disk and outflow masses, along with the accretion rate, are displayed in Fig. 9. First, the sink masses are almost equal in both runs before t ≃ 14 kyr and M_⋆ = 5 M_⊙. Even though their evolution differs between 14 and 20 kyr, their values remain similar and the divergence appears only at t ≃ 20 kyr, when M_⋆ = 12 M_⊙ (in the HY run). At that point, the radiative cavities appear in the HY run (see the bottom-right panel of Fig. 9). They appear in the FLD run after the massive star has reached 16   M_⊙. From this time on, the sink mass increases more slowly in the HY run, this can be seen on the accretion rate (top-right panel of Fig. 9). The final stellar mass is M_⋆ = 23.3 M_⊙ in the FLD run and 17.6 M_⊙ in the HY run.

As it is shown on the top-right panel of Fig. 9, the stars experience bursts of accretion separated by ~ 100 yr to a few kyr. These bursts are due to a low-mass companion (sink particle) being accreted by the most massive star. In each run, the main sink experiences accretion rates of M_⊙~10⁻⁴−10⁻²M_⊙ yr⁻¹, which is consistent with previous numerical studies (Klassen et al. 2016) and observations (review by Motte et al. 2018 and references therein). In total, eight companions are formed and accreted in each run. These accretion events contribute to a total of 6.7   M_⊙ in the FLD run and 3.9 M_⊙ in the HY run, hence about 28 and 22% of the final primary star mass, respectively. All sinks appeared after the primary mass was greater than 10   M_⊙ and were accreted in less than 3 kyr (except one, formed at large radius and which does not fall directly onto the primary sink). In each run, the most massive secondary is ~ 2   M_⊙ and gathers most of its mass when orbiting close to the primary, in the disk densest regions. Our merging criterion can lead to overestimating the mass of the primary star and affect the system multiplicity. Assessing the impact of the merging criterion would require a dedicaded study, which is beyond the scope of this paper.

Fig. 9 Time evolution of the main sink mass (top-left panel), accretion rate (top-right), disk mass (bottom-left), and outflow mass (bottom-right).

4.3.1 Disk properties

As shown in the bottom-left panel of Fig. 9, the disks obtained are massive (≈ 17 M_⊙ at the end of the simulation) and similar in mass in both runs. Observational constraints on the disk mass remain sparse (Motte et al. 2018) but the disk mass we obtain is consistent with the previous numerical work of Klassen et al. (2016).

We investigate the disk stability by computing the Toomre parameter Q defined by $$Q=\frac{c_{\text{s}}\kappa }{\pi G\Sigma },$$(18)

where c_s is the sound speed, κ is the epicyclic frequency and is equal to the rotation frequency for a Keplerian disk and Σ is the column density. We recall that the Toomre parameter computes the ratio of the thermal support and differential rotation support over gravitational fragmentation and that the disk is locally unstable if Q < 1. The gas in our simulation is initially in solid-body rotation but the disks formed indeed exhibit rotation curves consistent with Keplerian rotation, as shown in Fig. 10. As a result, the epicyclic frequency κ is equal to Ω, the Keplerian rotational frequency.

To calculate Q, we have taken the column density integrated over the x-axis (perpendicular to the disk). Moreover, the selection given by the criteria presented above gives a disk with a vertical structure. Therefore, we evaluate the Toomre Q parameter in the disk selection, then we average Q over the disk height. For completeness, we have also computed Q with c_s and κ evaluated in the disk midplane and have obtained very similar results.

We also take the radiation into account as an extra-support against fragmentation because the radiative pressure contributes to the sound speed (Mihalas 1984, Eq. (101.22) therein) $$c_{\text{s}}^2={\Gamma }_1\frac{P+P_{\text{r}}}{\rho },$$(19)

where P is the gas pressure, P_r is the radiative pressure; Γ₁ = 5∕3 for a non-radiating fluid (P_r = 0, pure hydrodynamical case), and Γ₁ ≃ 1.43 if P_g = P_r). Therefore, we argue that even for a disk in a strong radiation field and gas-radiation coupling (P_r ~P_g), Q only increases by a factor of ≃1.3 as compared to the pure hydrodynamical case. Figure 11 displays the local Toomre Q value taking the radiative support into accountbut the values of Q without the radiative support lead to the same conclusions.

As shown in Fig. 11, the disks obtained in both runs are Toomre unstable close to the massive star and in the spiral arms. This is consistent with the regular creation of sink particles in those spiral arms. Eight low-mass short-lived companions are generated in both runs. Even though the appearance of sink sparticles is quite resolution-dependent, we limit it with our refinement criterion based on the Jeans length.

The left panel of Fig. 12 shows the main star mass against the disk mass. Both generally increase with time but the disk also undergoes losses of mass as it feeds the main sink particle. Indeed, the main accretion mode in oursimulation is disk accretion, although the accretion bursts are due to the accretion of sink particles recently created in the Toomre unstable spiral arms of the disk. The accretion in our simulation is more stable than what is obtained in the work of Klassen et al. (2016), where the global disk instability leads to an increase of ≃ 10   M_⊙ in a few kyr in their 100 M_⊙ run.

Fig. 10 Radial rotational velocity profile in the disk cells for the HY run at t = 30   kyr and Keplerian profile computed with the main stellar mass. The slope of the velocity profile is consistent with Keplerian rotation.

Fig. 11 Density slices of the disk selection (left panels) and Toomre Q parameter (right panels) in a (2000 AU)2 region centered on the location of the most massive sink particle (left panels), in the FLD run (top) and the HY run (bottom), at t = 30 kyr.

Fig. 12 Primary star mass versus disk mass (left) and outflow mass versus star mass (right), for both FLD and HY runs.

4.3.2 Radiative cavities – outflows

As mentioned in Sect. 4.3: we define outflows as gas flowing away from the central star at a velocity greater than the escape velocity. In the FLD run, radiative cavities appear at t ≃ 22 kyr (bottom-right panel of Fig. 9). They develop earlier and at lower stellar mass (M_⋆ = 12 M_⊙, t ≃ 20 kyr) in the HY run than in the FLD run (M_⋆ = 16 M_⊙, right panel of Fig. 12). In both runs, the cavities grow symmetrically with respect to the disk plane until they reach an extent of ≃ 2000 AU in the FLD run and ≃3000 AU in the HY run at t ≃ 30 kyr (see Fig. 13). The right panels of Fig. 13 display a slice of the density within the outflow selection of cells. The gas velocity is also higher in the HY run, ≈ 25   km s⁻¹, against ≈ 15   km s⁻¹ in the FLD run. As displayed in Fig. 14, gas is pushed away by the radiative force, which locally exceeds gravity. It illustrates the flashlight effect: the radiative force dominates in the poles while the gravity, and hence the accretion, dominates in the disk plane. The consequence of the stronger radiative force is that the outflows in the HY run are able to transport higher density gas than in the FLD run (see right panels of Fig. 13), mainly because it spans a wider angle, particularly in the vicinity of the star (see Fig. 14). Indeed, the outflows displayed in Fig. 13 have masses M_o,HY ≃ 0.6 M_⊙ > M_o,FLD ≃ 0.06 M_⊙, as displayedon the bottom-right panel of Fig. 9. During almost all the simulations the outflows in the HY run are more massive than in the FLD run. In addition, the temporal evolution at t ≃ 30 kyr seems to show that this mass is still going to increase in the HY run but not in the FLD run. We finally note that the peak in the outflow mass at t ≃ 20 kyr is due to the launching of the outflows in a high-density medium close to the star, which therefore gives a higher outflow mass: the low-density cavity has not formed yet.

4.3.3 Rayleigh–Taylor instabilities

No Rayleigh–Taylor instabilities appear in the aforementioned runs (FLD and HY). They have been shown to contribute significantly to the star-disk system evolution and to the final mass of the star in several studies (Krumholz et al. 2009; Rosen et al. 2016, 2019). Discussions about the presence of these instabilities in massive star formation simulations lean on arguments of numerical resolution, since the smaller-scale modes are the most unstable (Jacquet & Krumholz 2011). Here we try to tackle this problem by cranking up the resolution to see if we get any.

We conduct a run whose spatial resolution permits one to resolve the seeds of radiative Rayleigh–Taylor instabilities. We call this run HY-RTi, based on the HY run restarted at the time when radiative cavities appear. We rely on the AMR framework to resolve the radiative cavities interfaces with a refinement strategy based on the gradient of the stellar radiation $$\frac{\nabla E_{\text{M1}}\Delta x}{E_{\text{M1}}}<10\%,$$(20)

where E_M1 is the M1 module radiative energy and Δx the cell width. This means that if the radiative energy of the M1 module changes by more than 10% between two adjacent cells, these cells are flagged for refinement. We add a second and a third conditions to flag a cell: E_M1 > E_thres = 3 × 10⁻¹² erg cm⁻³, and $\Vert x\Vert >1500$ AU (over and under the sink-disk plane). The last condition is applied once the cavities are developed beyond this height. The second and third criteria are necessary to not over-refining other regions than the interfaces of the cavities, which would explode the cost of the simulation.

The left panel of Fig. 15 shows the radiative cavities in the HY-RTi run with AMR level contours overplotted and the right panel of Fig. 15 displays the M1 radiative energy with respect to the radius. The zone within the contours “11” is at the AMR level 12, which is the highest level. Contours show that the radiative energy drop around R = 5000 AU is resolved to the AMR level 12, which is the highest level. As a result, the finest resolution is put on the cavity interfaces. Despite this refinement strategy, no radiative Rayleigh–Taylor instability has developed in any of our simulations. We explain below why this result is not numerical but physical.

We compare the typical advection time of the flow τ_adv and the growth time of the instability τ_instab; the condition forthe instability to develop is τ_adv > 3τ_instab (Foglizzo et al. 2006). First, the flow in the bubble has supersonic speeds and forms a shock of thickness H ≃ 300 AU (measured on the density profile) when it encounters the accretion flow. In the shock frame (whose velocity is 2 km s⁻¹), we measure a gas velocity of 10 km s⁻¹. Hence, the advection timescale of the gas in the shock is τ_adv ≃ 0.1 kyr.

We now compute the growth rate of the Rayleigh–Taylor instability for the shortest perturbations we can capture (of spatial scale λ = Δx_min = 20 AU, which is the fastest growing mode) using the equation (80) from Jacquet & Krumholz (2011). We obtain a growth rate of ω ≃ 4.5 kyr⁻¹, hence a growth timescale of τ_instab ≃ 0.2 kyr. This is longer than the advection timescale τ_adv, so the gas isadvected before the instability develops. Furthermore, the calculation in Jacquet & Krumholz (2011) is based on the adiabatic approximation which is valid when the cavity edge temperature is taken to be equal to the dust sublimation temperature (~ 1100 K). Numerically, we get a temperature of a few ~100 K at the cavity edge, thus the adiabatic approximation breaks down in our simulation. Physically, the cavity edge is mainly heated by stellar radiation, which is geometrically diluted in the optically-thin cavity. Therefore, it can be shown that the cavity edge should have a temperature of a few 100 K at a distance of ~ 3000 AU from a ~ 10⁵ L_⊙ source. Moreover, the cavity interior is optically-thin and thus is not adiabatic, as mentioned in Jacquet & Krumholz (2011). If compressed, the gas radiates away its energy instead of heating as an optically-thick gas (adiabatic) would.

For these reasons, we go one step further and relax the adiabatic approximation in the cavity interior. Hence, the entropy within the cavity cannot account for radiation. We compute the total entropy (gas plus radiation) as a function of the coupling between gas and radiation via the local optical depth τ $$s_{\text{tot}}=\frac{k_{\text{B}}}{m(\gamma -1)}\ln \left(P_{\text{g}}{\rho }^{-\gamma }\right)+\min (\tau ,1)\frac{4P_{\text{r}}}{\rho T},$$(21)

where k_B is Boltzmann’s constant, m is the molecular hydrogen mass, P_g is the gas pressure, and P_r is the radiation pressure. The maximum growth rate is given by the Brunt-Väisäla (or buoyancy) frequency, which is the oscillation frequency of a fluid particle in a stratified medium $$\omega =\sqrt{\frac{\gamma -1}{\gamma }{g}_{\text{eff}}\nabla S},$$(22)

where S is the total entropy s_tot normalized by k_B∕m and g_eff is the effective gravity g_eff = g − κF∕c, with κF∕c the radiative acceleration. We compute ω in our simulation at the bubble edge and get ω ≲ 10 kyr⁻¹, which gives τ_instab ≳ 0.1 kyr ≃ τ_adv. Therefore, no Rayleigh–Taylor instability should develop in our simulation.

Fig. 13 Left panels: density slices perpendicular to the disk in a (10000 AU)2 region. Right panels: density in outflow selections in a (8000 AU)2 region. Top panels: FLD run; bottom panels: HY run. t = 30 kyr. Figures are centered on the location of the most massive sink particle.

Fig. 14 Radiative force to gravitational force normalized ratio in the FLD run (left panel) and HY run (right panel) in a (10 000 AU)2 region perpendicular to the disk. t = 30 kyr. Figures are centered on the location of the most massive sink particle. Regions of outflows (see Fig. 13) are dominated by the radiative force.

Fig. 15 HY-RTi run at t~0.7τff. Left panel: density slice perpendicular to the disk in a (10 000 AU)2 region. Right panel: scatter plot of the M1 radiative energy against the radius. Contours show the AMR level. The cavity edges are zones of primary absorption for the stellar radiation and are resolved to the highest level (12).

5 Conclusions and discussion

We have implemented a new hybrid radiative transfer method in the AMR code RAMSES based on the flux-limited diffusion (FLD) module (Commerçon et al. 2011b) and M1 module in RAMSES-RT (Rosdahl et al. 2013), in order to treat accurately both the stellar irradiation and the diffuse component around a massive protostar. Our hybrid approach takes advantage of the M1 module fully tested in the optically-thin regime and of the FLD module in the optically-thick limit. Moreover, in contrast to the consideration of local photons inherent to the FLD approach, our method keeps the frequency information of the stellar photons propagated with RAMSES-RT, which leads to an improvement in the treatment of coupling with the dust-gas mixture via the temperature and the radiative force.

We tested this improvement in pure radiative transfer tests of a star irradiating a static disk structure (Pascucci et al. 2004, Pinte et al. 2009). Our results show that the hybrid method is very accurate is the optically-thin regime (≈ 2% maximal error, ≈62% with the FLD method alone), and more accurate than the FLD method in the optically-thick regime (≈ 25% instead of ≈ 36% with the FLD method). It is also capable of capturing partially the self-shielding in the disk mid-plane, because it is shielded from stellar radiation by the disk inner region. Therefore, the hybrid method is suited to determine accurately the gas temperature structure in different regimes of optical thickness. In addition, because stellar photons are treated apart from the photons emitted by the dusty disk, the associated radiative force is computed more consistently with the hybrid method and its value is about ~ 100 times greater than with pure FLD. This shows the need for such an hybrid method, beyond the diffusion approaches.

After testing our hybrid approach in pure radiative transfer cases, we have applied it to a radiation-hydrodynamical problem: the collapse of a massive prestellar core. Both runs lead to the formation of a massive star. The multi-dimensionality of our simulations leads to accretion via the flashlight effect: accretion occurs through a disk while radiation escapes via the poles (Yorke & Sonnhalter 2002). Low-mass sink particles are created in the Toomre unstable spiral arms of the disk, they move together with the fluid and are rapidly accreted onto the central massive star. At the end of the simulation (t ≃ 30 kyr ≃ 0.71τ_ff), the star mass is 23.3   M_⊙ in the FLD run and 17.6   M_⊙ in the HY run, showing no signs of decrease in the accretion. This difference is explained by the direct radiative pressure onto the disk, which lowers the accretion rate. When the star reaches 12 M_⊙ (HY run) or 16 M_⊙ (FLD run), radiative polar cavities develop because of the stellar radiative pressure. Radiative outflows in the HY run are ~ 50% more extended than in the FLD run.

Our method also contains a few assumptions and limitations we shall discuss. As shown above, the hybrid method is accurate within ≈25−65% in the optically-thick limit. In the prestellar core collapse problem, the accretion disk around the protostellar source is very optically-thick (≳ 10³) and the photon mean free path is barely resolved in AMR codes with current computational facilities. The disk midplane temperature is therefore affected by this error: it is generally overestimated, which increases the Jeans length and therefore stability. Yet, the hybrid method clearly performs better than the pure FLD approach in predicting the mid-plane temperature. Also, the temperature in the optically-thin cavities is computed more accurately with the hybrid method.

In this study, we have focused on a gray approach for both the direct and the diffuse radiations. However, the gray radiative transfer is not inherent to our model and multifrequency methods have been implemented in RAMSES-RT (Rosdahl et al. 2013) and in the FLD module (González et al. 2015). We have chosen to work with the gray methods to save computational time and memory. The multifrequency version of our hybrid method is beyond the scope of this paper. We have also focusedon the irradiation by one source with the M1 method. However, the several sinks produced in these simulations were rapidly accreted by the primary one, so this does not change our conclusions. Moreover, the irradiation can be generalized to several sources via one photon group per frequency band or per source, and we leave this to further work. Finally, ionization processes were not taken into account here. They can be relevant for massive starformation, but the opening of the ionized cavities and the disk photo-evaporation has been shown to occur toward the end of our simulation (Kuiper & Hosokawa 2018). These physics will be included in further works.

RAMSES includes magneto-hydrodynamics (Fromang et al. 2006), and the interplay between the radiation (modeled with FLD) and the magnetic field has been proven: magnetic braking enhances the accretion speed and hence the radiative shock energy release (Commerçon et al. 2011a). This energy heats the disk and prevents further fragmentation. The treatment of direct stellar irradiation within our hybrid method can modify this interplay because the disk self-shielding will be captured better than with the FLD method (if the mean-free path is resolved, see Sect. 3.2), and also via the launching of magnetic and radiative outflows. Thus, the hybrid method offers new perspectives for the radiation-magneto-hydrodynamics simulations of massive star formation.

Acknowledgements

Part of this work was supported by the CNRS “Programme National de Physique Stellaire” (PNPS). The numerical simulations presented here were run on the CEA machines Irfucoast and Alfvén. The visualisation of RAMSES data has been executed with the OSYRIS python package. We thank the referee for constructivecomments. R.M.R. acknowledges C. Pinte for MCFOST data, and J. Ramsey, R. Kuiper, and S. Fromang for useful discussions. R.M.R. finally thanks T. Foglizzo for his fruitful insight on Rayleigh–Taylor instabilities.

Appendix A Temperature structure with isotropic scattering

Figure 1 shows that the extinction opacity is dominated by scattering at high frequency, that is, where stellar irradiation dominates. Hence, we examine how the gray opacities are modified when including isotropic scattering in the FLD and M1 equations.

It can be shown that taking the scattering into account does not modify the first moment of the RT equation (the conservationof the radiative energy), because the source term involves the radiative energy and its redistribution is the same with and withoutscattering.

However, the first moment of the RT equation (see Eq. (5)) is modified because the coupling between the gas and the radiative flux is enhanced by the scattering. Therefore, the opacity in this equation is the extinction (absorption+scattering) opacity instead of the absorption opacity.

Similarly, the gray version of the flux evolution equation in the M1 module introduces the Planck mean opacity as defined in Eq. (3) but with the extinction opacities.

In fact, the inclusion of isotropic scattering has no explicit impact on the radiative energy repartition but it increases the absorption and redistribution of the radiative flux. To test the behavior of the hybrid approach with isotropic scattering we run the most optically-thick case of the setup from Pascucci et al. (2004), with τ = 100, following the configuration and numerical parameters of Sect. 3.1.

The left panel of Fig. A.1 shows that the temperature in the mid-plane of the disk is well reproduced by the hybrid approach, as compared to the result obtained with MCFOST. The maximal error is ≈ 20% with the hybrid against ≈40% with the FLD method alone. The right panel of Fig. A.1 emphasizes that most of the shielding effect in the disk mid-plane is captured by the hybrid approach, unlike for the FLD method. We compare both panels of Fig. A.1 with the left panel of Fig. 4 and right panel of Fig. 5, respectively, as they show the result of the same setup without scattering. We observe that the temperature structure is not much influenced by the treatment of scattering. First, as shown in Fig. 1, scattering only dominates at high-frequencies and therefore mainly at the first interaction of stellar photons with the medium, which corresponds to the disk inner edge. After this interaction, photons are reemitted at the local temperature, where the absorption opacity dominates. This explains why the temperature at the disk inner edge is ~ 400 K with isotropic scattering against ~350 K without scattering. Second, since opacities are frequency-averaged in our hybrid approach, taking isotropic scattering into account does not impact much the Planck and Rosseland mean opacities. Therefore, we do not consider scattering in the collapse calculations in Sect. 4.

Fig. A.1 Gas temperature profiles, following the test of Pascucci et al. (2004). We compare the gas temperature computed using MCFOST, the hybrid method (M1+FLD) and the FLD method alone in RAMSES, with isotropic scattering. T⋆ = 5800 K and the integrated optical depth in the disk mid-plane is τ = 100. Top: radial profile in the disk mid-plane. Bottom: vertical profile at a disk radius of 20 AU.

Appendix B Performance test

In this section, we compare the performance of the FLD and the hybrid methods. First, we run the test from Pascucci et al. (2004), with τ = 100 and T_⋆ = 5800 K to probe thescaling of each method. For this test, we choose a grid with two levels of refinement: 8 and 9, which leads to ~ 2 × 10⁷ cells. Figure B.1 shows the strong scaling results from 2 to 32 cores. We can see that the scaling properties of the FLD and the hybrid methods are very similar and close to the theoretical line. The departure from the theoretical line (speedup of ≃ 11 instead of 16 for 32 cores) is likely due to the high number of global communications which occur in the FLD conjugate gradient algorithm.

Fig. B.1 Strong scaling result for 2–32 cores in the test of Pascucci et al. (2004), with τ = 100, T⋆ = 5800 K. The ideal theoretical speedup is represented by the black line. We compare the strong scaling between the hybrid method (M1+FLD, blue squares) and the FLD method alone (orange circles). We normalize the speedup by that obtained with two cores.

The FLD implementation without hydrodynamics is implicit and therefore is not restricted by the CFL condition in our pure radiative transfer tests, in constrast to the M1 part of our hybrid method. Therefore, we do not compare the computational time in those tests. However, we look at the total CPU time in the collapse calculations presented in Sect. 4. The HY run took ≈ 5100 CPU hours, against ≈ 3900 for the FLD run, which consists in an additional time of about ≈30%. The time step in the HY run is first constrained by the M1 CFL condition when the primary sink forms, which is responsible for this difference of computational time: more steps were needed to reach the same physical time. As mentioned in Sect. 2, the FLD modifies the hydrodynamical CFL time step : the sound speed accounts for both thermal and radiative pressures. It decreases as the radiative pressure (hence energy) increases, while the M1 time step is fixed. As the central mass gains mass, its temperature and luminosity generally rise and so does the radiative energy. Therefore, the FLD time step decreases in both runs, but is still greater than the M1 time step in the HY run, at first. Then when the outflowsare launched, both runs are limited by the FLD time step, because both the radiative energy has become significant and the density is very low in the outflow (hence the modified sound speed increases). From this time on, thetime step is comparable in both runs. However, the number of iterations in the conjugate gradient is ~ 10% smaller in the HY run than in the FLD run, and so is the elapsed time per time step.

References

All Tables

All Figures

Fig. 1 Frequency-dependent opacities and blackbody spectra for Tdisk = 300 K and T⋆ = 5800 K. Opacities are absorption (blue pluses), scattering (red crosses) and extinction (black dots) coefficients for the dust-and-gas mixture used in the Pasccuci test. The table contains 61 frequency bins and dataare taken from Draine & Lee (1984). Apart from the broad opacity features at about 10 and 20 μm, which correspond to Si-O vibrational transitions, the opacity generally increases with the photon frequency. The opacity at stellar-like radiation frequencies is generally greater than at disk-like radiation frequencies.

In the text

	Fig. 2 Planck’s (blue dashed curve) and Rosseland’s (orange dot-dashed curve) mean opacities, as a function of temperaturein the Pascucci setup.
In the text

	Fig. 3 Radial gas temperature profiles in the mid-plane of the disk following the test of Pascucci et al. (2004) for τ = 0.1. We comparethe gas temperature computed using MCFOST (black dotted-line) and RADMC-3D (red dashed-line), the hybrid method (M1+FLD, blue dots) and the FLD method alone (orange dots) in RAMSES. Left: central star temperature T⋆,1 = 5800 K; right: T⋆,2 = 15 000 K.
In the text

	Fig. 4 Same as Fig. 3, but for τ = 100.
In the text

	Fig. 5 Vertical gas temperature profiles at a cylindrical radius of 20 AU, following the test of Pascucci et al. (2004). We compare the gas temperature computedusing  MCFOST, the hybrid method (M1+FLD) and FLD alone in RAMSES for T⋆,1, τ = 0.1 (left) and τ = 100 (right).
In the text

	Fig. 6 1000 AU disk edge-on slices of the radiative acceleration, following the test of Pascucci et al. (2004) obtained with RAMSES after stationarity is reached. Left: FLD run; right: hybrid run. Star and disk parameters: τ = 100 and T⋆,1. The hybrid radiative acceleration is about 100 times greater than the FLD one.
In the text

	Fig. 7 AMR level needed to resolve the mean free path (mfp) of photons (red dashed-line) and effective AMR level (black dots) in the disk midplane following the test of Pinte et al. (2009).
In the text

Fig. 8 Left: radial gas temperature profile in the mid-plane of the disk following the test of Pinte et al. (2009). Right: vertical gas temperature profile at a cylindrical radius of 0.2 AU in the disk. We compare the gas temperature computed using  RADMC-3D, the hybrid method (M1+FLD) and the FLD method alone in RAMSES. The integrated optical depth in the disk mid-plane is τ810 nm = 103 and the stellar temperature is T⋆ = 4000 K.

In the text

	Fig. 9 Time evolution of the main sink mass (top-left panel), accretion rate (top-right), disk mass (bottom-left), and outflow mass (bottom-right).
In the text

	Fig. 10 Radial rotational velocity profile in the disk cells for the HY run at t = 30   kyr and Keplerian profile computed with the main stellar mass. The slope of the velocity profile is consistent with Keplerian rotation.
In the text

	Fig. 11 Density slices of the disk selection (left panels) and Toomre Q parameter (right panels) in a (2000 AU)2 region centered on the location of the most massive sink particle (left panels), in the FLD run (top) and the HY run (bottom), at t = 30 kyr.
In the text

	Fig. 12 Primary star mass versus disk mass (left) and outflow mass versus star mass (right), for both FLD and HY runs.
In the text

	Fig. 13 Left panels: density slices perpendicular to the disk in a (10000 AU)2 region. Right panels: density in outflow selections in a (8000 AU)2 region. Top panels: FLD run; bottom panels: HY run. t = 30 kyr. Figures are centered on the location of the most massive sink particle.
In the text

	Fig. 14 Radiative force to gravitational force normalized ratio in the FLD run (left panel) and HY run (right panel) in a (10 000 AU)2 region perpendicular to the disk. t = 30 kyr. Figures are centered on the location of the most massive sink particle. Regions of outflows (see Fig. 13) are dominated by the radiative force.
In the text

	Fig. 15 HY-RTi run at t~0.7τff. Left panel: density slice perpendicular to the disk in a (10 000 AU)2 region. Right panel: scatter plot of the M1 radiative energy against the radius. Contours show the AMR level. The cavity edges are zones of primary absorption for the stellar radiation and are resolved to the highest level (12).
In the text

	Fig. A.1 Gas temperature profiles, following the test of Pascucci et al. (2004). We compare the gas temperature computed using MCFOST, the hybrid method (M1+FLD) and the FLD method alone in RAMSES, with isotropic scattering. T⋆ = 5800 K and the integrated optical depth in the disk mid-plane is τ = 100. Top: radial profile in the disk mid-plane. Bottom: vertical profile at a disk radius of 20 AU.
In the text

	Fig. B.1 Strong scaling result for 2–32 cores in the test of Pascucci et al. (2004), with τ = 100, T⋆ = 5800 K. The ideal theoretical speedup is represented by the black line. We compare the strong scaling between the hybrid method (M1+FLD, blue squares) and the FLD method alone (orange circles). We normalize the speedup by that obtained with two cores.
In the text