Issue 
A&A
Volume 529, May 2011



Article Number  A35  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201015880  
Published online  28 March 2011 
Radiation hydrodynamics with adaptive mesh refinement and application to prestellar core collapse
I. Methods
^{1}
Max Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
email: benoit@mpiahd.mpg.de
^{2}
Laboratoire AIM, CEA/DSM – CNRS – Université Paris Diderot, IRFU/SAp, 91191 GifsurYvette, France
^{3}
Laboratoire de radioastronomie (UMR 8112 CNRS), École Normale Supérieure et Observatoire de Paris, 24 rue Lhomond, 75231 Paris Cedex 05, France
^{4}
École Normale Supérieure de Lyon, Centre de recherche Astrophysique de Lyon (UMR 5574 CNRS), 46 allée d’Italie, 69364 Lyon Cedex 07, France
^{5}
Universität Zürich, Institut für Theoretische Physik, Winterthurerstrasse 190, 8057 Zürich, Switzerland
^{6}
School of Physics, University of Exeter, Exeter, UK EX4 4QL, UK
Received: 6 October 2010
Accepted: 14 February 2011
Context. Radiative transfer has a strong impact on the collapse and the fragmentation of prestellar dense cores.
Aims. We present the radiationhydrodynamics (RHD) solver we designed for the RAMSES code. The method is designed for astrophysical purposes, and in particular for protostellar collapse.
Methods. We present the solver, using the comoving frame to evaluate the radiative quantities. We use the popular fluxlimited diffusion approximation under the grey approximation (one group of photons). The solver is based on the secondorder Godunov scheme of RAMSES for its hyperbolic part and on an implicit scheme for the radiation diffusion and the coupling between radiation and matter.
Results. We report in detail our methodology to integrate the RHD solver into RAMSES. We successfully test the method in several conventional tests. For validation in 3D, we perform calculations of the collapse of an isolated 1 M_{⊙} prestellar dense core without rotation. We successfully compare the results with previous studies that used different models for radiation and hydrodynamics.
Conclusions. We have developed a full radiationhydrodynamics solver in the RAMSES code that handles adaptive mesh refinement grids. The method is a combination of an explicit scheme and an implicit scheme accurate to the secondorder in space. Our method is well suited for starformation purposes. Results of multidimensional densecorecollapse calculations with rotation are presented in a companion paper.
Key words: hydrodynamics / radiative transfer / methods: numerical / stars: formation / ISM: kinematics and dynamics / stars: lowmass
© ESO, 2011
1. Introduction
Within recent years, starformation calculations have undergone a rapid increase in the variety of the physical models that are included. The coupling between radiative transfer and hydrodynamics has been widely studied for many years and considering different regimes and frames (e.g. Mihalas & Mihalas 1984; Lowrie et al. 2001; Mihalas & Auer 2001; Krumholz et al. 2007b). Radiation hydrodynamics (RHD) methods have been developed in gridbased codes (Stone & Norman 1992; Hayes & Norman 2003; Krumholz et al. 2007a; Kuiper et al. 2010; Sekora & Stone 2010; Tomida et al. 2010) and also in smoothed particles hydrodynamics (SPH) codes (Boss et al. 2000; Whitehouse & Bate 2006; Stamatellos et al. 2007). Most of these studies use the popular fluxlimited diffusion approximation (FLD, e.g. Minerbo 1978; Levermore & Pomraning 1981) approximation to model the radiation transport.
In starformation calculations, the easiest method to take into account radiative transfer is to use a barotropic approximation, which reproduces approximately the thermal behaviour of the gas during the collapse. However, more accurate RHD calculations show that a barotropic equation of state (EOS) cannot account for realistic cooling and heating of the gas (e.g. Boss et al. 2000; Attwood et al. 2009; Commerçon et al., in prep., hereafter Paper II). Recently, using radiationmagnetohydrodynamics calculations, Commerçon et al. (2010) have shown that the barotropic approximation cannot properly account for the combined effects of magnetic field and radiative transfer in the first collapse and in the first core formation. On larger scales, radiative transfer has been found to greatly reduce the fragmentation because of the radiative feedback owing to accretion and protostellar evolution (Bate 2009; Offner et al. 2009).
In this study, we present a new RHD solver based on the FLD approximation, which we integrate in the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). The solver is consistently integrated in the secondorder predictorcorrector Godunov scheme of RAMSES, which we modify to account for the radiative pressure. But we add an implicit solver to handle the radiation diffusion and the coupling between matter and radiation, which involves physical processes on timescales much shorter than the hydrodynamical one. The FLD is easier to implement in an AMR code than more sophisticated methods that would require an additional equation on the first moment of the radiative transfer equation (e.g. M1 model, González et al. 2007). The extension to (ideal) MHD flows presents no particular difficulties and has already been used (Commerçon et al. 2010) based on the solver presented in Fromang et al. (2006) and Teyssier et al. (2006).
The paper is organized as follows: in Sect. 2 we recall the RHD equations in the comoving frame we use and we briefly present the FLD approximation. The RHD solver for the RAMSES code is presented in Sect. 3. In Sect. 4 the method is tested for wellknown test cases. As a final test, RHD dense core collapse calculations with a very high resolution are performed. Section 5 summarizes our work and the main results andr present our assessment of the method’s potential.
2. Radiation hydrodynamics in the fluxlimited diffusion approximation
2.1. Radiation hydrodynamics in the comoving frame
We consider the equations governing the evolution of an inviscid, radiating fluid, where radiative quantities are estimated in the comoving frame and are frequencyintegrated (Mihalas & Mihalas 1984) (1)where ρ is the material density, u is the velocity, P the thermal pressure, σ_{R} is the Rosseland mean opacity, F_{r} is the radiative flux, E the fluid total energy E = ρϵ + 1/2ρu^{2} (ϵ is the internal specific energy), σ_{P} is the Planck opacity, B = B(T) is the Planck function, E_{r} is the radiative energy and P_{r} is the radiation pressure. We see that the term σ_{R}F_{r}/c acts as a radiative force on the material. The material energy lost by emission is transferred into radiation, and radiative energy lost by material absorption is added to the material. To close this system, we need two closure relations: one for the gas and one for the radiation. In this work, we only consider an ideal gas closure relation for the material: P = (γ − 1)e = ρk_{B}T/μm_{H} where γ is the specific heats ratio, μ is the mean molecular weight, and e = ρc_{v}T is the gas internal energy. For the radiation, we use the FLD approximation to close the system of moment equations (Minerbo 1978; Levermore & Pomraning 1981). In this work, we consider only the simplified case of a grey material, where all frequencydependent quantities are integrated over frequency. We cannot use a frequencydependent model for our purpose because of CPU limitation.
In comparison with the laboratory frame formulation, Castor (1972) demonstrated that in the comoving frame, an additional advective flux of the radiation enthalpy is not taken into account. In the dynamic diffusion regime, where the optical depth τ ≫ 1 and (v/c)τ ≫ 1, this radiative flux can dominate the diffusion flux, emission or absorption. For an alternative mixed frame formulation, see Krumholz et al. (2007b). In the lowmass star domain, the main focus of this work, we do not expect to encounter dynamic diffusion situations.
2.1.1. The fluxlimited diffusion approximation
As mentioned earlier, we need a closure relation to solve the moment equations coupled to the hydrodynamics (closed by the perfect gas relation), and such a relation is of prime importance. Many possible choices for the closure relation exist. Among these models, the FLD approximation is one of the simplest ones and uses moment models of radiation transport (Minerbo 1978; Levermore & Pomraning 1981).
Under the fluxlimited diffusion approximation, the radiative flux is expressed directly as a function of the radiative energy and is proportional and colinear to the radiative energy gradient (Fick’s law). Under the grey approximation, we have (2)where λ = λ(R) is the flux limiter, and R =  ∇E_{r}  /(σ_{R}E_{r}). In this study, we retain the flux limiter that has been derived by Minerbo (1978), assuming the intensity as a piecewise linear function of solid angle (3)The flux limiter has the property that λ → 1/3 in optically thick regions and λ → 1/R in optically thin regions. We recover the proper value for diffusion in the optically thick regime, F = −c/(3σ_{R})∇E_{r}, and the flux is limited to cE_{r} in the optically thin regime. Under the FLD approximation, the radiative transfer equation is then replaced by a unique diffusion equation on the radiative energy (4)
3. A multidimensional radiation hydrodynamics solver for RAMSES
3.1. The AMR RAMSES code
We use the RAMSES code (Teyssier 2002), which integrates the equations of ideal magnetohydrodynamics (Fromang et al. 2006; Teyssier et al. 2006) using a secondorder Godunov finite volume scheme. The MHD equations are integrated using a MUSCL predictorcorrector scheme, originally presented in van Leer (1979). Fluxes at the cell interface are estimated with an approximated Riemann solver (LaxFriedrich, HLL, HLLD, etc.). For its AMR grid, RAMSES is based on a “treebased” AMR structure, the refinement is made on a cellbycell basis. Various refinements can be used (fluid variable gradients, instability wavelength, etc.).
The AMR code RAMSES has often been used for starformation purposes (Hennebelle & Fromang 2008; Hennebelle & Teyssier 2008; Hennebelle & Ciardi 2009; Commerçon et al. 2010). Commerçon et al. (2008) have thoroughly and successfully compared its results with standard SPH, showing a good agreement between the methods.
3.2. The Eulerian approach and properties of conservation laws
In Eulerian hydrodynamics, the mesh is fixed and gas density, velocity, and internal energy are primary variables. Eulerian methods fall into two groups: finite difference methods (e.g. the ZEUS code, Stone & Norman 1992; Turner & Stone 2001) and finite volume methods (e.g. the RAMSES code, Teyssier 2002). In the first group, flow variables are conceived as being samples at certain points in space and time. Partial derivatives are then computed from these sampled values and follow Euler equations. In the finite volume approach, flow variables correspond to average values over a finite volume – the cell – and obey the conservation laws in the integral form. Their evolution is determined through Godunov methods by calculating the flux of every conserved quantity across each cell interface.
For an inviscid, compressible flow, the Euler equations in their conservative form read (5)This system can be written in the general hyperbolic conservative form (6)where the vector U = (ρ,ρu,E) contains conservative variables, and the flux vector F(U) = (ρu,ρu ⊗ u + PI,u^{(}E + P^{)}) is a linear function of U.
In this paper, we use the secondorder Godunov method, but applied to the modified Euler equation system under the FLD approximation.
3.3. The conservative radiation hydrodynamics scheme
Let us rewrite the grey RHD equations under the FLD approximation within the comoving frame, taking into account the gravity terms (7)Note that we rewrite the opacity σ_{i} as κ_{i}ρ. The dimension of κ_{i} is cm^{2} g^{1}.
The basic idea is to build a solver for a radiative fluid, with an additional pressure owing to the radiation field: the radiative pressure. Following the Euler equations in their conservative form, the new conservative quantities are density ρ, momentum ρu, total energy E_{T} of the fluid (gas + photon) per unit volume, i.e. ρϵ + ρu^{2}/2 + E_{r}. Primitive hydrodynamical variables do not change for the fluid, but we add a fourth equation for the radiative energy.
In order to facilitate these equations in RAMSES and to minimize the number of changes with the pure hydrodynamical version, we decompose each term where the flux limiter λ appears as follows: λ = 1/3 + (λ − 1/3). We thus distinguish a diffusive part (Eddington approximation, P_{r} = 1/3E_{r}I) and a correction part. The computation of predicted states and fluxes in the MUSCL scheme is made under the Eddington approximation, which is then corrected in an additional corrective step. The RHD equations can be rewritten as (8)The new system ∂_{t}U + ∇F(U) = S(U) is composed of the modified hyperbolic left hand side (LHS) and the right hand side (RHS) source, corrective and coupling terms S(U) = S_{exp} + S_{imp}. The hyperbolic system, as well as the source and corrective terms S_{exp}, are integrated in time with an explicit scheme. The modified RHD hyperbolic system reads (9)This system is used in the predictorcorrector MUSCL temporal integration. To predict states, we consider the worst case, where the radiative pressure is the greatest (1/3E_{r}). For the conservative update (corrector step) we consider the LHS of system (8). The associated eigenvalues corresponding to the three waves are (10)Radiative pressure enlarges the span of solutions, since wave speeds are faster. Once again, with the Eddington approximation, we build the system for the case where the radiative pressure would be the greatest. Therefore, the waves propagate at a speed that is within the wave extrema.
The next step consists in correcting errors due to the Eddington approximation by integrating source terms S_{ne}(11)we here consider an isotropic radiative pressure tensor P_{r} = λE_{r}I. Other authors considered extensions to this closure relation using the Levermore (1984) FLD theory (Turner & Stone 2001; Krumholz et al. 2007b).
To ensure the stability of the explicit step, the CourantFriedrichLevy (CFL) stability condition used to estimate the timestep also takes into account the radiative pressure. The updated CFL condition is simply (12)
3.4. The implicit radiative scheme
The most demanding step in our timesplitting scheme is to deal with the diffusion term and the coupling term κ_{P}ρc(a_{R}T^{4} − E_{r}), which corresponds to S_{imp}. This update has to be made with an implicit scheme, because the time scales of these processes are much shorter than those of pure hydrodynamical processes. Two coupled equations are integrated implicitly (13)which give the implicit scheme on a uniform grid^{1}(14)where ρϵ = C_{v}T. The nonlinear term (T^{n + 1})^{4} makes this scheme difficult to invert. Yet it is much easier to solve implicitly a linear system. Assuming that changes of temperature are small within a time step, we can write (15)eventually, with (14a), we obtain as a function of and . Then can be directly injected in the radiative energy Eq. (14b), and is finally expressed as a function of , , and . The implicit scheme for the radiative energy in a cell of volume V_{i} in the xdirection becomes (16)The gas temperature within a cell is simply given by (17)We compute the Planck and Rosseland opacities and the flux limiter with a gas temperature value given before the implicit update (with index n) to preserve the linearity of the solver.
3.5. Implicit scheme integration with the conjugate gradient algorithm
Equation (16) is solved on a full grid made of N cells. It results in a system of N linear equations, which can be written as a linear system of equations (18)where x is a vector containing radiative energy values. The conjugate gradients (CG) method is one of the most popular nonstationary iterative methods for solving large symmetric systems of linear equations Ax = b. The CG method can be used if the matrix A to be inverted is square, symmetric, and positivedefinite. The CG is memoryefficient (no matrix storage) and runs quickly with sparse matrices. For a N × N matrix, the CG converges in less than N iterations. Basically, the CG method is a steepestgradientdescent method in which descent directions are updated at each iteration. Another advantage is that the CG method can be run easily on parallel machines.
To improve convergence of the CG or even to insure convergence if one deals with an illconditioned matrix A, we use a preconditioning matrix M that approximates A. M is also assumed to be symmetric and positive definite. In this work, we use a simple diagonal preconditioning matrix, which retains only the inverse of A diagonal elements. The convergence of the CG algorithm is estimated following two criteria: estimation of the norm L^{2} (criterion r^{(j)}/ r^{(0)} < ϵ) or estimation of the norm L^{∞} (maximum residual value max { r^{(j)} } /max { r^{(0)} } < ϵ). Values of ϵ typically range from 10^{8} to 10^{3}. In appendix we present an alternative method to the conjugate gradient, the SuperTime Stepping method, which can be used efficiently on uniform grids or in some particular cases.
3.6. Comparison to other schemes
Other RHD solvers based on the FLD approximation have been designed in gridbased codes. Among them, the ZEUS and ORION implementations are the most widely used and discussed. Compared to ZEUS (Stone & Norman 1992; Turner & Stone 2001; Hayes et al. 2006), our method is fundamentally different, although they also use the comoving frame to estimate the radiative quantities. ZEUS code is based on a finite difference scheme, using artificial viscosity and regular grids. Its nonconservative formulation can lead to spurious wave propagation when the resolution is not high enough, or if the radiative pressure dominates the characteristic velocity (the classical Burgers equation problem). All radiative terms such as the radiation transport and the radiative pressure work are integrated implicitly in ZEUS. The implicit scheme is based on a NewtonRaphson method, using GMRES or LU algorithms for the matrix inversion.
ORION is a patchedbased AMR code, which is less flexible than the treebased AMR (Krumholz et al. 2007a). Krumholz et al. (2007a) implemented the mixedframe RHD equations using a multigrid, multitimestep method to solve the implicit scheme for the radiation module (Howell & Greenough 2003). The hydrodynamics part of ORION uses a secondorder conservative Godunov scheme, with approximate Riemann solvers and very little artificial viscosity to treat shocks and discontinuities. Using the same idea as in this study, the diffusion and matterradiation coupling terms are integrated implicitly, while the radiative force and radiative pressure work are integrated explicitly. Contrary to our work, Krumholz et al. (2007a) do not take into account the radiative pressure in the flux estimate at the cell interface for the conservative update, which could also lead to an inaccurate wave speed propagation in radiationpressure dominated regions.
3.7. Implicit scheme on an AMR grid
For studies involving large dynamical ranges, such as star formation, it is necessary to extend our implicit scheme to AMR grids. The difficulty is to compute the correct fluxes and gradients at the interfaces between two cells. We need to carefully consider the energy balance on a given volume. Energy balances are performed on volumes overlapping two cells, which depends on whether the mesh is refined or not. If one considers the face of a cell on a level ℓ; three connecting configurations with other cells are possible (see Fig. 1):

