A&A 475, 37-49 (2007)
DOI: 10.1051/0004-6361:20077373
D. Stamatellos^{1} - A. P. Whitworth^{1} - T. Bisbas^{1} - S. Goodwin^{2}
1 - School of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, Wales, UK
2 - Department of Physics and Astronomy, The University of Sheffield,
Hicks Building, Hounsfield Road, Sheffield S3 7RH, UK
Received 28 February 2007 / Accepted 18 August 2007
Abstract
Aims. We introduce and test a new and highly efficient method for treating the thermal and radiative effects influencing the energy equation in SPH simulations of star formation.
Methods. The method uses the density, temperature and gravitational potential of each particle to estimate a mean optical depth, which then regulates the particle's heating and cooling. The method captures - at minimal computational cost - the effects of (i) the rotational and vibrational degrees of freedom of H_{2}; (ii) H_{2} dissociation and H
ionisation; (iii) opacity changes due to ice mantle melting, sublimation of dust, molecular lines, H^{-}, bound-free and free-free processes and electron scattering; (iv) external irradiation; and (v) thermal inertia.
Results. We use the new method to simulate the collapse of a
cloud of initially uniform density and temperature. At first, the collapse proceeds almost isothermally (
;
cf. Larson 2005, MNRAS, 359, 211). The cloud starts heating fast when the optical depth to the cloud centre reaches unity (
). The first core forms at
and steadily increases in mass. When the temperature at the centre reaches
,
molecular hydrogen starts to dissociate and the second collapse begins, leading to the formation of the second (protostellar) core. The results mimic closely the detailed calculations of Masunaga & Inutsuka (2000, ApJ, 531, 350). We also simulate (i) the collapse of a
cloud, which initially has uniform density and temperature, (ii) the collapse of a
rotating cloud, with an m=2 density perturbation and uniform initial temperature, and (iii) the smoothing of temperature fluctuations in a static, uniform density sphere. In all these tests the new algorithm reproduces the results of previous authors and/or known analytic solutions. The computational cost is comparable to a standard SPH simulation with a simple barotropic equation of state. The method is easy to implement, can be applied to both particle- and grid-based codes, and handles optical depths
.
Key words: stars: formation - methods: numerical - radiative transfer - hydrodynamics - ISM: clouds
Smoothed Particle Hydrodynamics (SPH) (Lucy 1977; Gingold & Monaghan 1977) is a Lagrangian method which invokes a large ensemble of particles to describe a fluid, by assigning properties such as mass, m_{i}, position, , and velocity, , to each particle, i. Intensive thermodynamic variables like density and pressure (and their derivatives) are estimated using local averages (for reviews see Benz 1990, 1991; Monaghan 1992, 2005).
Radiative transfer (RT) has only recently been included in SPH codes (Oxley & Woolfson 2003; Whitehouse & Bate 2004, 2006; Whitehouse et al. 2005; Viau et al. 2005; Mayer et al. 2007). These SPH-RT codes use different simplifying assumptions in order to by-pass the full treatment of multi-frequency radiative transfer in 3 dimensions (a task which is not possible with current computing resources), but they still tend to be computationally expensive. Indeed, even the treatment of full 3D radiative transfer on a single snapshot during the evolution of a simulation is computationally quite expensive (Stamatellos & Whitworth 2005; Stamatellos et al. 2005).
More often, in SPH simulations of star formation, it is standard practice to use a barotropic equation of state, i.e. to put (e.g. Bonnell 1994; Whitworth et al. 1995; Bate 1998). The form of is chosen to mimic the thermodynamics of star forming gas, as revealed by computations of the spherically symmetric collapse of a single, isolated protostar (e.g. Boss & Myhill 1992; Masunaga & Inutsuka 2000).
This is not an ideal situation. (a) A barotropic equation of state is unable to account for the fact that the thermal history of a protostar depends sensitively on its environment, geometry and mass; for example, low-mass protostars remain optically thin to their cooling radiation to higher densities than high-mass ones. Thus, the evolution of the density and temperature cannot be approximated by a single barotropic equation for every system. Even for the same system, the density and temperature evolution away from the centre of the cloud does not follow the corresponding evolution at the centre of the cloud (Whitehouse & Bate 2006). (b) A barotropic equation of state is unable to capture thermal inertia effects (i.e. situations where the evolution is controlled by the thermal timescale, rather than the dynamical one). Such effects appear to be critical at the stage when fragmentation occurs (e.g. Boss et al. 2000).
One of the ultimate goals of star formation simulations is to track the thermal history of star-forming gas. Strictly speaking, this requires a computational method which can treat properly, in 3 dimensions, the time-dependent radiation transport which controls the energy equation. However, this is computationally very expensive, significantly more expensive than the hydrodynamics. We have therefore developed a new algorithm which enables us to distinguish the thermal behaviours of protostars of different mass, in different environments, with different metallicities, and to capture thermal inertia effects, without treating in detail the associated radiation transport. The method uses the density, temperature and gravitational potential of each SPH particle (which are all calculated using the standard SPH formalism) to estimate a characteristic optical depth for each particle. This optical depth then regulates how each particle heats and cools.
The paper is organized as follows. In Sect. 2, we present the new algorithm we have developed to treat the energy equation within SPH, also describing the aspects of SPH that are required for a complete picture of the new method. In Sect. 3 we outline the properties we adopt for the gas and dust in a star-forming cloud, i.e. the composition, energy equation, equation of state, and opacity. In Sect. 4, we present a simulation of the collapse of a molecular cloud; we describe in detail the different stages of the collapse, and compare our results with those of Masunaga & Inutsuka (2000). In Sect. 5, three additional tests are presented: (i) the collapse of a molecular cloud (Boss & Myhill 1992; Whitehouse & Bate 2006), (ii) the collapse of a rotating molecular cloud, with an m=2 density perturbation (Boss & Bodenheimer 1979; Whitehouse & Bate 2006), and (iii) the smoothing of temperature fluctuations in a static, uniform-density sphere (Spiegel 1957). Finally, in Sect. 6, we summarise the method and the tests performed, and discuss the applicability of the new algorithm to simulations of astrophysical systems.
The key to the new method is to use an SPH particle's density, , temperature, T_{i}, and gravitational potential, , to estimate a mean optical depth, , for the SPH particle. This mean optical depth then regulates the SPH particle's radiative heating and cooling; in other words, it determines the extent to which the SPH particle is shielded from external radiation, and the extent to which the SPH particle's cooling radiation is trapped. (The gravitational potential is used here, purely because gravity is the only particle parameter which is already calculated by the SPH code but is not a local function of state. Therefore it should, in some very general sense, represent the larger-scale environment surrounding the SPH particle.)
Specifically, each SPH particle is treated as if it were embedded in a spherically-symmetric pseudo-cloud (its personal pseudo-cloud). The density and temperature profiles of the pseudo-cloud are modelled with a polytrope of index n=2, but the pseudo-cloud is not assumed to be in hydrostatic balance. (We will show later that the choice of n is not critical.)
The position of the SPH particle within its pseudo-cloud is not specified; instead we take a mass-weighted average over all possible positions (see Figs. 1 and 2). For any given position of the SPH particle within the pseudo-cloud, the central density, , and scale-length, , are chosen to reproduce the density and gravitational potential at the position of the SPH particle. Similarly, the pseudo-cloud's central temperature, is chosen to match the temperature at the position of the SPH particle. (Because the pseudo-cloud is not necessarily in hydrostatic equilibrium, we cannot - in general - write , where is the mean gas-particle mass and is Boltzmann's constant.)
The optical depth, , is then calculated by integrating out along a radial line from the given position to the edge of the pseudo-cloud, i.e. through the cooler and more diffuse outer parts of the pseudo-cloud. In this way, samples the different opacity regimes that are likely to surround the SPH particle. This accounts for the fact that, even if the opacity at the position of the SPH particle is low, its cooling radiation may be trapped by cooler more opaque material in the surroundings.
Finally, is obtained by taking a mass-weighted average over all possible positions within the pseudo-cloud.
We use the DRAGON SPH code (Goodwin et al. 2004a,b). DRAGON uses variable smoothing lengths, h_{i}, adjusted so that the number of neighbours is exactly
;
it is important to have a constant number of neighbours to minimise numerical diffusion (Attwood et al. 2007). An octal tree is used to collate neighbour lists and calculate gravitational accelerations, which are kernel-softened using particle smoothing lengths. Standard artificial viscosity is invoked in converging regions, and multiple particle time-steps are used.
Figure 1: Schematic representation of the pseudo-cloud around an SPH particle.The location of the SPH particle inside its pseudo-cloud is not specified. | |
Open with DEXTER |
Figure 2: Density and temperature profiles for a polytropic pseudo-cloud with n=2. The SPH particle could be located anywhere in the cloud (solid and dashed line circles). | |
Open with DEXTER |
The density at the position of SPH particle i is given by a sum over its neighbours, j,
(1) |
The equation of motion for SPH particle i is
(3) |
(4) |
The second term on the righthand side of Eq. (2) is the gravitational acceleration experienced by particle i, given by
The energy equation for particle i is
Suppose that SPH particle i is embedded at radius
,
in a pseudo-cloud with central density
,
scale-length
,
and polytropic index n (see Figs. 1 and 2).
is thus a dimensionless radius, and
and
are chosen so as to reproduce - at this radius - the actual density and gravitational potential of the SPH particle, i.e.
= | (10) |
= | (11) |
(12) |
If we fix n (and hence the forms of
and ), and we pick an arbitrary value for
(modulo that it must be within the pseudo-cloud, i.e.
), then we obtain
= | (13) | ||
= | (14) |
= | T_{i}, | (15) | |
= | (16) |
= | |||
= | (17) |
(19) |
We calculate the pseudo-mean optical depth in the same way as the pseudo-mean column-density. If the Rosseland-mean opacity is a function of density and temperature only,
,
the radial optical depth from radius
to the boundary of the pseudo-cloud is
= | |||
= | |||
(20) |
(22) |
The net radiative heating for SPH particle i is given by
The positive term in Eq. (24) - the one involving
- represents radiative heating due to the background radiation field with effective temperature
.
This term ensures that the SPH particle does not cool radiatively below
.
In a simulation which includes stars - either pre-existing, or formed as an outcome of the simulation; and with luminosities
and positions
- we set
(25) |
(i) If
,
we are in the optically thin cooling regime and Eq. (24) approximates to
(26) |
(ii) If
,
we are in the optically thick cooling regime and Eq. (24) approximates to
To see that this is just the diffusion approximation, suppose that the pseudo-cloud has pseudo-mass
and pseudo-radius
.
Equation (27) then reduces to
(28) |
In order to avoid very short time-steps, we use the following scheme to update the internal energy, u_{i}. From SPH we know the net compressive plus viscous heating rate,
(33) |
On the other hand, if
,
Eq. (32) approximates to
(34) |
The method is easy to implement. In practice we do the following, for each SPH particle i, at each timestep:
Although the method is very efficient, it evidently has limitations, in particular:
Although metals make essential contributions to the opacity, they make very little contribution to the equation of state. Therefore, for the purpose of treating the gas-phase chemistry, we assume that the gas is 70% hydrogen and 30% helium by mass: X=0.7, Y=0.3, Z=0 . At low temperatures, hydrogen is molecular, but as the temperature increases it becomes dissociated and then ionised. At low temperatures, helium is neutral atomic, but as the temperature increass it becomes ionised, first to He^{+}, and then to He^{++}. The relative abundances of these constituents depend on the density, , and the temperature, T_{i}, and are calculated using Saha equations (e.g. Black & Bodenheimer 1975), with the simplifying assumption that the dissociation of H_{2} is complete before ionization of H begins; and similarly, that the ionization of He is complete before the ionization of He^{+} begins.
Figure 3: The variation of the mean molecular weight with density and temperature. Isopycnic curves are plotted from to , every two orders of magnitude (bottom to top). | |
Open with DEXTER |
If we define
to be the degree of dissociation of hydrogen,
to be the degree of ionization of hydrogen,
to be the degree of single ionisation of helium, and
to be the degree of double ionisation of helium, then the mean molecular weight is given by
(35) |
For densities up to
the ideal gas approximation holds, and hence the gas pressure is
(36) |
Figure 4: The variation of the specific internal energy with density and temperature. Isopycnic curves are plotted from , every two orders of magnitude (from top to bottom). The rotational degrees of freedom of H_{2} are excited around , H_{2} starts to be dissociated around , and H starts to be ionised around . | |
Open with DEXTER |
The specific internal energy (energy per unit mass) of an SPH particle i is the sum of contributions from molecular, atomic and ionised hydrogen, atomic, singly ionised and doubly ionised helium, and the associated dissociation and ionisation energies,
(37) |
= | (38) | ||
= | (39) | ||
= | (40) | ||
= | (41) | ||
= | (42) | ||
= | (43) | ||
= | (44) |
(45) |
In the present work we do not distinguish between the Rosseland-mean and Planck-mean opacities; we use the parametrisation proposed by Bell & Lin (1994) for both, i.e.
The opacity at low temperatures is dominated by icy dust grains. At K the ices evaporate and the opacity is due to metal grains up to K, when the metal grains start to evaporate. The opacity drops considerably in the temperature range from K to K, as it is too hot for dust to exist, and too cool for H^{-} to contribute, so the opacity is mainly due to molecules; this region of low opacity is sometimes referred to as the opacity gap. The opacity starts to increase again above K due to H^{-} absorption and then decreases again above K, when free-free transitions take over. At very high temperatures, electron scattering delivers an approximately constant opacity.
At low temperatures, , the Bell & Lin parametrisation agrees well with the Rosseland-mean dust opacity calculated by Preibisch et al. (1993). Similarly, at high temperatures, , it agrees well with the Rosseland-mean gas opacities calculated by Alexander & Ferguson (1994) and Iglesias & Rogers (1996).
Equation (46) gives local Rosseland- and Planck-mean opacities. To calculate the pseudo-mean opacity used in Eq. (24), we have to convolve this opacity with polytropic density and temperature profiles according to Eq. (23). In Fig. 6 we present the pseudo-mean opacity computed in this way, using a polytropic index n=2. We reiterate that the choice of n affects the computed pseudo-mean opacities only weakly.
Table 1: Opacity law parameters (from Bell & Lin 1994).
Figure 5: The variation of the local Rosseland-mean opacity with density and temperature. Isopycnic curves are plotted from to , every two orders of magnitude (from bottom to top). The opacity gap is evident at temperatures K, over a wide range of densities. | |
Open with DEXTER |
The first test of our new method for treating the energy equation in SPH is to simulate the collapse of a molecular cloud, which initially is spherical with radius and uniform density . We set the background radiation temperature to . This problem has been investigated by Masunaga & Inutsuka (2000) using a code that treats the hydrodynamics in 1 dimension (i.e. assuming spherical symmetry) and the radiative transfer exactly in 3 dimensions (i.e. by solving the angle-dependent and frequency dependent radiation transfer equation). It therefore constitutes a stiff test for our new method to reproduce their results. For the simulation presented here we use SPH particles.
The evolution of the cloud is followed up to density and temperature K, i.e. a difference from the initial conditions of about 16 orders of magnitude in density, and 3 orders of magnitude in temperature.
As long as the density is below
,
the temperature increases slowly with increasing density. For
we can approximate this with
(47) |
The cloud core starts heating more rapidly when it becomes optically thick. This continues until the the temperature reaches
at density
,
when the rotational degrees of freedom of H_{2} start to get excited, and hence the temperature increases more slowly with density (see the kink around
in Fig. 7). As the temperature increases further, the thermal pressure starts to decelerate the contraction, and the first core is formed (Larson 1969;
Masunaga et al. 1998; Masunaga & Inutsuka 2000;
Whitehouse & Bate 2006) at
,
where
is the free-fall time at the start of collapse, i.e.
yr (see Fig. 8).
Figure 6: The variation with density and temperature of the pseudo-mean opacity. Isopycnic curves are plotted as in Fig. 5. For comparison the local opacity at density is also plotted (dashed line). | |
Open with DEXTER |
The first core grows in mass, contracts, and heats up until the temperature reaches , when H_{2} starts to dissociate. Consequently the compressional energy delivered by contraction does not all go to heat the core; instead some of it goes into dissociating H_{2}, and the second collapse starts. This second collapse proceeds until almost all of the molecular hydrogen at the centre has been dissociated. When the density reaches and the temperature rises above , the collapse again decelerates and the second core (i.e. the protostar) is formed (Larson 1969; Masunaga & Inutsuka 2000) at . At first, the second core pulsates (see Fig. 8 and Larson 1969), but eventually it settles down into quasistatic contraction. Due to computational constraints, the evolution is not followed further.
The evolution of the core in our simulation is very similar to that obtained by Masunaga & Inutsuka (2000), as shown in Fig. 7. Differences at densities are attributable to the different opacities we use. The timescales in our simulation are also very similar to those obtained by Masunaga & Inutsuka (2000), as is shown in Fig. 8, where the evolution of the density at the centre of the cloud is plotted against time. The times computed from our simulation fit well with the times from the Masunaga & Inutsuka (2000) simulation if we synchronise the two simulations at density , to avoid discrepancies due to small differences in the initial conditions at the onset of the collapse.
In order to describe the evolution away from the centre of the cloud, and to make a more detailed comparison with the results of the Masunaga & Inutsuka (2000) simulation, we focus on eight representative instants during the cloud evolution. Critical parameters at these instants are listed in Table 2. The central densities have been chosen so as to match those used by Masunaga & Inutsuka (2000) for the same purpose (see their Fig. 1 and their Table 1).
Figure 7: Evolution of the central density and central temperature of the collapsing cloud (thick red line). For comparison, the dashed black line shows the Masunaga & Inutsuka (2000) simulation, and the dotted black lines delineate the different opacity regimes (see Fig. 5). The thin solid lines are the loci for different degrees of H_{2} dissociation (bottom set of blue lines, respectively) and of H ionisation (top green line, x=0.01). The results of our model are very close to the simulation of Masunaga & Inutsuka (2000). Differences at high densities are attributable to our using different opacities. | |
Open with DEXTER |
Figure 8: Evolution of the central density with time, given in units of the initial free fall time . The times of the formation of the first and second core are also marked. The red squares correspond to the Masunaga & Inutsuka (2000) simulation. | |
Open with DEXTER |
Figure 9 shows the run of temperature against density at the instants defined in Table 2. During the early stages of the collapse the variation of temperature with density mimics the evolution of the central temperature and density (see the green points in Fig. 9; points representing previous instants are overlapped). At later stages, the regions around the centre start heating at lower densities, as in the simulation of Whitehouse & Bate (2006).
Table 2: Colour code, elapsed time (t, measured from the beginning of the simulation and in units of the initial freefall time), central density ( ), and central temperature ( ), for the instants illustrated in Figs. 9 to 12.
In Figs. 10 through 12 we present the density, temperature and radial infall velocity profiles at different instants during the evolution of the cloud. These profiles are very similar to those reported by Masunaga & Inutsuka (2000). The velocity profiles (Fig. 12) clearly show the formation of an accretion shock at the boundary of the first core at radius . There is also an accretion shock at the boundary of the second core, initially at a radius of , but later expanding to (cf. Larson 1969).
Figure 9: The run of temperature against density at different instants during the cloud evolution (instants 2 to 8 in Table 2). The region around the centre of the cloud heats at lower densities than the centre of the cloud. The density increases as time evolves (black, green, cyan, magenta, black, green, red). For reference, we also plot the evolution of the temperature and density at the centre of the cloud in our simulation (lower solid black line) and in the Masunaga & Inutsuka (2000) simulation (dashed black line). | |
Open with DEXTER |
We repeat this simulation using different numbers of SPH particles, , to check for convergence. We use (Fig. 13). The run of central temperature against central density is almost identical for different numbers of particles, and the results are fully converged up to densities with . Currently available supercomputing facilities allow SPH simulations of star formation with up to particles, and so the convergence condition quoted above can easily be met.
Figure 10: Density profiles at different instants during the cloud evolution (instants 2 to 8 in Table 2). The colour coding is the same as in Fig. 9. | |
Open with DEXTER |
Figure 11: Temperature profiles at different instants during the cloud evolution (instants 2 to 8 in Table 2). The colour coding is the same as in Fig. 9. | |
Open with DEXTER |
Figure 12: Radial infall velocity profiles at different instants during the cloud evolution (instants 2 to 8 in Table 2). The colour coding is the same as in Fig. 9. | |
Open with DEXTER |
Figure 13: The run of central temperature against central density during the evolution of a molecular cloud, simulated using different numbers of SPH particles: . The results converge for . The black dashed line represents the Masunaga & Inutsuka (2000) simulation. | |
Open with DEXTER |
The second test of our new method is to simulate the evolution of a cloud with uniform initial density , uniform initial temperature , and initial radius , as originally investigated by Boss & Myhill (1992). This problem has recently been revisited by Whitehouse & Bate (2006), using SPH with flux-limited diffusion. In our simulation we have SPH particles. The evolution of the cloud is very similar to the evolution already described in Sect. 4. Figure 14 compares the run of central temperature against central density which we obtain, with that obtained by Whitehouse & Bate (2006). There are three small differences. (i) In our simulation, the cloud starts heating before it becomes optically thick (as reported also by Masunaga & Inutsuka (2000) in a similar test). In contrast, the Whitehouse & Bate (2006) cloud remains strictly isothermal during this phase. (ii) In the Whitehouse & Bate (2006) simulation the centre of the cloud heats up more rapidly at densities , i.e. at lower densities than in our simulation. However, Masunaga & Inutsuka (2000) in a similar test also find lower temperatures than Whitehouse & Bate (2006) in this regime. (iii) There are differences at densities , which are attributable to the use of different opacities. These differences are small and the overall evolution of the cloud is similar in the two simulations.
The third test of our new method is to simulate the evolution of a rotating cloud; the cloud initially rotates as a solid body with angular velocity , and so the ratio of rotational-to-gravitational energy is ; the density includes an m=2 perturbation, i.e. , where is the azimuthal angle in the plane perpendicular to the axis of rotation; the initial radius of the cloud is and the temperature is initially uniform at ; this problem was originally investigated by Boss & Bodenheimer (1979), and has recently been revisited by Whitehouse & Bate (2006). In our simulation we have SPH particles. Figure 15 compares the evolution of the density and temperature of the densest part of the cloud in our simulation, with that obtained by Whitehouse & Bate (2006). The differences are again only small.
The result of the collapse is a binary with separation . In Fig. 16 we plot the density and the temperature on the xy-plane at three instants during the evolution. The components of the binary system are connected by a bar which subsequently fragments. If the collapse were isothermal, this bar should not show any tendency to fragment (e.g. Truelove et al. 1998; Klein et al. 1999; Kitsionas & Whitworth 2002). However, Bate & Burkert (1997) show that if the gas is allowed to heat, but the heating happens at sufficiently high densities, , then the bar does fragment. In our simulation, the gas in the bar starts to heat up rapidly only when the density reaches . Therefore the fragmentation of the bar is consistent with the predictions of Bate & Burkert (1997). However, we note that Whitehouse & Bate (2006) do not report any bar fragmentation.
Figure 14: Evolution of the central density and central temperature for the Boss & Myhill (1992) test problem. The dashed line corresponds to the Whitehouse & Bate (2006) simulation, and the dotted lines define the different opacity regimes (see Fig. 5). The results of our model are very close to the simulation of Whitehouse & Bate (2006). Differences at densities are attributable to the use of different opacities. | |
Open with DEXTER |
Figure 15: Evolution of the density and temperature at the densest part of a collapsing, rotating molecular cloud. The dashed line corresponds to the Whitehouse & Bate (2006) simulation of the same problem. The results of our simulation are very close to those of Whitehouse & Bate (2006). Differences at densities are attributable to different opacities. | |
Open with DEXTER |
Figure 16: Three instants from our simulation of the Boss & Bodenheimer test problem, at , , and ( left to right). We plot the logarithmic density ( top) and the logarithmic temperature ( bottom), on the xy-plane, i.e. the plane perpendicular to the rotation axis. The central densities and central temperatures of the southern condensation are , K, , K, , K. The corresponding central densities and central temperatures of the northern condensation are , K, , K, , K. | |
Open with DEXTER |
Finally, we test the time-dependence of our new method by simulating the relaxation of temperature fluctuations in a static sphere with uniform density
and radius
.
We assume an equilibrium temperature of
and an initial temperature perturbation of the form
(Masunaga et al. 1998; Spiegel 1957), where
is the amplitude of the perturbation and
is its characteristic wavenumber. Masunaga et al. (1998) have shown that at subsequent times the temperature should be
(50) |
In Figs. 17 to 20 we present the radial temperature profiles at different instants during the relaxation towards equilibrium for spheres having different optical depths, . The simulation results approximate well to Eq. (48).
The relaxation rates are also reproduced well. In Fig. 21 we present the dispersion relation, i.e. the relaxation rate for different values of the ratio . We plot the relaxation rates for all the SPH particles that represent the uniform density sphere, calculated at seven different snapshots during the temperature relaxation (dots that saturate to form lines). The theoretical values (calculated from Eq. (49)) are also plotted (red line and squares).
We have developed a new method to treat the influence of radiative transfer on the energy equation in SPH simulations of star formation. The method uses the density, temperature and gravitational potential of each particle to make an educated estimate of the mean optical depth which regulates its heating and cooling. It can treat both the optically thin and optically thick regimes.
The energy equation takes account of heating by compression or cooling by expansion (i.e.
work); viscous dissipation; external irradiation; and radiative cooling. In situations where the thermal timescale is much longer than the dynamical timescale, the resulting thermal inertia effects are captured properly. Conversely, where the thermal timescale is much shorter than the dynamical timescale, we avoid very short timesteps, essentially by assuming thermal equilibrium.
Figure 17: Thermal relaxation of a static, uniform sphere of optical depth ; seven instants are plotted, showing the temperature relaxing to its equilibrium value . | |
Open with DEXTER |
Figure 18: Same as Fig. 17, but for a sphere of optical depth . | |
Open with DEXTER |
Figure 19: Same as Fig. 17 but for a sphere of optical depth . | |
Open with DEXTER |
Figure 20: Same as Fig. 17 but for a sphere of optical depth . | |
Open with DEXTER |
The equation of state and the internal energy take account of (i) the rotational and vibrational degrees of freedom of H_{2}, and (ii) the different chemical states of hydrogen and helium (cf. Black & Bodenheimer 1975; Boley et al. 2007)
A simple parametrisation of the frequency averaged opacity is used (Bell & Lin 1994), which reproduces the basic features of more sophisticated opacity models (e.g. Preibisch et al. 1993; Alexander & Ferguson 1994). This parametrisation accounts for the effects of dust sublimation, molecules, H^{-}, free-free transitions, and electron scattering
To test the SPH-RT method we have examined the collapse of a molecular cloud of initially uniform density and uniform temperature. The collapse proceeds almost isothermally until the density in the centre rises above . Then the temperature rises rapidly, the thermal pressure decelerates the collapse, and the first core is formed. The first core grows in mass, and contracts quasistatically. When its temperature reaches , the H_{2} starts dissociating and the second collapse starts, resulting in the formation of the second core, i.e. the protostar. Our method reproduces well the results of the detailed simulation of Masunaga & Inutsuka (2000): the first and the second cores form at similar densities, having similar sizes, and at similar times after the start of the collapse.
We have also performed the Boss & Myhill (1992) and Boss & Bodenheimer (1979) tests, and obtained results very similar to those of Whitehouse & Bate (2006). Finally, we have performed the thermal relaxation test of Masunaga et al. (1998). The geometries treated in this paper establish the fidelity of the method in treating both spherical and flattened geometries. Furthermore, the method also performs well on the Hubeny (1990) test, which deals with equilibrium discs, and hence it can also be applied to disc simulations. We will discuss the Hubeny (1990) test, and applications of this method to discs in a forthcoming paper (Stamatellos & Whitworth 2007).
The new SPH-RT method performs very well, and most importantly it is very efficient. The computational time is almost the same as (only 3% longer than) an SPH simulation using a barotropic equation of state. The method is inherently three dimensional, and so it can be used to treat a variety of astrophysical systems, where the radiative processes and thermal inertia effects are important. We will report on applications of the method in future publications.
Figure 21: Dispersion relation for the thermal relaxation mode. The relaxation rates (in units of ) of all SPH particles at seven different instants are plotted against the ratio (dots that saturate to form lines). The theoretical values (calculated from Eq. (49)) are also plotted (red line and squares). | |
Open with DEXTER |
Acknowledgements
We would like to thank the referee, S. Inutsuka, for his suggestions which helped to improve the original manuscript. The computations reported here were performed using the UK Astrophysical Fluids Facility (UKAFF). We thank S. Whitehouse and M. Bate for kindly providing data from Whitehouse & Bate (2006), C. Clarke, M. Bate, I. Bonnell and R. Wünsch for useful discussions, and D. Price for the use of SPLASH (Price 2007). We also acknowledge support by PPARC grant PPA/G/O/2002/00497.