© ESO, 2011

1. Introduction

Within recent years, star-formation 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 grid-based 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 flux-limited diffusion approximation (FLD, e.g. Minerbo 1978; Levermore & Pomraning 1981) approximation to model the radiation transport.

In star-formation 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 radiation-magnetohydrodynamics 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 second-order predictor-corrector 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 well-known 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 flux-limited 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 frequency-integrated (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² (ϵ 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 σ_RF_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_BT/μm_H where γ is the specific heats ratio, μ is the mean molecular weight, and e = ρc_vT 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 frequency-dependent quantities are integrated over frequency. We cannot use a frequency-dependent 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 low-mass star domain, the main focus of this work, we do not expect to encounter dynamic diffusion situations.

2.1.1. The flux-limited 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 flux-limited 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 | /(σ_RE_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 second-order Godunov finite volume scheme. The MHD equations are integrated using a MUSCL predictor-corrector scheme, originally presented in van Leer (1979). Fluxes at the cell interface are estimated with an approximated Riemann solver (Lax-Friedrich, HLL, HLLD, etc.). For its AMR grid, RAMSES is based on a “tree-based” AMR structure, the refinement is made on a cell-by-cell basis. Various refinements can be used (fluid variable gradients, instability wavelength, etc.).

The AMR code RAMSES has often been used for star-formation 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 second-order 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² 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 + 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_rI) 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 ∂_tU + ∇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 predictor-corrector 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_rI. 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 Courant-Friedrich-Levy (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 time-splitting scheme is to deal with the diffusion term and the coupling term κ_Pρc(a_RT⁴ − 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¹(14)where ρϵ = C_vT. The nonlinear term (T^n + 1)⁴ 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 x-direction 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 non-stationary 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 positive-definite. The CG is memory-efficient (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 steepest-gradient-descent 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 ill-conditioned 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² (criterion  ||r^(j)||/|| r⁽⁰⁾||  < ϵ) or estimation of the norm L^∞ (maximum residual value max { r^(j) } /max { r⁽⁰⁾ }  < ϵ). Values of ϵ typically range from 10^-8 to 10^-3. In appendix  we present an alternative method to the conjugate gradient, the Super-Time 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 grid-based 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 non-conservative 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 Newton-Raphson method, using GMRES or LU algorithms for the matrix inversion.

ORION is a patched-based AMR code, which is less flexible than the tree-based AMR (Krumholz et al. 2007a). Krumholz et al. (2007a) implemented the mixed-frame RHD equations using a multi-grid, multi-timestep method to solve the implicit scheme for the radiation module (Howell & Greenough 2003). The hydrodynamics part of ORION uses a second-order 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 matter-radiation 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 radiation-pressure 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 matrix-vectors 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 high-mass star formation, Yorke & Sonnhalter (2002) show that using a frequency-dependent 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 diffusion-coupling 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 mutli-time 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 round-off 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 double-double 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 second-order 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 second-order 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₀ = 1 × 10⁵. To simplify, we choose ρκ_R = 1 and a constant time step. We apply Von Neumann boundary conditions, i.e. zero-gradient. The analytical solution in a p-dimensional 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): L1 norm of the error as a function of h = 1/Δx. The dotted line shows the evolution of the error as a function of h2 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₁ norm of the relative error as a function of h = 1/Δx. The L₁ 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² (dotted line), which indicates that our method is second-order 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.

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 high-energy 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 high-radiative 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₁ 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 (second-order 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², 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) L1 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 h2.

4.3. Matter-radiation coupling test

Another conventional test is the matter-radiation 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¹⁰ erg cm^-3 and e = 10² erg cm^-3. In both tests, the following quantities are taken constant: the radiative energy E_r = 1 × 10¹² 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_RT⁴) in Eq. (15).

Fig. 4 Matter-radiation coupling test. The radiative energy is kept constant, Er = 1 × 1012 erg cm-3, whereas the initial gas energies are out of thermal balance (e = 102 erg cm-3 and e = 1010 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 × 104 s. Right: temperature profiles for a supercritical shock with piston velocity v = 20 km s-1, at time t = 7.5 × 103 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 super-critical radiative shocks.

Initial conditions are as follows: uniform density ρ₀ = 7.78 × 10^-10 g cm^-3 and temperature T₀ = 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 super-critical 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 super-critical 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, non-equilibrium radiative precursor is obvious. As expected, pre- and post-shock 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₂ ~ 810 K. Numerical calculations give T₂ ~ 825 K at time t = 3.8 × 10⁴ 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 time-splitting method to integrate the RHD equations.

4.5. 3D dense-core-collapse calculations without rotation

In this section, we perform calculations of a 1 M_⊙ dense-core 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 star-formation purposes.

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₀ = 10 K, and has a uniform density ρ₀ = 1.38 × 10^-18 g cm^-3. The ratio α between initial thermal energy and gravitational energy is α ~ 0.50. The initial radius is R₀ = 7.07 × 10¹⁶ cm. The theoretical free-fall 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 low-mass 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.012tff.

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 low-mass 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 Larson-Penston 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 radiation-hydrodynamics solver using the flux-limited diffusion approximation, which is integrated in the AMR RAMSES code. Our solver uses a time-splitting 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 second-order accurate in space, even when AMR is used. We also performed collapse calculations of a non-rotating 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 dense-core 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 radiation-magnetohydrodynamics flows in Commerçon et al. (2010).

The next step following this work will be to tune our solver for adaptive time-stepping 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 time-step scheme is then suitable. Another improvement is to use a multi-group 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 Max-Planck-Institut 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/2007-2013 Grant Agreement No. 247060). B.C. also thanks Neal Turner for useful discussions.

References

Appendix A: The super-time stepping versus the conjugate gradient

In this appendix we present the super-time 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 super-time stepping

The STS is a very simple and effective way to speed up explicit time-stepping 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 time-step. 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 second-order differences on a uniform mesh, we have λ_max = 4χΔx² (Δt_expl = Δx²/2χ).

In the STS method, the restrictive stability condition is relaxed by requiring the stability at the end of a cycle of N_sts time-steps instead of requiring stability at the end of each time step Δt. It leads to a Runge-Kutta-like method with N_sts stages. Following Alexiades et al. (1996), we introduce a superstep consisting of N_stssubsteps τ₁,τ₂,...,τ_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 τ₁,τ₂,...,τ_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 re-estimated 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.05istep, 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 free-fall 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 star-formation calculations. Indeed, the diffusion time is very short compared to the dynamical time estimated as the free-fall 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 star-formation 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 non-trivial solution would be to couple the CG and the STS methods.

All Tables

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): L1 norm of the error as a function of h = 1/Δx. The dotted line shows the evolution of the error as a function of h2 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) L1 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 h2.

In the text

	Fig. 4 Matter-radiation coupling test. The radiative energy is kept constant, Er = 1 × 1012 erg cm-3, whereas the initial gas energies are out of thermal balance (e = 102 erg cm-3 and e = 1010 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 × 104 s. Right: temperature profiles for a supercritical shock with piston velocity v = 20 km s-1, at time t = 7.5 × 103 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.012tff.
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.05istep, 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 free-fall times for collapse calculations. The diffusion time is estimated as , where l is the local Jeans length.
In the text