Configuration 1: the neighbouring cell is at the same levelℓ: cells 1 and 3;

Configuration 2: the neighbouring cell is at level ℓ − 1: cells 1 and i;

Configuration 3: 2 neighbouring cells exist at level ℓ + 1: cell i with cells 1 and 2.
Last but not least, the diffusion routine is called only once per coarse step (no multiple timestepping) and scans the full grid, from the finer level to the coarse level. In order to optimise matrixvectors products, we choose to avoid dealing with configuration 3. Hence, when cells at level ℓ + 1 are monitored, values for cells at level ℓ are updated. Configurations 2 and 3 are then performed at the same time. Depending on the configuration, gradients and flux estimates are different. In the following, we will focus on cell 1, of size Δx × Δx.
Fig. 1 Example of AMR grid configuration. 
3.7.1. Gradient estimate
Gradients ∇E_{r} are estimated between the two neighbouring cells centre
3.7.2. Flux estimate
Let S^{ℓ} be the surface of the interface of a cell at level ℓ and the flux across this surface between two cells i and j. The energy rate F × S that is exchanged at this interface is Because access to neighbouring finer cells is not straightforward, we see from configuration 3 all the interest of updating quantities at level ℓ − 1 when scanning grid at level ℓ.
3.8. Limits of the methods
The first drawback is the use of the FLD approximation that implies isotropy of the radiation field. Anisotropies in the transparent regime are not well processed with the FLD, contrary to more accurate models like M1 (González et al. 2007) or VETF (Hayes & Norman 2003). A second limitation comes from the grey opacity assumption, which could limit the accretion on the protostars (see Zinnecker & Yorke 2007,and references therein). In highmass star formation, Yorke & Sonnhalter (2002) show that using a frequencydependent radiative transfer model enhances the flashlight effect and helps to accrete more mass onto the central protostar.
From a technical point of view, our method works only for unique time stepping, i.e. all levels evolve with the same time step. We do not take advantage of the multiple time stepping possibility. As a compromise, we investigate the possibility to evolve finer levels with their own time steps and perform a diffusioncoupling step every 2, 4 or more finer time steps. As a result, we find that performing the diffusion step only every 2 or 4 fine time steps gives correct results. The frequency of the implicit solver calls is left to the user convenience, by use of the mutlitime stepping of RAMSES. For instance, for a grid of levels ranging form level ℓ_{min} to ℓ_{max}, a unique timestep can be used for levels ranging form level ℓ_{min} to ℓ_{i}. In that case, only the levels finer that ℓ_{i} will use not updated radiative quantities in the Godunov solver. A future development would be to use a multigrid solver or preconditioner for parabolic equations (Howell & Greenough 2003).
Another difficulty comes from the residual norm and scalar estimates in the CG algorithm. For a large grid with a large number of cells, the dot product can be dominated by roundoff errors, owing to estimates close to machine precision. This becomes even worse in parallel calculations. The usual MPI function MPI_SUM fails with a large number of processors, the results of any sum becoming a function of number of processors. This dramatically affects the number of iterations. We implemented a new MPI function that performs summation in doubledouble precision following He & Ding (2000), using the Knuth (1997) method.
Eventually, our method involves only immediate neighbouring cells, whatever their refinement level. As a consequence, our method is only first order at the border between levels. This could give rise to a loss of accuracy in diffusion problems, because gradient estimates are not secondorder accurate when neighbouring cells are at finer levels (see configuration 3, Popinet 2003). However, the tests we performed ascertain that the method is still roughly secondorder accurate. The errors are only confined to surfaces much smaller than the total volume.
4. Radiation hydrodynamics solver tests
4.1. 1D test: linear diffusion
We only consider the radiative energy diffusion equation in this test, without either hydrodynamics or coupling with the gas. The equation to integrate is simply (19)Consider a box of length L = 1. The initial radiative energy corresponds to a delta function, namely it is equal to 1 everywhere in the box, except at the centre, where it equals E_{r,L/2}Δx = E_{0} = 1 × 10^{5}. To simplify, we choose ρκ_{R} = 1 and a constant time step. We apply Von Neumann boundary conditions, i.e. zerogradient. The analytical solution in a pdimensional problem is given by (20)where χ = c/(3ρκ_{R}).
Fig. 2 a): Comparison between numerical solution (squares) and analytical solution (red line) at time t = 1 × 10^{12} for the calculations with 16 cells. b): L_{1} norm of the error as a function of h = 1/Δx. The dotted line shows the evolution of the error as a function of h^{2} and the dashed line the evolution of the error as a function of h. 
Figure 2a shows results at time t = 1 × 10^{12} for a resolution of N = 16 cells. The numerical solution is very close to the analytical one, even with this small number of cells. In Fig. 2b we show the evolution of the L_{1} norm of the relative error as a function of h = 1/Δx. The L_{1} norm is estimated as (21)where E_{r,i} is the numerical value of the radiative energy at position x_{i} and time t, and E_{r,a}(x_{i},t) the corresponding analytic value. The error clearly grows as h^{2} (dotted line), which indicates that our method is secondorder accurate.
In Table 1 we report the CPU time, the total number of iterations, and the number of time steps for various numerical resolutions. At low resolution, the number of time steps increases linearly with the number of cells, as well as the CPU time. The number of iterations per time step is constant, i.e. the convergence of CG does not depend on the dimension of the problem, but on the nature of the problem.
CPU time, total number of time steps, and number of iterations per time step for various numbers of cells N.
4.2. 1D test: nonlinear diffusion
In this second test, we consider an initial discontinuity in a box with different initial radiative energy states: E_{r} = 4 on the left and E_{r} = 0.5 on the right. We apply Von Neumann boundary conditions. We integrate the same equation as in the previous test, but with a Rosseland opacity as a nonlinear function of the radiative energy, i.e. . Last, we allow refinement with a criterion based on the radiative energy gradient. In each region where ∇E_{r}/E_{r} > 3%, the grid is refined.
Figure 3a shows the radiative energy profiles at time t = 1.4 × 10^{2} for calculations run with a coarse grid of 16 cells and a maximum effective resolution of 512 cells (squares), and for calculations run with 2048 cells, taken to be the “exact” solution (red curve). Because of the nonlinear opacity, the diffusion is more efficient in the highenergy region. The mean opacity at cell interface is computed using an arithmetic average, which is more adapted for the case of nonlinear opacity. The levels are finer (higher resolution) in highradiative energy gradient regions. Note that we checked that we obtained similar results in a 2D plane parallel case and in a 2D case with an initial step function that maked an angle π/4 with the computational box axis. This validates our routine in the x and y directions.
In Fig. 3b we show the evolution of the L_{1} norm of the error as a function of the mesh spacing. The uniform grid points (diamonds) correspond to a calculations run with a number of cells ranging from 16 to 512 (i.e. ℓ = 4 to ℓ = 9) . When AMR is used (squares), the error is plotted as function of the minimum grid spacing, corresponding to effective resolutions ranging from 32 to 512 cells. The coarse level remains unchanged, ℓ_{min} = 4 (i.e. 16 coarse cells). The advantage of the AMR is clear, the error remains identical compared to the uniform grid case, but the number of cells is greatly reduced (25 cells with ℓ_{max} = 5, 38 cells with ℓ_{max} = 6, 59 cells with ℓ_{max} = 7, 83 cells with ℓ_{max} = 8 and 110 cells with ℓ_{max} = 9). The AMR implementation works well (secondorder accuracy), and does not suffer from the fact that our scheme is only first order in space at the level interface, which validates our scheme used to estimate the gradients at the cell interface in Sect. 3.7. Finally, as in the previous test, the error increases with h^{2}, even if the diffusion problem is nonlinear for radiation.
Fig. 3 Nonlinear diffusion of an initial step function with AMR, the refinement criterion based on radiative energy gradients. a) Radiative energy profiles at time t = 1.4 × 10^{2} (square – numerical solution, “exact” solution in red, run with 2048 cells). The AMR levels (right axis) are plotted in blue. b) L_{1} norm of the error as a function of h = 1/Δx, without AMR (diamond) and with AMR (squares), up to an effective resolution of 512 cells (the error is plotted as a function of the minimum mesh spacing, corresponding to the maximum resolution). The dotted line shows the evolution of the error as a function of h^{2}. 
4.3. Matterradiation coupling test
Another conventional test is the matterradiation coupling. Consider a static, uniform, absorbing fluid initially out of thermal balance, in which the radiation energy E_{r} dominates and is constant. An analytic solution can be obtained for the time evolution of the gas energy e, by solving the ordinary differential equation (Turner & Stone 2001) (22)We performed two tests, with two initial gas energies, e = 10^{10} erg cm^{3} and e = 10^{2} erg cm^{3}. In both tests, the following quantities are taken constant: the radiative energy E_{r} = 1 × 10^{12} erg cm^{3}, the opacity σ = 4 × 10^{8} cm^{1}, the density ρ = 10^{7} g cm^{3}, the mean molecular weight μ = 0.6, and the adiabatic index γ = 5/3. Figure 4 shows the evolution in time of the gas energy for the analytic solution (red line) and the numerical solution (squares). In the first calculations, where the initial gas temperature is greater than the radiative temperature, we used a variable time step Δt that increases with time, starting from 10^{20} s. This good sampling gives very good results. In the second case, we used a constant time step Δt = 10^{12} s. Although the sampling is bad at early times and longer than the cooling time, numerical solutions always match the analytic one. This validates our linearization of the emission term (a_{R}T^{4}) in Eq. (15).
Fig. 4 Matterradiation coupling test. The radiative energy is kept constant, E_{r} = 1 × 10^{12} erg cm^{3}, whereas the initial gas energies are out of thermal balance (e = 10^{2} erg cm^{3} and e = 10^{10} erg cm^{3}). Numerical (square) and analytic (red curve) evolutions of the gas energy are given as a function of time. 
Fig. 5 Left: temperature profiles for a subcritical shock with piston velocity v = 6 km s^{1}, at time t = 3.8 × 10^{4} s. Right: temperature profiles for a supercritical shock with piston velocity v = 20 km s^{1}, at time t = 7.5 × 10^{3} s. In both cases, the temperatures are displayed as a function of z = x − vt. The squares represent the gas temperature and the diamonds the radiative temperature. The red curves represent the gas and radiative temperatures obtained with a calculation using 2048 cells, which we take as the “exact” solution. The AMR levels (blue line – right axis) are overplotted. 
4.4. 1D full RHD tests: radiative shocks
Testing the numerical method for radiative shock calculations is a last important step that every code attempting to integrate RHD equations should perform (Hayes & Norman 2003; Whitehouse & Bate 2006; González et al. 2007). Following the Ensman (1994) initial conditions, we tested our routine for sub and supercritical radiative shocks.
Initial conditions are as follows: uniform density ρ_{0} = 7.78 × 10^{10} g cm^{3} and temperature T_{0} = 10 K. The box length is L = 7 × 10^{10} cm, the opacity is constant (σ = 3.1 × 10^{10} cm^{1}), μ = 1 and γ = 7/5. The 1D homogeneous medium moves with a uniform speed (piston speed) from right to left and the left boundary is a wall. The shock is generated at this boundary and travels backwards. The piston velocity varies, producing sub or supercritical radiative shocks. The AMR is used, and the refinement criterion is based on the density and radiative energy gradients (30%), the grid has 32 coarse cells and we use five levels of refinement. We use the Minerbo flux limiter. The time step is given by the hydrodynamics CFL for the explicit and implicit schemes.
Figure 5 shows the gas and radiative temperatures for sub and supercritical radiative shocks, as a function of z = x − vt, where v is the piston’s velocity. The AMR is used in both calculations. The squares represent the gas temperature and the diamonds the radiative temperature. The red curves represent the gas and radiative temperatures obtained with a calculation using 2048 cells, which we take as the “exact” solution. The subcritical shock is obtained with a piston velocity v = 6 km s^{1}, whereas the supercritical shock is obtained with v = 20 km s^{1}. In both tests, the occurrence of an extended, nonequilibrium radiative precursor is obvious. As expected, pre and postshock gas temperatures are equal in the supercritical case.
For the subcritical case, the postshock gas temperature is given by (Ensman 1994; Mihalas & Mihalas 1984) (23)where ℛ = k/μm_{H} is the perfect gas constant. For our initial setup, this analytic estimate gives T_{2} ~ 810 K. Numerical calculations give T_{2} ~ 825 K at time t = 3.8 × 10^{4} s, which agrees with the analytic estimate comparable to values obtained with more accurate methods (González et al. 2007). The characteristic temperature T_{ − } ~ 275 K immediately in front of the shock agrees very well with the analytic estimate (Mihalas & Mihalas 1984) (24)This means that in front of the shock, the gas internal energy flux flowing downstream is equal to the radiative flux flowing upstream. The entire radiative energy is absorbed upstream and contributes to heat the upstream gas. Similarly, the spike temperature T_{ + } ~ 1038 K also agrees well with the analytic estimate of Mihalas & Mihalas (1984)(25)We note that the AMR scheme enables us to accurately describe the gas temperature spike at the shock. The medium around the spike is optically thin, and the numerical resolution in this region is therefore of crucial importance. The spike’s amplitude varies according to the model used for radiation and to the effective numerical resolution. Thanks to the AMR scheme, the spike’s amplitude is larger in the supercritical case, but not as large as those obtained with M1 or VTEF models (Hayes & Norman 2003; González et al. 2007). However, this last test shows the ability of our timesplitting method to integrate the RHD equations.
4.5. 3D densecorecollapse calculations without rotation
In this section, we perform calculations of a 1 M_{⊙} densecore collapse without rotation, using our grey FLD solver. We compare our FLD results for a model without initial rotation with the ones obtained by Masunaga et al. (1998) and with our results obtained with a 1D code (see Commerçon et al. 2011). We also qualitatively compare our results with the pioneered ones of Larson (1969) and Winkler & Newman (1980). This latter test provides a validation in 3D for starformation purposes.
Summary of first core properties at time t = 1.012 t_{ff} and ρ_{c} = 2.7 × 10^{11} g cm^{3}.
4.5.1. Initial conditions
To facilitate the comparison with other studies, we used the same initial conditions as in Commerçon et al. (2008) and in Paper II and the Lax Friedrich Riemann solver. We chose highly gravitationally unstable initial conditions. The initial sphere is isothermal, T_{0} = 10 K, and has a uniform density ρ_{0} = 1.38 × 10^{18} g cm^{3}. The ratio α between initial thermal energy and gravitational energy is α ~ 0.50. The initial radius is R_{0} = 7.07 × 10^{16} cm. The theoretical freefall time is t_{ff} = 57 kyr. The initial isothermal sound speed is c_{s0} ~ 0.19 km s^{1} and γ = 5/3. The outer region of the sphere is at the same temperature as the core temperature, but is 100 times less dense. The sphere radius is equal to a quarter of the box length to minimize border effects.
We use the set of opacities given by Semenov et al. (2003) for low temperature (<1000 K), which we compute as a function of the gas temperature and density. For each cell we perform a bilinear interpolation on the mixed opacities table. Below 1500 K the opacities are dominated by grain (silicate, iron, troilite, etc.). Semenov et al. (2003) take into account the dependence of the evaporation temperatures of ice, silicates, and iron on the gas density. We here use spherical composite aggregate particles for the grain structure and topology and a normal iron content in the silicates, Fe/(Fe + Mg) = 0.3.
4.5.2. Results
To resolve the Jeans length, we use N_{J} = 10 (i.e. 10 points per Jeans length). Masunaga et al. (1998) showed that the first core properties are independent of the initial conditions for lowmass star formation. We can then compare our results with those obtained by Masunaga et al. (1998), even if we use different initial conditions. We also compare our results with those obtained using a 1D spherical code (Audit et al. 2002) in Commerçon et al. (2011).
Fig. 6 Profiles of density, radial velocity, temperature, optical depth, and integrated mass as a function of the radius and the temperature as a function of density in the 3D computational domain. All values are computed at time t = 1.012t_{ff}. 
Table 2 summarizes the first core properties obtained at time t = 1.012 t_{ff} with RAMSES. The first core radius and mass are qualitatively similar to the results obtained in other 1D Lagrangean calculations (see Masunaga et al. 1998; Commerçon et al. 2011), even though we use a completely different hydrodynamical scheme (e.g. no artificial viscosity, Eulerian, etc.). The first core radius is however a factor 2 greater than the one found in Larson (1969) and Winkler & Newman (1980), who used simplified dust opacity models. Since the first core is mainly set by the opacity, this explains the differences. We define the first core radius as the radius at which the infall velocity is maximal. The accretion rate on the first core is typical of lowmass star formation, ~ 10^{5} M_{⊙}/yr. We note that the value α_{acc} ~ 24 is relatively high, with α_{acc} defined as (where c_{s0} the isothermal sound speed). This indicates that our collapse model is closer to the dynamical LarsonPenston collapse solution (Larson 1969; Penston 1969) than to the classical SIS model of Shu (1977), for which α_{acc} ~ 0.975.
In Fig. 6 we show the profiles of density, radial velocity, temperature, optical depth, and integrated mass as a function of the radius and the temperature as a function of the density in the computational domain at time t = 1.012 t_{ff}. All quantities are mean values in the equatorial plane. In the density profiles, all cells of the calculations have been displayed (blue points). The spread in the density distribution is very small. The spherical symmetry is thus well conserved in the 3D calculations with RAMSES. We compare these profiles with those obtained in Commerçon et al. (2011). The density jump between the first core border and the centre is of the same order of magnitude as for the 1D spherical case. The infall velocity at the shock is also comparable (~2 km s^{1}). The accretion shock takes place around τ ~ 5−10, in the optically thick region. We do not see a jump in temperature through the accretion shock, which is a supercritical radiative shock. Eventually, we see from the temperature versus density plot that the thermal behaviour of the gas is not perfectly adiabatic in the central core. The first core is not fully adiabatic and is able to decrease its entropy level by radiating in the upstream material. The slight kink in the curve at T ~ 80 K (log(T) ~ 1.7) corresponds to ice evaporation in the opacity table. The opacity decreases abruptly, this is the reason why the cooling is more efficient in that region.
5. Summary and perspectives
We have developed a full radiationhydrodynamics solver using the fluxlimited diffusion approximation, which is integrated in the AMR RAMSES code. Our solver uses a timesplitting integrator scheme and combines explicit and implicit methods. Each step of the integration is detailed in this work. The method was successfully tested in several conventional tests in 1D and 2D. We demonstrated that our method is secondorder accurate in space, even when AMR is used. We also performed collapse calculations of a nonrotating dense core and successfully compared our results with those obtained by Masunaga et al. (1998) and Commerçon et al. (2011), which are based on different methods in 1D spherical codes. Our method has thus been demonstrated to be robust and well suited for star formation. In Paper II we present detailed RHD calculations with a very high resolution of densecore collapse in rotation. We showed that our method enables us to accurately handle the heating and cooling processes. Last but not least, we extended our method to the radiationmagnetohydrodynamics flows in Commerçon et al. (2010).
The next step following this work will be to tune our solver for adaptive timestepping to make use of all benefits of the AMR in RAMSES. For example, the next stages of the collapse, the second collapse and the second core formation, require a huge amount of numerical resolution and the dynamical timescale becomes much shorter. An adaptive timestep scheme is then suitable. Another improvement is to use a multigroup approach in the radiation solver. Some attempts have been presented in the literature (e.g. Shestakov & Offner 2008), but the computational cost remains too high nowadays compared to the grey model.
Acknowledgments
The calculations were performed at CEA on the DAPHPC cluster. The research of B.C. is supported by the postdoctoral fellowships from MaxPlanckInstitut für Astronomie. The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/20072013 Grant Agreement No. 247060). B.C. also thanks Neal Turner for useful discussions.
References
 Alexiades, V., Amiez, G., & Gremaud, P.A. 1996, Com. Num. Meth. Eng, 12, 12 [Google Scholar]
 Attwood, R. E., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2009, A&A, 495, 201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Audit, E., Charrier, P., Chièze, J., & Dubroca, B. 2002, [arXiv:0206281] [Google Scholar]
 Bate, M. R. 2009, MNRAS, 392, 1363 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P., Fisher, R. T., Klein, R. I., & McKee, C. F. 2000, ApJ, 528, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I. 1972, ApJ, 178, 779 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Audit, E., Chabrier, G., & Chièze, J. 2011, A&A, accepted [arXiv:1102.2921] [Google Scholar]
 Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188 [NASA ADS] [CrossRef] [Google Scholar]
 He, Y., & Ding, C. H. Q. 2000, in ICS ’00: Proceedings of the 14th international conference on Supercomputing (New York, NY, USA: ACM), 225 [Google Scholar]
 Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Ciardi, A. 2009, A&A, 506, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Howell, L. H., & Greenough, J. A. 2003, J. Comput. Phys., 184, 53 [Google Scholar]
 Knuth, D. E. 1997, The art of computer programming, 3rd edn., seminumerical algorithms (Boston, MA, USA: AddisonWesley Longman Publishing Co., Inc.), 2 [Google Scholar]
 Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, ApJ, 656, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007b, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spec. Radiat. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Lowrie, R. B., Mihalas, D., & Morel, J. E. 2001, J. Quant. Spec. Radiat. Transf., 69, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Masunaga, H., Miyama, S. M., & Inutsuka, S.I. 1998, ApJ, 495, 346 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mihalas, D., & Auer, L. H. 2001, J. Quantit. Spec. Radiat. Transf., 71, 61 [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics, ed. D. Mihalas, & B. W. Mihalas [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541 [NASA ADS] [CrossRef] [Google Scholar]
 Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131 [NASA ADS] [CrossRef] [Google Scholar]
 O’Sullivan, S., & Downes, T. P. 2006, MNRAS, 366, 1329 [NASA ADS] [CrossRef] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Popinet, S. 2003, J. Comput. Phys., 190, 572 [NASA ADS] [CrossRef] [Google Scholar]
 Sekora, M. D., & Stone, J. M. 2010, J. Comput. Phys., 229, 6819 [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]
 Shestakov, A. I., & Offner, S. S. 2008, J. Comput. Phys., 227, 2154 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S. 2007, A&A, 475, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2010, ApJ, 714, L58 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1979, J. Comput. Phys., 32, 101 [Google Scholar]
 Whitehouse, S. C., & Bate, M. R. 2006, MNRAS, 367, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Winkler, K., & Newman, M. J. 1980, ApJ, 238, 311 [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]
Appendix A: The supertime stepping versus the conjugate gradient
In this appendix we present the supertime stepping (STS) method. It is used to solve parabolic equation systems, like the conjugate gradient (CG) we used previously. We implement the STS scheme into RAMSES. We compare the CG and the STS methods for the particular case of the 1D linear diffusion test presented in Sect. 4.1.
A.1. The supertime stepping
The STS is a very simple and effective way to speed up explicit timestepping schemes for parabolic problems. The method has been recently rediscovered in Alexiades et al. (1996), but it remains relatively unknown in computational astrophysics (Mignone et al. 2007; O’Sullivan & Downes 2006). The STS frees the explicit scheme from stability restrictions on the timestep. It can be very powerful in some cases and is easy to implement compared to implicit methods that involve matrix inversions.
The STS is designed for time dependent problem such as (A.1)where A is a square, symmetric positive definite matrix. Equation (A.1) is rewritten with the corresponding standard explicit scheme (A.2)the explicit scheme is subject to the restrictive stability condition (A.3)where ρ(·) denotes the spectral radius. The equivalent CFL condition is (A.4)where λ_{max} stands for the highest eigenvalue of A. For the 1D heat equation ∂u/∂t = χΔu, discretized by standard secondorder differences on a uniform mesh, we have λ_{max} = 4χΔx^{2} (Δt_{expl} = Δx^{2}/2χ).
In the STS method, the restrictive stability condition is relaxed by requiring the stability at the end of a cycle of N_{sts} timesteps instead of requiring stability at the end of each time step Δt. It leads to a RungeKuttalike method with N_{sts} stages. Following Alexiades et al. (1996), we introduce a superstep consisting of N_{sts}substeps τ_{1},τ_{2},...,τ_{Nsts}. The idea is to ensure stability over the superstep ΔT, while trying to maximize its duration. The inner values, estimated after each τ_{j}, should only be considered as intermediate calculations. Only the values at the end of the superstep approximate the solution of the problem.
The new algorithm can be written as (A.5)and the corresponding stability condition is (A.6)In order to find ΔT as high as possible, the properties of Chebyshev polynomials are exploited, providing a set of optimal values for the substeps given by (A.7)where ν_{sts} is a damping factor that should satisfy 0 < ν_{sts} < λ_{min}/λ_{max}. The superstep ΔT is given by (A.8)Note that as ν_{sts} → 0. The method is unstable in the limit ν_{sts} = 0. The STS method is thus almost N_{sts} times faster than the standard explicit scheme. When ΔT is taken to be the advective (CFL) time step Δt while coupling with the hydrodynamics, the STS requires only approximately (Δt/Δt_{expl})^{1/2} iterations rather than Δt/Δt_{expl} with an explicit scheme.
Fig. A.1 Comparison of the numerical solutions using STS or CG with the analytic one (black line) at time t = 1 × 10^{13} s. The color curves depict numerical solutions obtained with timestep Δt equals to 1 × 10^{13} (blue), 5 × 10^{14} (red), 1 × 10^{14} (green), 1 × 10^{15} (yellow) and 1 × 10^{16} (cyan). 
A.2. The STS implementation for the FLD equation
The STS scheme replaces the implicit radiative scheme presented in Sect. 3.4. Equations of system (14) written with an explicit scheme become (A.9)the explicit time step Δt_{expl} is estimated using values at time n. The next step consists of estimating values of N_{sts} and ν_{sts}, the latter depending on the spectral properties of A. However, as mentioned in Alexiades et al. (1996), it is not required to have a precise knowledge of the spectral properties for the method to be robust. N_{sts} and ν_{sts} are thus arbitrary chosen by the user. Instead of executing one time step of length Δt_{expl}, one executes supersteps of length ΔT. N_{sts} substeps τ_{1},τ_{2},...,τ_{Nsts} are thus performed without outputing until the end of each superstep. When the STS is coupled to the hydrodynamics solver, the cycle is repeated until the time step, given by the hydrodynamical CFL condition, is reached. Superstep ΔT and substeps τ_{i} are reestimated at the end of each cycle.
A.3. Comparison with the conjugate gradient method
To compare the STS with the CG algorithm we used throughout, we consider the test case presented in Sect. 4.1. The equation to integrate is simply (A.10)The initial setup is identical to those in Sect. 4.1. It consists of an initial pulse of radiative energy in the middle of the box. We present here calculations made with either the STS method or the CG algorithm. In both cases, CG and STS are applied over an arbitrary time step Δt that simulates the time step that would be given by the hydro CFL. All calculations were performed on a grid made of 1024 cells. In the STS calculations, for each value of Δt, calculations have been performed using various values of N_{sts} and ν_{sts}. For the CG method, only the convergence criterion ϵ_{conv} changes.
Figure A.1 shows the radiative energy profiles at time t = 1 × 10^{13} s for all calculations we performed. In all panels, the analytic solution is plotted (black line). The two upper plots give results for the CG method. For Δt ≥ 10^{14}, the accuracy is very limited. We also see that for Δt ≥ 10^{13}, the diffusion wave does not propagate at the right speed. The total energy is conserved, but the diffusion wave does not have the correct extent. But the STS results are much more accurate, except for N_{sts} = 20 and ν_{sts} = 1 × 10^{6}. By construction, STS is expected to be more accurate. The stability is poor when N_{sts} = 20 and ν_{sts} = 1 × 10^{6} because ν_{sts} is close to the stability limit (see Alexiades et al. 1996).
Fig. A.2 Comparison of calculations done using STS or CG and a variable time step given by Δt = 1 × 10^{16} ∗ 1.05^{istep}, where istep is the index of the number of global (hydro) time steps. Results are given at time t = 1 × 10^{13} s. 
In Table A.1 we give the CPU time and N_{iter}, which corresponds to the number of iterations for the CG and to the number of substeps for the STS. The number of operations per iteration in the CG and per substep in the STS are equivalent, since it involves the same number of cells (1024). The CPU time spent with the STS is ten times smaller than the one of the CG method. The STS also requires often twice less iterations than the CG. The bottom lines give the results for calculations made with a variable time step, which increases with time. The corresponding profiles are plotted in Fig. A.2. The STS remains more accurate in this case than the CG, which is quite accurate over more than three orders of magnitude. The CG gives good results, because, thanks to the variable time steps, the diffusion wave propagates at a correct speed. Indeed, at t = 0, the gradient of radiative energy is steep and the diffusion wave speed is very high. Using an initial short time step Δt = 10^{16} enables us to be closer to the CFL condition associated to the diffusion wave speed. Then, radiative energy gradients and the former CFL condition relax and the time step can increase with time. This relaxation on the integration time step enables us to maximise the accuracy of implicit methods using a subcycling scheme based on the diffusion wave speed propagation. However, this speed remains quite difficult to estimate.
Fig. A.3 Contours in the equatorial plane of the ratio between diffusion and freefall times for collapse calculations. The diffusion time is estimated as , where l is the local Jeans length. 
Eventually, we must conclude by pointing out that even if the STS method is well adapted for this problem, it remains very limited for starformation calculations. Indeed, the diffusion time is very short compared to the dynamical time estimated as the freefall time (see Fig. A.3) and then, the STS requires too many substeps. The convergence of the CG depends on the nature of the problem and is not affected by strong differences between the diffusion and the dynamical times. Moreover, we never encounter these steep gradients in the radiative energy distribution in starformation calculations. The STS could be efficient only within the fragments, where the diffusion time is very long. This is the reason why we only use the CG method throughout. An alternative but nontrivial solution would be to couple the CG and the STS methods.
All Tables
CPU time, total number of time steps, and number of iterations per time step for various numbers of cells N.
Summary of first core properties at time t = 1.012 t_{ff} and ρ_{c} = 2.7 × 10^{11} g cm^{3}.
All Figures
Fig. 1 Example of AMR grid configuration. 

In the text 
Fig. 2 a): Comparison between numerical solution (squares) and analytical solution (red line) at time t = 1 × 10^{12} for the calculations with 16 cells. b): L_{1} norm of the error as a function of h = 1/Δx. The dotted line shows the evolution of the error as a function of h^{2} and the dashed line the evolution of the error as a function of h. 

