Issue 
A&A
Volume 635, March 2020



Article Number  A42  
Number of page(s)  17  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201936605  
Published online  06 March 2020 
A new hybrid radiative transfer method for massive star formation
^{1}
AIM, CEA, CNRS, Univ. ParisSaclay, Univ. Paris Diderot,
Sorbonne Paris Cité,
91191
GifsurYvette,
France
email: raphael.mignonrisse@cea.fr
^{2}
Centre de Recherche Astrophysique de Lyon UMR5574, ENS de Lyon, Univ. Lyon1, CNRS, Université de Lyon,
69007
Lyon, France
^{3}
Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574,
69230
SaintGenisLaval, France
Received:
28
August
2019
Accepted:
28
January
2020
Context. Frequencydependent and hybrid approaches for the treatment of stellar irradiation are of primary importance in numerical simulations of massive star formation.
Aims. We seek to compare outflow and accretion mechanisms in star formation simulations. We investigate the accuracy of a hybrid radiative transfer method using the gray M1 closure relation for protostellar irradiation and gray fluxlimited diffusion (FLD) for photons emitted everywhere else.
Methods. We have coupled the FLD module of the adaptivemesh refinement code RAMSES with RAMSESRT, which is based on the M1 closure relation and the reduced speedoflightapproximation. Our hybrid (M1+FLD) method takes an average opacity at the stellar temperature for the M1 module, instead of the local environmental radiation field. Due to their construction, the opacities are consistent with the photon origin. We have tested this approach in radiative transfer tests of disks irradiated by a star for three levels of optical thickness and compared the temperature structure with the radiative transfer codes RADMC3D and MCFOST. We applied it to a radiationhydrodynamical simulation of massive star formation.
Results. Our tests validate our hybrid approach for determining the temperature structure of an irradiated disk in the opticallythin (2% maximal error) and moderately opticallythick (error smaller than 25%) regimes. The most opticallythick test shows the limitation of our hybrid approach with a maximal error of 65% in the disk midplane against 94% with the FLD method. The opticallythick setups highlight the ability of the hybrid method to partially capture the selfshielding in the disk while the FLD alone cannot. The radiative acceleration is ≈100 times greater with the hybrid method than with the FLD. The hybrid method consistently leads to about + 50% more extended and widerangle radiative outflows in the massive star formation simulation. We obtain a 17.6 M_{⊙} star at t ≃ 0.7τ_{ff}, while the accretion phase is still ongoing, with a mean accretion rate of ≃7 × 10^{−4} M_{⊙} yr^{−1}. Finally, despite the use of refinement to resolve the radiative cavities, no Rayleigh–Taylor instability appears in our simulations, and we justify their absence by physical arguments based on the entropy gradient.
Key words: stars: formation / stars: massive / stars: protostars / radiative transfer / hydrodynamics / methods: numerical
© R. MignonRisse et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Massive stars shape the dynamical and chemical evolution of galaxies because of their powerful feedback in radiation, winds, explosions in supernova, and metalenrichment. However, their formation remains a longstanding problem. Observationally, massive stars are embedded in dense clouds; they form on timescales much shorter than their lowmass 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 initialmass function (IMF) is builtup naturally. On the other hand, the turbulent core model (McKee & Tan 2003) is an extension of the lowmass 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^{−4}−10^{−3} M_{⊙} yr^{−1}; 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 wellsuited in the opticallythin regime. Current and past studies of massive star formation have mainly focused on its radiative aspects. In a 1D sphericallysymmetric 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_{⊙}. Twodimensional 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).
MonteCarlo 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 fluxlimited diffusion (FLD) and the M1 methods for radiationhydrodynamics (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 opticallythick regime while the use of a fluxlimiter corrects the propagation speed in the opticallythin regime (Levermore & Pomraning 1981). In addition, the FLD method does not permit to capture shadows behind very opticallythick 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 opticallythin limit. Recent approaches treat stellar irradiation in a more consistent way with raytracing (Kuiper et al. 2010b; Kim et al. 2017), longcharacteristics (Rosen et al. 2017), MonteCarlo radiative transfer (Haworth & Harries 2012; Harries et al. 2017, who took advantage of the independency between photon packets to parallelize efficiently the MonteCarlo 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).
Multidimensional 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 freefall times. Most of these works put lowerlimit 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 radiationpressuredominated cavities that feed the stardisk 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 frequencydependent 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 diskaccretion is sufficient to reach masses consistent with the massive stars observed. Our hybrid method, implemented in the Cartesian adaptivemesh 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 fluxlimited 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 fluxlimited diffusion method for the dust emission.
2.1 Coupling fluxlimited 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) (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 dusttogas 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 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 (2)
where E_{r} is the frequencyintegrated radiative energy, λ is the fluxlimiter and is built to recover the right propagation speed in opticallythin 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 (3)
where the temperature derivatives appear from the chain rule applied to ∇B. The a T^{4} 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 (5)
where is the rate of radiative energy injected from stellar sources. F_{r} and are the frequencyintegrated radiative flux and pressure respectively.
One main asset is that the directionality of the photons beam is wellretained. 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 opticallythick medium and M1 for irradiation. Both FLD and M1 methods described above can involve several groups of photons or only one with frequencyaveraged 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 UVlike energy and relatively high opacities, and photons emitted by the dust, which have a IRlike 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 (6)
where 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_{v}T 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 (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 fluxaveraged opacity – which means that momentum is transferred each time a photon is absorbed – so the radiative acceleration is given by (8)
On the other hand, the gray radiative acceleration in the FLD approximation is given by (9)
Taking the asymptotic values of λ into account (see Levermore 1984), we get in the limits of low optical depth and 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 adaptivemesh 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 volumeaveraged over the cell) and a secondorder Godunov method is used to evolve hydrodynamical variables. This code includes theFLD (Commerçon et al. 2011b, C11 hereafter) and the M1 method within RAMSESRT (Rosdahl et al. 2013^{1}, R13 hereafter), both coupled to the hydrodynamics. For the FLD method, a timesplitting approach is performed. A predictorcorrector 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 secondorder Godunov scheme of RAMSES and fluxes are estimated with an approximate Riemann solver (LaxFriedrich, HLL, HLLD, etc.). The diffusion and radiationmatter coupling are handled in the implicit part of the timesplitting scheme. The diffusion part of the FLD solver is secondorder 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 CourantFriedrichLewy (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 radiationmatter 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 firstorder Godunov scheme, thus it obeys the CFL condition. The trick used here is the reduced speedoflight 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 MonteCarlo radiative transfer codes. We explore three levels of optical thickness integrated along the disk midplane: τ = 0.1, τ = 100 and τ = 10^{3}. The parameters and results are summarized in Table 1. We refer the reader to Appendix B for performance tests.
3.1 Opticallythin and moderately opticallythick 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 opticallythin and moderately opticallythick regimes. In particular, we compare our results once the temperature structure is converged with respect to time with those obtained with MonteCarlo RT codes such as RADMC3D (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 (10)
where ρ_{0} is the density normalization and is linked to the only freeparameter, the integrated opticaldepth throughout the midplane of the disk, . The two functions f_{1} and f_{2} are given by (11)
where the flaring function is (13)
In this setup, r_{d} = r_{out}∕2 = 500 AU and z_{d} = r_{out}∕8 = 125 AU are the scaleradius and the scaleheight. 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 opticallythin and moderately opticallythick regimes, respectively. The dusttogas 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^{−3}. Frequencydependent 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^{−23} g cm^{−3}. The same density floor is applied between the star and the disk edge to mimic the vacuum that RADMC3D 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 freepath 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^{−4} (see Ramsey & Dullemond 2015).
Results from pure radiative transfer tests.
Fig. 1 Frequencydependent opacities and blackbody spectra for T_{disk} = 300 K and T_{⋆} = 5800 K. Opacities are absorption (blue pluses), scattering (red crosses) and extinction (black dots) coefficients for the dustandgas 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 SiO vibrational transitions, the opacity generally increases with the photon frequency. The opacity at stellarlike radiation frequencies is generally greater than at disklike radiation frequencies. 
Fig. 2 Planck’s (blue dashed curve) and Rosseland’s (orange dotdashed 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 RADMC3D 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 midplane against the xaxis for the most opticallythin 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 midplane, 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 RADMC3D) 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 wellsuited for opticallythin media. In the hybrid method, since there is not much absorption because of the low opticaldepth, the M1 module is mainly at work and is adapted to opticallythin 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 midplane for the moderately opticallythick 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 opticallythin 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 frequencydependence of RADMC3D 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 opticallythick disk produces selfshielding in the midplane 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 RADMC3D 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 opticallythin case, τ = 0.1. No selfshielding 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 midplane 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 opticallythin and moderately opticallythick regime. The hybrid method is very accurate in the opticallythin regime (less than ≈ 2%). In the moderately opticallythick regime, the hybrid method gives a nonnegligible error (up to ≈ 31% for a 15 000 K star) in the transition between opticallythin 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 selfshielding in the disk midplane while the FLD approximation does not.
Fig. 3 Radial gas temperature profiles in the midplane of the disk following the test of Pascucci et al. (2004) for τ = 0.1. We comparethe gas temperature computed using MCFOST (black dottedline) and RADMC3D (red dashedline), 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. 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 opticallythick 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 opticallythin limit it is proportional to the radiative energy and in the opticallythick 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 opticallythin 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 opticallythin 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 Opticallythick 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^{3}. The flared disk density profile ρ(r, z) is analytically given by (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^{−3}.
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 AMRgrid codes than for cylindricalgrid 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 opticallythin 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 logarithmicallyscaled 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^{−4} 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 RADMC3D.
The left panel of Fig. 8 plots the radial temperature profile in the disk midplane obtained with RAMSES with the FLD and hybrid methods versus RADMC3D. 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 RADMC3D 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 RADMC3D but with selfshielding 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 opticallythick 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 frequencydependent 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 edgeon 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 dashedline) 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 midplane 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 RADMC3D, the hybrid method (M1+FLD) and the FLD method alone in RAMSES. The integrated optical depth in the disk midplane is τ_{810 nm} = 10^{3} 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 fluxlimited 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 fluxlimited diffusion module alone (which we call the FLD run) or coupled to RAMSESRT (the HY run) within our hybrid approach. The opacities were originally used in the lowmass star formation calculations of Vaytet et al. (2013) which include frequencydependent 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 dusttogas ratio that decreases with temperature, and a sublimation temperature that increases with the density. The profile of the dusttogas mass ratio is given by (15)
where is the initial dusttogas mass ratio, and the evaporation temperature is given by (16)
with g = 2000 K cm^{3} g^{−1}, β =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^{2} g^{−1} 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 sphericallysymmetric and ρ(r)∝ r^{−1.5}. The freefall time is then (17)
where G is the gravitational constant and 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 solidbody rotation around the xaxis with rotational to gravitational energy of ≈4%, typical of observed cores (Goodman et al. 1993). The initial dusttogas mass ratio is .
The base resolution is level 7 (128^{3}) and the finest resolution is 4096^{3} (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^{−14} g.cm^{−3} 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 peaktosaddle 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^{3}, 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 premain sequence evolution models of Kuiper & Yorke (2013) and depend on their timeaveraged 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 cellbasis 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 rotationallysupported structure (i.e., not thermally supported): , 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 lowdensity spiral arms, a gas number density threshold is set: n > 1 × 10^{9} cm^{−3};

The gas is not on the verge of collapsing radially: v_{ϕ} > f_{thres}v_{r}, where v_{r} is the radial velocity;

The vertical structure is in hydrostatic equilibrium: v_{ϕ} > f_{thres}v_{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 , 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 bottomright 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 (topright 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 topright panel of Fig. 9, the stars experience bursts of accretion separated by ~ 100 yr to a few kyr. These bursts are due to a lowmass companion (sink particle) being accreted by the most massive star. In each run, the main sink experiences accretion rates of M_{⊙}~10^{−4}−10^{−2}M_{⊙} yr^{−1}, 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 (topleft panel), accretion rate (topright), disk mass (bottomleft), and outflow mass (bottomright). 
4.3.1 Disk properties
As shown in the bottomleft 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 (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 solidbody 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 xaxis (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 extrasupport against fragmentation because the radiative pressure contributes to the sound speed (Mihalas 1984, Eq. (101.22) therein) (19)
where P is the gas pressure, P_{r} is the radiative pressure; Γ_{1} = 5∕3 for a nonradiating fluid (P_{r} = 0, pure hydrodynamical case), and Γ_{1} ≃ 1.43 if P_{g} = P_{r}). Therefore, we argue that even for a disk in a strong radiation field and gasradiation 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 lowmass shortlived companions are generated in both runs. Even though the appearance of sink sparticles is quite resolutiondependent, 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 (bottomright 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^{−1}, against ≈ 15 km s^{−1} 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 bottomright 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 highdensity medium close to the star, which therefore gives a higher outflow mass: the lowdensity 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 stardisk 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 smallerscale 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 HYRTi, 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 (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^{−12} erg cm^{−3}, and AU (over and under the sinkdisk plane). The last condition is applied once the cavities are developed beyond this height. The second and third criteria are necessary to not overrefining 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 HYRTi 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^{−1}), we measure a gas velocity of 10 km s^{−1}. 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^{−1}, 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 opticallythin 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^{5} L_{⊙} source. Moreover, the cavity interior is opticallythin and thus is not adiabatic, as mentioned in Jacquet & Krumholz (2011). If compressed, the gas radiates away its energy instead of heating as an opticallythick 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 τ (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 BruntVäisäla (or buoyancy) frequency, which is the oscillation frequency of a fluid particle in a stratified medium (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^{−1}, 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 HYRTi 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 fluxlimited diffusion (FLD) module (Commerçon et al. 2011b) and M1 module in RAMSESRT (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 opticallythin regime and of the FLD module in the opticallythick 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 RAMSESRT, which leads to an improvement in the treatment of coupling with the dustgas 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 opticallythin regime (≈ 2% maximal error, ≈62% with the FLD method alone), and more accurate than the FLD method in the opticallythick regime (≈ 25% instead of ≈ 36% with the FLD method). It is also capable of capturing partially the selfshielding in the disk midplane, 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 radiationhydrodynamical problem: the collapse of a massive prestellar core. Both runs lead to the formation of a massive star. The multidimensionality of our simulations leads to accretion via the flashlight effect: accretion occurs through a disk while radiation escapes via the poles (Yorke & Sonnhalter 2002). Lowmass 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 opticallythick limit. In the prestellar core collapse problem, the accretion disk around the protostellar source is very opticallythick (≳ 10^{3}) 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 midplane temperature. Also, the temperature in the opticallythin 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 RAMSESRT (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 photoevaporation has been shown to occur toward the end of our simulation (Kuiper & Hosokawa 2018). These physics will be included in further works.
RAMSES includes magnetohydrodynamics (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 selfshielding will be captured better than with the FLD method (if the meanfree 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 radiationmagnetohydrodynamics 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 opticallythick 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 midplane 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 midplane 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 highfrequencies 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 frequencyaveraged 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 midplane is τ = 100. Top: radial profile in the disk midplane. 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^{7} 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
 Aubert, D., & Teyssier, R. 2008, MNRAS, 387, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Beuther, H., Churchwell, E. B., McKee, C. F., & Tan, J. C. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 165 [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015 [NASA ADS] [CrossRef] [Google Scholar]
 Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011b, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012, A&A, 545, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Draine, B. T. 2003, ApJ, 598, 1026 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [Google Scholar]
 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
 Foglizzo, T., Scheck, L., & Janka, H. T. 2006, ApJ, 652, 1436 [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gnedin, N. Y., & Abel, T. 2001, New Astron., 6, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Goddi, C., Ginsburg, A., Maud, L., Zhang, Q., & Zapata, L. 2018, ArXiv eprints [arXiv:1805.05364] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528 [NASA ADS] [CrossRef] [Google Scholar]
 Harries, T. J., Douglas, T. A., & Ali, A. 2017, MNRAS, 471, 4111 [NASA ADS] [CrossRef] [Google Scholar]
 Haworth, T. J., & Harries, T. J. 2012, MNRAS, 420, 562 [Google Scholar]
 Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F., & Grudić, M. Y. 2019, MNRAS, 483, 4187 [NASA ADS] [CrossRef] [Google Scholar]
 Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E., & Krumholz, M. R. 2011, ApJ, 730, 116 [Google Scholar]
 Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kannan, R., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 485, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J.G., Kim, W.T., Ostriker, E. C., & Skinner, M. A. 2017, ApJ, 851, 93 [Google Scholar]
 Klassen, M., Pudritz, R. E., Kuiper, R., Peters, T., & Banerjee, R. 2016, ApJ, 823, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R. 2016, ArXiv eprints [arXiv:1511.03457] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kuiper, R., & Hosokawa, T. 2018, A&A, 616, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., Flock, M., & Klahr, H. 2008, Proc. Int. Astron. Union, 4, 103 [CrossRef] [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ, 722, 1556 [CrossRef] [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010b, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010c, Proc. Int. Astron. Union, 6, 215 [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2012, A&A, 537, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2014, in The Labyrinth of Star Formation (Berlin: Springer), 36, 379 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Larson, R. B., & Starrfield, S. 1971, A&A, 13, 190 [NASA ADS] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spectr. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford: Oxford University Press) [Google Scholar]
 Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41 [Google Scholar]
 Owen, J. E., Ercolano, B., & Clarke, C. J. 2014, Ap&SS, 36, 127 [Google Scholar]
 Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [Google Scholar]
 Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2006, [arXiv:astroph/0603592] [Google Scholar]
 Ramsey, J. P., & Dullemond, C. P. 2015, A&A, 574, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188 [NASA ADS] [CrossRef] [Google Scholar]
 Rosdahl, J., Schaye, J., Teyssier, R., & Agertz, O. 2015, MNRAS, 451, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Rosen, A. L., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2016, MNRAS, 463, 2553 [NASA ADS] [CrossRef] [Google Scholar]
 Rosen, A., Krumholz, M., Oishi, J., Lee, A., & Klein, R. 2017, J. Comput. Phys., 330, 924 [Google Scholar]
 Rosen, A. L., Li, P. S., Zhang, Q., & Burkhart, B. 2019, ApJ, 887, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 149 [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2010, ApJ, 714, L58 [NASA ADS] [CrossRef] [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [NASA ADS] [CrossRef] [Google Scholar]
 Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Yorke, H. W.,& Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Frequencydependent opacities and blackbody spectra for T_{disk} = 300 K and T_{⋆} = 5800 K. Opacities are absorption (blue pluses), scattering (red crosses) and extinction (black dots) coefficients for the dustandgas 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 SiO vibrational transitions, the opacity generally increases with the photon frequency. The opacity at stellarlike radiation frequencies is generally greater than at disklike radiation frequencies. 

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

In the text 
Fig. 3 Radial gas temperature profiles in the midplane of the disk following the test of Pascucci et al. (2004) for τ = 0.1. We comparethe gas temperature computed using MCFOST (black dottedline) and RADMC3D (red dashedline), 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 edgeon 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 dashedline) 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 midplane 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 RADMC3D, the hybrid method (M1+FLD) and the FLD method alone in RAMSES. The integrated optical depth in the disk midplane is τ_{810 nm} = 10^{3} and the stellar temperature is T_{⋆} = 4000 K. 

In the text 
Fig. 9 Time evolution of the main sink mass (topleft panel), accretion rate (topright), disk mass (bottomleft), and outflow mass (bottomright). 

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 HYRTi 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 midplane is τ = 100. Top: radial profile in the disk midplane. 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.