In the text 
Fig. 3 Nonlinear diffusion of an initial step function with AMR, the refinement criterion based on radiative energy gradients. a) Radiative energy profiles at time t = 1.4 × 10^{2} (square – numerical solution, “exact” solution in red, run with 2048 cells). The AMR levels (right axis) are plotted in blue. b) L_{1} norm of the error as a function of h = 1/Δx, without AMR (diamond) and with AMR (squares), up to an effective resolution of 512 cells (the error is plotted as a function of the minimum mesh spacing, corresponding to the maximum resolution). The dotted line shows the evolution of the error as a function of h^{2}. 

In the text 
Fig. 4 Matterradiation coupling test. The radiative energy is kept constant, E_{r} = 1 × 10^{12} erg cm^{3}, whereas the initial gas energies are out of thermal balance (e = 10^{2} erg cm^{3} and e = 10^{10} erg cm^{3}). Numerical (square) and analytic (red curve) evolutions of the gas energy are given as a function of time. 

In the text 
Fig. 5 Left: temperature profiles for a subcritical shock with piston velocity v = 6 km s^{1}, at time t = 3.8 × 10^{4} s. Right: temperature profiles for a supercritical shock with piston velocity v = 20 km s^{1}, at time t = 7.5 × 10^{3} s. In both cases, the temperatures are displayed as a function of z = x − vt. The squares represent the gas temperature and the diamonds the radiative temperature. The red curves represent the gas and radiative temperatures obtained with a calculation using 2048 cells, which we take as the “exact” solution. The AMR levels (blue line – right axis) are overplotted. 

In the text 
Fig. 6 Profiles of density, radial velocity, temperature, optical depth, and integrated mass as a function of the radius and the temperature as a function of density in the 3D computational domain. All values are computed at time t = 1.012t_{ff}. 

In the text 
Fig. A.1 Comparison of the numerical solutions using STS or CG with the analytic one (black line) at time t = 1 × 10^{13} s. The color curves depict numerical solutions obtained with timestep Δt equals to 1 × 10^{13} (blue), 5 × 10^{14} (red), 1 × 10^{14} (green), 1 × 10^{15} (yellow) and 1 × 10^{16} (cyan). 

In the text 
Fig. A.2 Comparison of calculations done using STS or CG and a variable time step given by Δt = 1 × 10^{16} ∗ 1.05^{istep}, where istep is the index of the number of global (hydro) time steps. Results are given at time t = 1 × 10^{13} s. 

In the text 
Fig. A.3 Contours in the equatorial plane of the ratio between diffusion and freefall times for collapse calculations. The diffusion time is estimated as , where l is the local Jeans length. 

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.