text:
# Monte Carlo radiative transfer in molecular cloud cores

## 1 Introduction

Understanding the structure of pre–protostellar cores is an essential step towards an understanding of protostellar evolution. The density distribution immediately prior to the onset of gravitational collapse defines the initial conditions for collapse and hence one has a strong motivation to attempt to derive this density distribution from the observations. Since these objects are most commonly observed by means of their millimeter dust emission, one has a strong interest in understanding the temperature distribution of the dust in the pre–protostellar core and its influence on the emergent intensity distribution at millimeter-submillimeter wavelengths. In this article, we present a radiative transport code which has the interpretation of mm-submm maps of pre–protostellar cores as its main goal.

The approach which we have adopted is influenced by the fact that, as a rule, the density structure of these objects shows evidence for large deviations from spherical symmetry. More precisely, the maps of the millimeter dust emission show clear departures from circular symmetry (see e.g. André et al. 2000; Caselli et al. 2002a,b; Tafalla et al. 2002). In some cases there is also evidence for polarization and hence for a magnetic field (Ward-Thompson et al. 2000). It seems likely to us that such behavior is caused by magnetic fields of energy densities sufficiently large to influence the core structure (see also Shu et al. 1987) and cause flattening along field lines (Basu 2000, Jones & Basu 2002). We can attempt to infer the density distribution on the basis of mm-submm continuum maps if we can derive the temperature distribution in regions where departures from spherical symmetry are important. A first attempt in this direction was made by Zucconi et al. (2001, hereafter ZWG, see also Evans et al. 2001 for the spherically symmetric case) but these authors neglected heating of dust grains due to re-absorption of photons emitted by the cloud itself. Here we present a Monte–Carlo code without this limitation but which is capable of handling a variety of geometries. We expect that this type of analysis will be especially useful as a tool for the interpretation of the high class data we expect to come from future instruments like ALMA and HERSCHEL.

Some high quality data are already available and are consistent with a picture in which the dust temperature of pre-protostellar cores decreases towards the center (Ward-Thompson et al. 2002, Bianchi et al. 2003). ISO observations in the mid-IR (Bacmann et al. 2000) and in the far-IR (Ward-Thompson et al. 2002) suggest dust temperature gradients consistent with heating from the external interstellar radiation field (ISRF). The available models predict a factor of $`2`$ increase in temperature from center to edge (ZWG, Evans et al. 2001, André et al. 2003), with the gradient dependent on the cloud structure. The predicted emission of several cores is in reasonable accord with the data, and implies that the assumed equilibrium configurations are either unstable or maintained by a magnetic field.

Additional observational constraints on the thermal structure of prestellar cores are given by spatially resolved measurements of the gas temperature from NH3 observations. Low-mass cores show fairly uniform gas temperature (e.g. Tafalla et al. 2002 for L1517B and L1498), whereas massive quiescent cores in Orion show significant temperature drops from edge to center (Li et al. 2003). Thus another check on our understanding of core structure is possible if radiative transfer models are also able to predict the gas temperature distribution resulting from the balance of the relevant heating and cooling mechanisms. We therefore have also developed the capability to predict gas temperature distributions for our model cores.

While techniques have been developed for the study of radiative transfer in cold cores in one dimension (e.g. Leung 1976, Rowan-Robinson 1980, Ivezic et al. 1997), their extension to three-dimensional geometries is not straightforward. A Monte Carlo technique has the advantage that it deals easily with any general geometry, independently of the density distribution and the anisotropy of the incident radiation field. The increasingly faster workstations available render this technique very useful, and fully 3-D Monte Carlo radiative transfer codes have been developed recently (e.g. Wolf et al. 1999, Niccolini et al. 2003). However, in these latter studies the radiative sources are internal, and this is a computationally distinct problem from a core heated from the outside for two main reasons. First, as we are interested in cold cores, it is reasonable to expect that they will be optically thin to their own emission, and therefore, for methods that use iteration, the temperature distribution quickly relaxes to the final value; alternatively, for methods that follow individual packets until they exit the domain, each individual packet will interact with the dust often only once. Second, as packets are launched from the outside, a wide range of spatial scales demands a prohibitively large photon statistics, as the probability of a packet being launched in the direction of the innermost cells becomes increasingly small.

In general, our paper confirms and supports the results recently obtained for embedded prestellar cores by Stamatellos & Whitworth (2003) with an independent but similar Monte Carlo radiative transfer code. While Stamatellos & Whitworth (2003) focus on the important effect of the ambient medium in which (spherical) cores are embedded, our paper mostly emphasizes the consequences of deviations from spherical symmetry in the density distribution and anisotropies of the ISRF. The two papers represent therefore complementary attempts to model radiative transfer in molecular cloud cores in less idealized situations than those considered so far.

The structure of this paper is as follows. In Sect. 2, we set the general problem of radiative transfer in a dusty cloud. In Sect. 3, we describe the Monte Carlo method we have used. In Sect. 4, we present our choice of opacities and interstellar radiation field and briefly discuss their effect in the results. In esections 5 and 6 we present our results for two- and three-dimensional models, respectively. Finally, in Sect. 7 we summarize our conclusions.

## 2 Radiative transfer in a dusty cloud

In local thermodynamic equilibrium and ignoring the effects of scattering of radiation by dust grains, the radiation field in a dusty cloud satisfies the equation of transfer

$$\widehat{𝐧}I_\nu (𝐫,\widehat{𝐧})=\rho (𝐫)\kappa _\nu (𝐫)\{B_\nu [T_\mathrm{d}(𝐫)]I_\nu (𝐫,\widehat{𝐧})\},$$
(1)
where $`I_\nu (𝐫,\widehat{𝐧})`$ is the specific intensity of radiation (in Jy ster-1) at position $`𝐫`$ in the direction $`\widehat{𝐧}`$ at frequency $`\nu `$, $`\rho `$ is the total (gas plus dust) density, $`\kappa _\nu `$ the total absorption opacity, $`T_\mathrm{d}`$ is the dust temperature, and $`B_\nu `$ is the Planck function,

$$B_\nu (T_\mathrm{d})=\frac{2h\nu ^3}{c^2}\left[\mathrm{exp}\left(\frac{h\nu }{kT_\mathrm{d}}\right)1\right]^1.$$
(2)
The dust temperature distribution $`T_\mathrm{d}(𝐫)`$ is obtained by solving the equation of balance of emitted and absorbed radiation,

$$_0^{\mathrm{}}\kappa _\nu B_\nu [T_\mathrm{d}(𝐫)]𝑑\nu =_0^{\mathrm{}}\kappa _\nu J_\nu (𝐫)𝑑\nu ,$$
(3)
where

$$J_\nu (𝐫)=\frac{1}{4\pi }I_\nu (𝐫,\widehat{𝐧})𝑑\mathrm{\Omega }$$
(4)
is the mean intensity of radiation (in Jy).

ZWG made the symplifying assumption that the dust grains are only heated by the incident interstellar radiation field $`J_\nu ^{\mathrm{ISRF}}`$ (assumed isotropic), neglecting re-emission of radiation absorbed by the grains themselves. That is to say, they assumed that the cloud is optically thin to its own radiation. Under this hypothesis, the equation of radiative transfer has the solution

$$J_\nu (𝐫)=\frac{J_\nu ^{\mathrm{ISRF}}}{4\pi }\mathrm{exp}[\tau _\nu (𝐫,\widehat{𝐧})]𝑑\mathrm{\Omega },$$
(5)
where

$$\tau _\nu (𝐫,\widehat{𝐧})=_r^{\widehat{𝐧}R}\rho (\widehat{𝐫}r^{})\kappa _\nu (\widehat{𝐫}r^{})𝑑r^{},$$
(6)
is the optical depth from a point $`𝐫`$ inside the cloud to a point $`\widehat{𝐧}R`$ on the cloud boundary, in the direction $`\widehat{𝐧}`$. The calculation of the radiation field inside the cloud is thus reduced to a large number of integrations along different directions (rays) defined at each grid point. This simplification allowed ZWG to model axisymmetric non-spherical density distributions, but has the disadvantage of slightly underestimating the dust temperature, particularly at the cold center of the cloud.

Our general approach is the following. First, for a given ISRF we compute the radiation field at any point inside the cloud core and solve the equation of balance of radiation emitted and absorbed by dust grains to obtain the dust temperature distribution $`T_\mathrm{d}(𝐫)`$. Second, we compute the flux emitted in a given view direction at different wavelengths and we compare the results with observed spectral energy distributions and monochromatic maps for specific objects.

In several cases we have also considered the effect of a spherical envelope surrounding the cloud core on the resulting temperature distribution. The rationale for this is that these objects are always embedded in a photon–dominated region (PDR) which transforms the incident optical–UV field into mid and far IR radiation. Often in practice, this is radiation from transiently heated particles, either small grains or polycyclic aromatic hydrocarbons, emitting in the 3–30 $`\mu `$m range. Since we wish to focus on the temperature structure of the high density core interior, we have decided to assume that the radiation field incident on embedded cores is cut off in the optical-UV wavelength range and that the radiation absorbed in the PDR is reradiated in the infrared. We thus incorporate the photons reradiated by the PDR in the external field. The error involved here can be estimated by varying the form of the ISRF and appears to be small (though see the discussion of André et al. 2003). Our focus here is on the effects of the core geometry upon temperature structure and we place minor emphasis on the spectral energy distribution of the incident field.

We stress that our approach is general: the density distribution, the optical properties of the dust grains, the intensity and spectral energy of the radiation source(s) are input parameters that can be freely specified and can be adapted to model various astrophysical situations.

## 3 The Monte Carlo method

In the Monte Carlo technique, the computational domain is divided in a large number of cells of mass $`m_i`$ that absorb and emit radiation. The energy that enters the computational domain is divided in $`N_\gamma `$ monochromatic packets of equal energy, that are launched stochastically from the boundary and followed until they exit the cloud. These packets may be absorbed by the dust particles, resulting in a temperature increase of the cells where they are absorbed, and are immediately reemitted to enforce radiative equilibrium. If the cell $`i`$ absorbs $`N_i`$ packets, the energy absorbed by the cell per unit time is

$$\frac{dE_i^{\mathrm{abs}}}{dt}=\left(\frac{N_i}{N_\gamma }\right)L^{\mathrm{ISRF}},$$
(7)
where $`L^{\mathrm{ISRF}}`$ is the luminosity of the ISRF at the cloud’s surface, obtained integrating the mean intensity of the (anisotropic) ISRF over the cloud’s surface and over frequency,

$$L^{\mathrm{ISRF}}=_0^{\mathrm{}}𝑑\nu 𝑑\mathrm{\Omega }J_\nu ^{\mathrm{ISRF}}(\widehat{𝐧})\widehat{𝐧}𝑑𝐒.$$
(8)
Assuming LTE conditions, the energy emitted by the $`i`$-th cell per unit time is

$$\frac{dE_i^{\mathrm{em}}}{dt}=4\pi m_i_0^{\mathrm{}}\kappa _\nu B_\nu (T_{\mathrm{d},i})𝑑\nu $$
(9)
where $`T_{\mathrm{d},i}`$ is the dust temperature of the $`i`$-th cell, obtained by equating the absorbed and emitted energies per unit time,

$$_0^{\mathrm{}}\kappa _\nu B_\nu (T_{\mathrm{d},i})𝑑\nu =\left(\frac{N_i}{N_\gamma }\right)\frac{L^{\mathrm{ISRF}}}{4\pi m_i}.$$
(10)
The integral on the left-hand side of this equation can be tabulated for a grid of values of the dust temperature and opacity. The solution of Eq. (10) can then be easily obtained by iteration and interpolation. In the limit to the continuum, $`N_\gamma \mathrm{}`$ and $`m_i0`$, Eq. (10) is equivalent to Eq. (3), but differs from the approximation adopted in ZWG, Eq. (5), because in $`N_i`$ are included the packets that may have been absorbed more than once in the cloud, and subsequently re-emitted.

We handled re-radiation (and, therefore, energy conservation) using the method devised by Bjorkman & Wood (2001): once a packet is absorbed, it is immediately re-emitted (to conserve energy) at a new frequency determined by the local dust temperature. To achieve this, one starts by noting that if a cell has emissivity $`\kappa _\nu B_\nu (T_{\mathrm{d},i}\mathrm{\Delta }T_{\mathrm{d},i})`$ prior to absorbing a wave packet, then after packet absorption its temperature increases by an amount $`\mathrm{\Delta }T_{\mathrm{d},i}`$, and the cell emissivity becomes $`\kappa _\nu B_\nu (T_{\mathrm{d},i})`$. Thus, the increment in the cell emissivity is

$`\mathrm{\Delta }j_\nu `$ $`=`$ $`\kappa _\nu [B_\nu (T_{\mathrm{d},i})B_\nu (T_{\mathrm{d},i}\mathrm{\Delta }T_{\mathrm{d},i})]`$ (11)
$`\kappa _\nu \left({\displaystyle \frac{dB_\nu }{dT_\mathrm{d}}}\right)_{T=T_{\mathrm{d},i}}\mathrm{\Delta }T_{\mathrm{d},i}.`$
A photon packet is immediately reemitted with a frequency obtained from a probability distribution having the same spectral shape as $`\mathrm{\Delta }j_\nu `$,

$$\frac{dP_i}{d\nu }=\frac{\kappa _\nu }{K}\left(\frac{dB_\nu }{dT_\mathrm{d}}\right)_{T=T_{\mathrm{d},i}},$$
(12)
where $`K=_0^{\mathrm{}}\kappa _\nu (dB_\nu /dT_\mathrm{d})𝑑\nu `$ is a normalization constant. The re-emitted photon packet is then followed until a new interaction occurs, and the procedure is repeated until the photon leaves the cloud. This method ensures energy conservation in a statistical sense. Due to the low temperature of the cores, energy is generally re-radiated at long wavelengths, to which the core has a low optical depth, and the overall effect of re-emission on the dust temperature is more significant in the core’s cold interior but is generally very small.

### 3.1 Tests of the code

We have tested our Monte Carlo code in different ways. Here we present the results of two complementary tests, one checking the accuracy of the packet propagation procedure and one checking the equality of emitted and absorbed energy.

To test the packet propagation procedure, we have computed dust temperatures profiles for spherically symmetric clouds, for which results can also be obtained with the ray integration method of ZWG. Specifically we have assumed the density profile of a Bonnor-Ebert (BE) sphere, an isothermal, pressure-bounded, equilibrium model, that reproduces the centrally flattened density distribution and the rapid density decrease at large radii observed in most cloud cores. Fig. 1 shows the dust temperature as function of radius for a BE with radius $`R=0.1`$ pc, central density $`n(\mathrm{H}_2)=4.4\times 10^6`$ cm-3, and isothermal gas temperature $`T_\mathrm{g}=12`$ K (note that the gas temperature here is merely for the purpose of specifying the assumed density distribution). As shown in Fig. 1, the result obtained using the ray integration method (solid line), is identical to the Monte Carlo result using the same approximation of optically thin re-emission (dotted line). We also notice that when re-emission is taken into consideration (dashed line), there is only a small increase of the central temperature, even in this case, when the extinction through the core center is $`A_{2.2\mu \mathrm{m}}25`$, and the central dust temperature is little more than 6 K.

To test energy conservation, we have computed the energy emitted per unit time by an arbitrary area element on the cloud’s surface by integrating $`I_\nu `$ (the solution of eq. 1) and $`I_\nu ^{\mathrm{ISRF}}`$ over solid angles and frequencies, finding a very small discrepancy between the two values (less than 1 %). This test shows that for the given number and size of cells, the accuracy of Eq. (10) and (12), that contain approximations introduced by the discretization of the computational volume, is very good.

## 4 Dust opacities and the interstellar radiation field

### 4.1 The interstellar radiation field

Following ZWG and Evans et al. (2001), we adopt as a reference the ISRF given by Black (1994), which is an average for the solar neighbourhood. Scaling the intensity of the ISRF by a factor $`G_0`$ leads to a change in the dust temperature at all radii by a factor of $`G_0^{1/(4+\beta )}`$, where $`\beta 2`$ is the power law index for the long wavelength opacity. Actually the radiation incident on molecular cloud cores may differ not only in intensity but also in spectral shape, depending on the degree of embedding of the cores in the ambient cloud and on the possible vicinity to hot stars. It appears for example to be the case that many cores associated with the $`\rho `$ Oph star forming regions are subjected to an incident ISRF roughly an order of magnitude larger than the solar neighbourhood ISRF. An example of what can happen is given by a recent study by André et al. (2003) of a particular core, heated by a ISRF higher than the Black field by an order of magnitude in both the mid-IR and far-IR. This resulted in a higher central temperature and a sharper gradient at the edge than in our models for a similar density distribution. Thus observers attempting to infer density distributions from maps of the mm–submm dust emission must bear in mind the fact that both the radiation field and the geometry can play a role.

In this paper however, our objective is to study the effects of geometry rather than those of the spectral characteristics or intensity of the incident radiation field. Geometry in this context can mean the geometry of the density distribution or that of the radiation field. We thus also consider a case where we use our Monte Carlo model to consider an anisotropic incident ISRF as exists in many galactic reflection nebulae.

### 4.2 Dust opacities

We adopt the dust opacities tabulated by Ossenkopf & Henning (1994, hereafter OH) in the wavelength range 1 $`\mu `$m–1.3 mm. Recent work by Bianchi et al. (2003) and by Kramer et al. (2003) has shown that these theoretically derived opacities are consistent with the observed ratio of mm dust continuum intensity and near infrared extinction. It is also clear from the study by OH that one can expect the dust opacity to vary between the lower density outer layers of a core and the high density interior where most of the heavy element content of the core is in the form of ice of various sorts. We neglect such effects in this study mainly because, as shown below, the variation between different OH dust models is relatively mild.

We label the opacity according to the column number in Table 1 of OH: OH1, standard MRN distribution; OH4, MRN with thin ice mantles; OH5, MRN with thin ice mantles, after $`10^5`$ yr of coagulation at density $`10^6`$ cm-3; OH8, MRN with thick ice mantles, after $`10^5`$ yr of coagulation at density $`10^6`$ cm-3.

The OH opacities are tabulated for wavelengths between 1 $`\mu `$m and 1.3 mm. Within this range we have obtained values of the opacity at arbitrary wavelengths by four-point interpolation and above 1.3 mm, we have used a power law extrapolation based on the last two tabulated values. We neglect scattering completely as we are interested in the transport of infrared photons.

The effect of the various opacities on the temperature profile of a cloud is shown in Fig. 2, for the reference case of a BE sphere with central density $`n(\mathrm{H}_2)=4.4\times 10^6`$ cm-3, radius 0.1 pc, gas temperature $`T_\mathrm{g}=12`$ K, embedded in a spherical envelope with total extinction $`A_{2.2\mu \mathrm{m}}^{\mathrm{env}}=0.1`$. In this particular example, the gas-to-dust ratio was varied for each opacity prescription in order to keep fixed the total extinction at 2.2 $`\mu `$m.

As Fig. 2 shows, temperature variations due to different choices of the opacity are generally small. The opacities OH1, OH4 and OH5 have a similar dependence on wavelength, and the resulting temperature difference is equivalent to scaling the intensity of the ISRF by about $`10\%`$. The OH8 opacity, having a slightly different dependence on wavelength, results in a the temperature profile slightly different from the other three cases. In all remaining models presented in this study we have used OH5.

## 5 Two-dimensional models

In this section we consider two particular situations that require a two-dimensional modeling of the radiative transfer. In the first example, the external radiation field is isotropic but the density distribution is not spherically symmetric. In the second example, the density profile is spherically symmetric but the incident radiation field is anisotropic, as the cloud is heated by an external stellar source. Both situations are easily handled by our Monte Carlo code.

### 5.1 Singular isothermal toroids

As a first application of our Monte Carlo code, we consider the radiative transfer in a singular isothermal toroid (Li & Shu 1996), a scale-free, axisymmetric equilibrium configuration of an isothermal cloud under the influence of sef-gravity, gas pressure and magnetic forces. These toroids are characterized by a single parameter, $`H_0`$, representing the fractional amount of support provided by the magnetic field. We choose the particular value $`H_0=0.5`$ to model a cloud with moderate axial ratio ($`2.5`$) intermediate between a spherical unmagnetized cloud and a magnetically dominated disklike configuration.

It should be noticed that, in all models presented, we are simply computing the dust temperature of pre-defined core density distributions, and therefore the resulting core is not the solution of the equilibrium equation if gas-dust coupling is assumed. Also, the practical definition of the boundary as an isobaric surface is only possible for spherical clouds, as in this case the isopycnic and isobaric surfaces necessarily coincide. In all other geometries, we take the boundary as the isopycnic surface corresponding to a specified value of the density.

In Fig. 3 are shown both the density distribution of the toroid and the resulting temperature distribution. The boundary of the configuration is defined as the isopycnic surface with $`n(\mathrm{H}_2)=10^4`$ cm-3. We choose this value of the density because it represents the threshold value that characterizes typical dense (“ammonia”) cores. For a density profile characterized by a $`r^2`$ law, the visual extinction of the cloud material with density below this threshold value corresponds to about $`A_V2`$–3. Notice that the temperature varies over the boundary surface by a factor of about two from the origin to the maximum radius, and this difference is due to an increase of grain exposure to the ISRF as the radius increases, an effect known also for disks of low aspect ratio (e.g. Spagna et al. 1991), and is particularly enhanced in a toroid. It is also evident from Fig. 3 that the temperature decreases outward both in the horizontal and vertical directions.

In Fig. 3, we also show synthetic (unconvolved with an observing beam) maps of the expected dust emission at 350 $`\mu `$m and 1.3 mm. In particular at 1.3 mm, the “observed contours” reflect fairly faithfully the projected column density and hence the density distribution. At 350 $`\mu `$m however, the intensity distribution is much more extended than at 1.3 mm. The half-power dimensions are $`1.5^{\prime \prime }\times 3.7^{\prime \prime }`$ at 1.3 mm (in angular size at a distance of 140 parsec) but $`18^{\prime \prime }\times 25^{\prime \prime }`$ at 350 $`\mu `$m. Thus the effect of heating from the exterior is mainly observable as an increase in size with frequency.

### 5.2 Anisotropic Radiation Field

While the above discussion concerns cores of non-spherical geometry subject to an isotropic incident radiation field, it should be noted that one expects often to find cases where the incident radiation field is anisotropic. Most obviously, this is the case in galactic reflection nebulae where the radiation field from single early type stars may be dominant. We therefore here (in a slight parenthesis to the discussion elsewhere in this paper) demonstrate the application of our code to such a case and compute the expected temperature distribution in a core subject to radiation from a nearby B star in addition to the isotropic ISRF (note earlier work in this field by Natta et al. 1981).

The vicinity of a star to a prestellar core results in an anisotropic radiation field, leading to a 2-D radiation transfer problem in the case of a spherically symmetric core. It should be noted that this is a case where the approximation of ZWG (Eq. 5) may not hold, as the expected increase in dust temperature may make the core no longer optically thin to its own radiation.

To study the effect of such a field, we start by considering a system composed of a prestellar core represented as a BE sphere, with the physical parameters given in section 3.1, and a B3 star at a distance of 0.15 pc from the edge of the core. The large range of spatial scales involved makes it difficult for a Monte Carlo program to solve this problem consistently, by emitting packets from the star, and following them through the ISM and the prestellar core. We avoid this problem by assuming that the stellar radiation is reprocessed by a tenuous but spatially extended photodissociation region (PDR) surrounding the core, and use a fit to the spectral emission of the PDR computed by Désert et al. (1990). In this way, the situation is reduced to the case treated in section 3.1, with the difference that the incident radiation is characterized by a PDR spectrum, instead of the standard ISRF, and a variable intensity over the surface of the core. We further assume that the total luminosity over any area element on the cloud’s surface is the same as that obtained assuming no attenuation between the star and envelope. Thus, the luminosity of a surface element in spherical coordinates is

$$dL(\theta )\frac{1}{r^2(\theta )}\mathrm{cos}\theta _\mathrm{i}\mathrm{sin}\theta d\theta ,$$
(13)
where $`\theta `$ is the angle between the point on the surface of the core and the star as seen from the core’s center, $`r(\theta )`$ is the distance of that point to the star, and $`\theta _\mathrm{i}`$ is the angle of incidence of the radiation. The proportionality constant is obtained by setting the total incident luminosity equal to $`(\mathrm{\Omega }/4\pi )L_{}`$, where $`\mathrm{\Omega }`$ is the solid angle of the core as seen from the star.

The resulting temperature distribution is shown in the top panel of Fig. 4. As expected, the cloud is hotter on the side facing the star, and the temperature decreases both radially and azimuthally, with the lowest temperature close to the centre of the core, which is almost twice higher than the minimum temperature of the Black ISRF heated core with the same density profile (Fig. 1). The ratio $`L_{}/L^{\mathrm{ISRF}}`$ is $`20`$, and the inverse square dependence of $`dL`$ (see Eq. 13) that makes this ratio over a surface element be much larger close to the star, has the consequence that the isotropic ISRF is negligible almost everywhere. The exception to this is on the core’s side opposite to the star, where the shielding provided by the cloud’s center suffices to make the isotropic ISRF the dominant heating field.

In Fig. 4 we show the simulated emission maps at 1.3 mm, 350$`\mu `$m, 100$`\mu `$m and 60$`\mu `$m. The 1.3 mm map shows slight distortions due to the temperature structure but clearly allows a reasonable determination of the mass distribution in the BE sphere. As one goes to shorter wavelengths however, the distortions become extreme and at 60 microns, one sees essentially the heated surface of the cloud. Thus maps at different wavelengths may allow the identification of a heating source and show where the star is along the line of sight relative to the cloud. In such situations, there will often be mid–IR data available which will additionally define the geometry of the “envelope” or PDR layer hypothesised in our analysis.

## 6 Three-dimensional models

Reverting to the case of an isotropic radiation field, we now consider the radiative transfer in a fully 3-D density distribution built on the results of Galli et al. (2001) for magnetized disks with uniform mass-to-flux ratio $`M(\mathrm{\Phi })/\mathrm{\Phi }`$ (isopedic). A computationally convenient characteristic of these models is the existence of non-axisymmetric but analytical solutions of the equations of magnetostatics, corresponding to surface density distributions $`\mathrm{\Sigma }`$ separable in polar coordinates $`(\varpi ,\phi )`$:

$$\mathrm{\Sigma }(\varpi ,\phi )=\frac{\mathrm{\Theta }a^2}{2\pi ϵG\varpi }S(\phi ),$$
(14)
where

$$S(\phi )=\frac{\sqrt{1e^2}}{1+e\mathrm{cos}\phi },$$
(15)
with $`0<e<1`$. Here, $`\mathrm{\Theta }1`$ and $`ϵ1`$ represent the increase of gas pressure and the reduction of the gravitational constant due to magnetic pressure and magnetic tension, respectively, related to the nondimensional mass-to-flux ratio $`\lambda =2\pi G^{1/2}M(\mathrm{\Phi })/\mathrm{\Phi }`$ by the relations (Shu & Li 1997)

$$ϵ=1\frac{1}{\lambda ^2},\mathrm{\Theta }=\frac{\lambda ^2+3}{\lambda ^2+1}.$$
(16)
We construct, in a rather ad hoc way, a 3-D model assuming a vertical stratification of gas in hydrostatic equilibrium and imposing the condition that the column density of the resulting density distribution reproduces Eq. (14), obtaining

$$\rho (\varpi ,\phi ,z)=\frac{a^2\mathrm{\Theta }^2}{8\pi ϵ^2G\varpi ^2}S^2(\phi )\mathrm{sech}^2\left[\frac{z}{H(\varpi ,\phi )}\right],$$
(17)
where $`H(\varpi ,\phi )=a^2/\pi G\mathrm{\Sigma }(\varpi ,\phi )`$.

As a specific example, we adopt the parameters derived by Galli et al. (2001) to match the thermal dust emission map obtained by Ward-Thompson et al. (1999) for the starless core L1544. We assume an effective sound speed $`a=0.21`$ km s-1, a nondimensional mass-to-flux ratio $`\lambda =2`$, an eccentricity $`e=0.54`$, and we define the core boundary by the isopycnic surface $`n(\mathrm{H}_2)=10^4`$ cm-3.

For this model we have also computed the gas temperature. To do this, we first notice that the energy deposited in the gas by cosmic ray ionization is negligible when compared to the energy absorbed by the dust, so that the energy transfer between gas and dust will not significantly affect the grain temperature. Therefore, it is only necessary to solve the equation of thermal equilibrium of the gas,

$$\mathrm{\Gamma }_{\mathrm{cr}}\mathrm{\Lambda }_\mathrm{g}\mathrm{\Lambda }_{\mathrm{gd}}=0,$$
(18)
with $`\mathrm{\Gamma }_{\mathrm{cr}}`$ the cosmic-ray heating rate, $`\mathrm{\Lambda }_\mathrm{g}`$ the gas cooling rate by molecular and atomic transitions, and $`\mathrm{\Lambda }_{\mathrm{gd}}`$ the gas-dust energy transfer rate. For these quantities, we assume the values given by Goldsmith (2001). We also assume a CO depletion factor of the form $`f_{\mathrm{dep}}=\mathrm{exp}(n/n_{\mathrm{dep}})`$, where the critical density for CO depletion is taken to be $`n_{\mathrm{dep}}=5.5\times 10^4`$ cm-3. Fig. 5 shows the density, dust temperature and gas temperature of our model, with the left and right panels showing the $`y=0`$ and the $`z=0`$ plane, respectively.

This 3-D model and the 2-D toroidal model presented in the previous section show qualitatively similar temperature distributions, but the mismatch between isopycnic and dust isothermal surfaces is enhanced in the 3-D model, particularly in the $`y=0`$ plane where the geometry is more complex.

The gas temperature, as it depends on both the dust temperature and on the gas density, shows a different distribution from either (Fig. 5). At the dense center, the temperatures of the gas and the dust are well coupled, but they are decoupled at the lower density boundary. This results in a non-monotonic behaviour of the gas temperature, which reaches a maximum at intermediate density (an effect that also occurs in 1D geometry, as shown in Galli et al. 2002). It is also evident, from Fig. (5) that the gas temperature has a slightly higher value in the right side of the singularity than at the left, which is a consequence of the different dust temperature values at the density where the turnover in the gas temperature occurs on the two sides of the singularity.

Fig. 6 shows simulated dust emission maps at 1.3 mm and 350 $`\mu `$m for the side view and top view. The flattening towards the center due to the dust temperature gradient is evident from the comparison between the emission at the two wavelengths. The vertical density stratification has a fairly strong effect on the shape of the map seen from the side, but is not noticable from the top. In general, one notes a slight “bottleneck” deviation from a purely elliptical shape of the singular isothermal disk model present in Galli et al. (2001). We note also that the “cometary shape” of the edge-on 1.3 mm map is rather similar to that observed towards L1544 and L63 by Ward-Thompson et al. (1999).

## 7 Conclusions

In this paper we have presented a Monte Carlo model for radiative transfer applicable to the study of the physical conditions in starless molecular cloud cores. The method described is able to compute accurately and efficiently the dust and gas temperature distribution in clouds of specified but arbitrary density profiles heated by the ISRF, an external stellar source, or a combination thereof. In addition, the method allows the computation of synthetic maps of the emission at near-infrared and submillimeter wavelengths for arbitrary viewing angles, and to predict the emergent spectral energy distribution. We have successfully tested our code by detailed comparison with one-dimensional calculations in spherical symmetry (ZWG, Galli et al. 2002) and have shown that the assumption of effective optical thinness is a good one for prestellar cores. As preliminary applications of our code, we have considered the radiative transfer in model clouds of increasing degrees of asymmetry and in a spherical cloud core heated by a neighboring star.

With our method, we confirm the results of ZWG, Evans et al. (2001), and Stamatellos & Whitworth (2003) that the dust temperature in cloud models with physical characteristics typical of dense cores varies monotonically between a minimum value at the center (6–7 K) and a maximum value near the cloud’s edge (14–15 K) for the standard ISRF.

In general, maps at submillimeter wavelengths closely reproduce the column density distribution, even for the most non-isothermal configurations. In the Rayleigh-Jeans regime, in fact, temperature differences of a factor $`2`$ along the line-of-sight are washed out by the much larger variation of the density. However, the decrease of the dust temperature at the cloud center has the consequence that submillimeter maps sample preferentially the low-density envelope of prestellar cores, and are therefore unable to distinguish centrally peaked (pivotal) configurations from more flattened density profiles. We also note that our three-dimensional model based on the work of Galli et al. (2001) produces in natural fashion the cometary structure observed towards many prestellar cores.

For a spherical core heated by an external stellar source close to the cloud’s surface, the temperature distribution is non-spherical, but the submillimeter emission generally reproduces the core’s column density profile. Deviations from circular symmetry become extremely apparent on the other hand shortward of 100 $`\mu `$m where one samples mainly the heated surface layer.

###### Acknowledgements.

We thank the referee S. P. Goodwin for useful comments and A. Natta for constructive criticism. We acknowledge financial support from the EC Research Training Network “The Formation and Evolution of Young Stellar Clusters” (HPRN-CT-2000-00155). JG acknowledges support from the scholarship SFRH/BD/6108/2001 awarded by the Fundação para a Ciência e Tecnologia (Portugal).


format: text

text:
# Simulations of Large-Scale Structure in the New Millennium

## 1. Introduction: evolution of cosmological simulations

Cosmological $`N`$-body simulations started in late 1970’s, and since then have played an important part in describing and understanding the nonlinear gravitational clustering in the universe. As far as I know, the cosmological N-body simulation in a comoving periodic cube, which is quite conventional now, was performed for the first time by Miyoshi & Kihara (1975) using $`N=400`$ particles. Figure 1 plots the evolution of the number of particles employed in cosmological N-body simulations. Here I consider only the “high-resolution” simulations including Particle-Particle, Particle-Particle–Particle-Mesh, and tree algorithms which are published in refereed journals (excluding, e.g., conference proceedings). I found that the evolution is well fitted by

$`N=400\times 10^{0.215(\mathrm{Year}1975)},`$ (1)
where the amplitude is normalized to the work of Miyoshi & Kihara (1975). Just for comparison, the total number of CDM particles of mass $`m_{\mathrm{CDM}}`$ in a box of the universe of one side $`L`$ is

$`N={\displaystyle \frac{\mathrm{\Omega }_{\mathrm{CDM}}\rho _{\mathrm{cr}}L^3}{m_{\mathrm{CDM}}}}10^{83}\left({\displaystyle \frac{\mathrm{\Omega }_{\mathrm{CDM}}}{0.23}}\right)\left({\displaystyle \frac{L}{1h^1\mathrm{Gpc}}}\right)^3\left({\displaystyle \frac{1\mathrm{k}\mathrm{e}\mathrm{V}}{m_{\mathrm{CDM}}}}\right)\left({\displaystyle \frac{0.71}{h}}\right).`$ (2)
If I simply extrapolate equation (1) and adopt the WMAP parameters (Spergel et al. 2003), then the number of particles that one can simulate in a $`(1h^1\mathrm{Gpc})^3`$ box will reach the real number of CDM particles in December 2348 and February 2386 for $`m_{\mathrm{CDM}}=1\mathrm{k}\mathrm{e}\mathrm{V}`$ and $`10^5\mathrm{eV}`$, respectively. I have not yet checked the above arithmetic, but the exact number should not change the basic conclusion; simulations in the new millennium will be unbelievably realistic.

## 2. 1970’s: simulating nonlinear gravitational evolution

It is not easy to identify who attempted seriously for the first time the numerical simulation of large-scale structure in the universe. Still I believe that a paper by Miyoshi & Kihara (1975) is truly pioneering, and let me briefly mention it here. They carried out a series of cosmological $`N`$-body experiments with $`N=400`$ galaxies (=particles) in an expanding universe. The simulation was performed in a comoving cube with a periodic boundary condition (Fig.2). As the title of the paper “Development of the correlation of galaxies in an expanding universe” clearly indicates, they were interested in understanding why the observed galaxies in the universe exhibit a characteristic correlation function of $`g(r)=(r_0/r)^s`$. In fact, Totsuji & Kihara (1969) had already found that $`r_0=4.7h^1`$Mpc and $`s=1.8`$ is a reasonable fit to the clustering of galaxies in the Shane – Wirtanen catalogue. One of the main conclusions of Miyoshi & Kihara (1975) is that “The power-type correlation function $`g(r)=(r_0/r)^s`$ with $`s2`$ is stable in shape; it is generated from motionless galaxies distributed at random and also from a system with weak initial correlation ”.

Several years after Totsuji & Kihara (1969) published the paper, Peebles (1974) and Groth & Peebles (1977) reached the same conclusion independently, which has motivated several cosmological $`N`$-body simulations all over the world. Among others, Aarseth, Gott & Turner (1979) conducted a series of careful and systematic simulations to explore nonlinear gravitational clustering. Those simulations in 1970’s assume that galaxy distribution is well traced by simulations particles. In fact, the above papers spent many pages in an argument to justify the assumption, and then attempted to understand the nonlinear gravitational evolution and to quantitatively describe the large-scale structure in the computer on the basis of two-point correlation functions. In this sense, I would say that the simulations in late 1970’s are more physics-oriented rather than astronomy. Also it is interesting to mention that Prof. Kihara was a solid state physicist in the University of Tokyo and I speculate that this is why he was able to accomplish truly pioneering work from such an interdisciplinary point of view.

## 3. 1980’s: introducing galaxy biasing

The primary goal of simulations of large-scale structure in 1980’s was to predict observable galaxy distribution from dark matter clustering. Necessarily one had to start distinguishing galaxies and simulation particles (designed to represent dark matter in the universe), i.e., to introduce the notion of galaxy biasing according to the current terminology. Furthermore a variety of astronomical and/or observational effects (selection function, redshift-space distortion, etc.) had to be incorporated towards more realistic comparison with galaxy redshift survey data which became available those days. Davis, Efstathiou, Frenk and White (1985) is the most influential and seminal paper in cosmological $`N`$-body simulations in my view. While their work is quite pioneering in many aspects, the most important message that they were able to show in a quite convincing fashion is that simulations of large-scale structure can provide numerous realistic and testable predictions of dark matter scenarios against the observational data from luminous galaxy samples. Considering the fact that they used only $`N=32^3`$ particles and thus had to identify the present epoch as when they advance the simulation merely by a factor of 1.4 relative to the initial epoch, this presents a convincing case that the most important is not the quality of simulations but those who interpret the result.

## 4. 1990’s: more accurate and realistic modeling of galaxy clustering

I started to work on cosmological N-body simulations around 1987, and it has been my major research topic for the next several years. At that time I often asked myself if purely N-body simulations would continue to advance our understanding of galaxy clustering significantly. My personal answer was “No. Without proper inclusion of hydrodynamics, radiative processes, star formation and feedback, it is unlikely to proceed further”, so I moved to more analytical and/or observational researches. Although I still do not think that my thought was terribly wrong, I have to admit that my decision was premature; purely N-body simulations in 1990’s turned out to be so successful and they achieved quite important contributions in (at least) three basic aspects; (i) accurate modeling of nonlinear two-point correlation functions, (ii) abundance and biasing of dark matter halos, and (iii) density profiles of dark halos, which are separately described below. Thus I returned to simulation work again in late 1990’s.

### 4.1. Nonlinear two-point correlation function of dark matter

The first breakthrough came from the discovery of the amazing scaling property in the two-point correlation functions (Hamilton et al. 1991). They found that the two-point correlation functions in N-body simulations can be well approximated by a universal fitting formula which empirically interpolates the linear regime and the nonlinear stable solution. Their remarkable insight was then elaborated and improved later (e.g., Peacock & Dodds 1996; Smith et al. 2003), and the resulting accurate fitting formulae have been applied in a variety of cosmological analyses.

Figure 3 plots two-point correlations of dark matter from N-body simulations (Hamana, Colombi & Suto 2001a). The symbols indicate the averages over the five realizations from simulations in real space (open circles) and in redshift space (solid triangles), and the quoted error-bars represent the standard deviation among them. The results for all particles (left panels) agree very well with the theoretical predictions (solid lines) which combine the Peacock-Dodds formula and the light-cone effect (e.g., Yamamoto & Suto 1999; Suto et al. 1999). The scales where the simulation data in real space become smaller than the corresponding theoretical predictions simply reflect the force resolution of the simulations. The result is fairly robust against the selection effect; Figure 4 indicates that the simulation results and the predictions are still in good agreement even after incorporating the realistic selection functions.

One of the main purposes of N-body simulations in 1970’s and 1980’s was to compute the nonlinear two-point correlation functions which were unlikely to predict analytically with a reasonable accuracy. In the light of this, it is interesting to note that as long as two-point correlation functions of dark matter are concerned, one does not have to run N-body simulations owing to the significant progress in semi-analytical modeling achieved on the basis of previous N-body simulations.

### 4.2. Biasing of dark matter halos

The second remarkable progress where N-body simulations have played a major role in 1990’s is related to the statistics of dark halos, in particular their mass function and spatial biasing. The standard picture of structure formation predicts that the luminous objects form in a gravitational potential of dark matter halos. Therefore, a detailed understanding of halo clustering is the natural next step beyond the description of clustering of dark matter particles. Both the extended Press-Schechter theory and high-resolution N-body simulations have made significant contributions in constructing a semi-analytical framework for halo clustering.

For a specific example, let me show our recent mass-, scale-, and time-dependent halo bias model (Hamana et al. 2001b):

$`b_{\mathrm{halo}}(M,R,z)`$ $`=`$ $`b_{\mathrm{ST}}(M,z)\left[1.0+b_{\mathrm{ST}}(M,z)\sigma _M(R,z)\right]^{0.15},`$ (3)
$`b_{\mathrm{ST}}(M,z)`$ $`=`$ $`1+{\displaystyle \frac{\nu 1}{\delta _c(z)}}+{\displaystyle \frac{0.6}{\delta _c(z)(1+0.9\nu ^{0.3})}},`$ (4)
which generalizes the previous work including Mo & White (1996), Jing (1998), Sheth & Tormen (1999) and Taruya & Suto (2000). The above biasing parameter is adopted for $`R>2R_{\mathrm{vir}}(M,z)`$, where $`R_{\mathrm{vir}}(M,z)`$ is the virial radius of the halo of mass $`M`$ at $`z`$, while we set $`b_{\mathrm{halo}}(M,R,z)=0`$ for $`R<2R_{\mathrm{vir}}(M,z)`$ in order to incorporate the halo exclusion effect approximately. In the above expressions, $`\sigma _M(R,z)`$ is the mass variance smoothed over the top-hat radius $`R(3M/4\pi \rho _0)^{1/3}`$, $`\rho _0`$ is the mean density, $`\delta _c(z)=3(12\pi )^{2/3}/20D(z)`$, $`D(z)`$ is the linear growth rate of mass fluctuations, and $`\nu =[\delta _c(z)/\sigma _M(R,z=0)]^2`$.

Figure 5 shows the comparison of the two-point correlation functions of dark matter halos between the above semi-analytical predictions and the simulation data (Hamana et al. 2001b). For this purpose, we analyze “light-cone output” of the Hubble Volume $`\mathrm{\Lambda }`$CDM simulation (Evrard et al. 2002). The two-point correlation functions on the light-cone plotted in Figure 5 correspond to halos with $`M>5.0\times 10^{13}h^1M_{}`$, $`M>2.2\times 10^{13}h^1M_{}`$ and dark matter from top to bottom. The range of redshift is $`0<z<1`$ (Left panel) and $`0.5<z<2`$ (Right panel). Predictions in redshift and real spaces are plotted in dashed and solid lines, while simulation data in redshift and real spaces are shown in filled triangles and open circles, respectively.

Our model and simulation data show quite good agreement for dark halos at scales larger than $`5h^1`$Mpc. Below that scale, they start to deviate slightly in a complicated fashion depending on the mass of halo and the redshift range. Nevertheless the clustering of clusters on scales below $`5h^1`$Mpc is difficult to determine observationally anyway, and our model predictions differ from the simulation data only by $`20`$ percent at most. This illustrates the fact that the clustering not only of dark matter but also of dark halos, at least as far as their two-point statistics is concerned, can be described well semi-analytically without running expensive N-body simulations at all.

### 4.3. Density profiles of dark halos

The third, and perhaps the most useful in cosmological applications, result out of N-body simulations in 1990’s is the discovery of the universal density profile of dark halos.

The study of the density profiles of cosmological self-gravitating systems or dark halos has a long history. Navarro, Frenk & White (1995, 1996, 1997) found that all simulated density profiles can be well fitted to the following simple model (now generally referred to as the NFW profile):

$$\rho (r)\frac{1}{(r/r_\mathrm{s})(1+r/r_\mathrm{s})^2}$$
(5)
by an appropriate choice of the scaling radius $`r_\mathrm{s}=r_\mathrm{s}(M)`$ as a function of the halo mass $`M`$. Subsequent higher-resolution simulations (Fukushige & Makino 1997, 2001; Moore et al. 1998; Jing & Suto 2000) have indicated that the inner slope of density halos is steeper than the NFW value, and the current consensus among most N-body simulators is given by

$$\rho (r)\frac{1}{(r/r_\mathrm{s})^\alpha (1+r/r_\mathrm{s})^{3\alpha }}$$
(6)
with $`\alpha 1.5`$ rather than the NFW value, $`\alpha =1`$ for $`r>0.01r_{\mathrm{vir}}`$.

Actually it is rather surprising that the fairly accurate scaling relation applies after the spherical average despite the fact that the departure from the spherical symmetry is quite visible in almost all simulated halos. A more realistic modeling of dark matter halos beyond the spherical approximation is important in understanding various observed properties of galaxy clusters and non-linear clustering of dark matter. In particular, the non-sphericity of dark halos is supposed to play a central role in the X-ray morphologies of clusters, in the cosmological parameter determination via the Sunyaev-Zel’dovich effect and in the prediction of the cluster weak lensing and the gravitational arc statistics (Bartelmann et al. 1998; Meneghetti et al. 2000, 2001).

Recently Jing & Suto (2002) presented a detailed non-spherical modeling of dark matter halos on the basis of a combined analysis of the high-resolution halo simulations (12 halos with $`N10^6`$ particles within their virial radius) and the large cosmological simulations (5 realizations with $`N=512^3`$ particles in a $`100h^1`$Mpc boxsize). The density profiles of those simulated halos are well approximated by a sequence of the concentric triaxial distribution with their axis directions being fairly aligned:

$`R^2(\rho _\mathrm{s})={\displaystyle \frac{X^2}{a^2(\rho _\mathrm{s})}}+{\displaystyle \frac{Y^2}{b^2(\rho _\mathrm{s})}}+{\displaystyle \frac{Z^2}{c^2(\rho _\mathrm{s})}}.`$ (7)
The origin of the coordinates is set at the center of mass of each surface, and the principal vectors $`𝐚`$, $`𝐛`$ and $`𝐜`$ ($`abc`$) are computed by diagonalizing the inertial tensor of particles in the isodensity surface $`\rho =\rho _\mathrm{s}`$. Figure 6 plots the density profiles measured in this way for individual halos as a function of $`R`$, which indicates that equation (6) is still a good approximation if the spherical radius $`r`$ is replaced by equation (7). The application of the above triaxial modeling for the X-ray, Sunyaev-Zel’dovich, and lensing data is in progress (Lee & Suto 2003, 2004; Oguri, Lee & Suto 2003).

## 5. 2000’s: galaxies in cosmological hydrodynamic simulations

Although serious attempts to create galaxies phenomenologically but directly from cosmological hydrodynamic simulations have been initiated in early 1990’s (e.g., Cen & Ostriker 1992; Katz, Hernquist & Weinberg 1992), those resulting simulated galaxies are far from realistic and there are still plenty of room for improvement. Thus this is one of the most important, and yet quite realistic, goals for the simulations in the new millennium, or hopefully in this decade.

Let me show the result of Yoshikawa et al. (2001) for an example of such approaches. They apply cosmological smoothed particle hydrodynamic simulations in a spatially-flat $`\mathrm{\Lambda }`$-dominated CDM model with particular attention to the comparison of the biasing of dark halos and simulated galaxies. Figure 7 illustrates the distribution of dark matter particles, gas particles, dark halos and galaxies at $`z=0`$. Clearly galaxies are more strongly clustered than dark halos. In order to quantify the effect, we define the following biasing parameter:

$$b_{\xi ,i}(r)\sqrt{\frac{\xi _{ii}(r)}{\xi _{\mathrm{mm}}(r)}},$$
(8)
where $`\xi _{ii}(r)`$ and $`\xi _{\mathrm{mm}}(r)`$ are two-point correlation functions of objects $`i`$ and of dark matter, respectively. Furthermore for each galaxy identified at $`z=0`$, we define its formation redshift $`z_\mathrm{f}`$ by the epoch when half of its cooled gas particles satisfy our criteria of galaxy formation. Roughly speaking, $`z_\mathrm{f}`$ corresponds to the median formation redshift of stars in the present-day galaxies. We divide all simulated galaxies at $`z=0`$ into two populations (the young population with $`z_\mathrm{f}<1.7`$ and the old population with $`z_\mathrm{f}>1.7`$) so as to approximate the observed number ratio of $`3/1`$ for late-type and early-type galaxies.

The difference of the clustering amplitude can be also quantified by their two-point correlation functions at $`z=0`$ as plotted in Figure 8. The old population indeed clusters more strongly than the mass, while the young population is anti-biased. The relative bias between the two populations $`b_{\xi ,\mathrm{g}}^{\mathrm{rel}}\sqrt{\xi _{\mathrm{old}}/\xi _{\mathrm{young}}}`$ ranges $`1.5`$ and 2 for $`1h^1\text{Mpc}<r<20h^1\text{Mpc}`$, where $`\xi _{\mathrm{young}}`$ and $`\xi _{\mathrm{old}}`$ are the two-point correlation functions of the young and old populations. It is interesting to note that even this crude approach is able to explain the morphological-dependence of bias, although still in a rather quantitative manner, derived later by Kayo et al. (2004) for SDSS galaxies. With the still on-going rapid progress of observational exploration (e.g., Lahav & Suto 2003 for a recent review on galaxy redshift survey), understanding galaxy biasing as a function of galaxy properties is definitely one of the unsolved important questions in observational cosmology, and the present result indicates that the formation epoch of galaxies plays a crucial role in the morphological segregation.

## 6. Distribution of dark baryons

Finally let me briefly mention yet another possibility of tracing the large-scale structure of the universe using the oxygen emission lines. It is widely accepted that our universe is dominated by dark components; 23 percent of dark matter, and 73 percent of dark energy (e.g., Spergel et al. 2003). Furthermore, as Fukugita, Hogan & Peebles (1998) pointed out earlier, even most of the remaining 4 percent of the cosmic baryons has evaded the direct detection so far, i.e., most of the baryons is indeed dark. Subsequent numerical simulations (e.g., Cen & Ostriker 1999a, 1999b; Davé et al. 2001) indeed suggest that approximately 30 to 50 percent of total baryons at $`z=0`$ take the form of the warm-hot intergalactic medium (WHIM) with $`10^5\mathrm{K}<T<10^7\mathrm{K}`$ which does not exhibit strong observational signature.

Figure 9 compares the distribution of WHIM ($`10^5\mathrm{K}<T<10^7\mathrm{K}`$), hot intracluster gas ($`T>10^7\mathrm{K}`$), and dark matter particles from cosmological smoothed-particle hydrodynamic simulations (Yoshikawa et al. 2001). Clearly WHIM traces the large-scale filamentary structure of mass distribution more faithfully than the hot gas which preferentially resides in clusters that form around the knot-like intersections of those filamentary regions. In order to carry out a direct and homogeneous survey of elusive dark baryons, we propose a dedicated soft-X-ray mission, DIOS (Diffuse Intergalactic Oxygen Surveyor; see Fig. 10). The detectability of WHIM through Oviii and Ovii emission lines via DIOS was examined in detail by Yoshikawa et al. (2003) assuming a detector which has a large throughput $`S_{\mathrm{eff}}\mathrm{\Omega }=10^2`$ cm<sup>2</sup> deg<sup>2</sup> and a high energy resolution $`\mathrm{\Delta }E=2`$ eV. Their results are summarized in Figure 11; they first create a light-cone output from the hydrodynamic simulation up to $`z=0.3`$, compute the bolometric X-ray surface intensity map, select several target fields and finally compute the corresponding spectra relevant for the DIOS survey. The high-spectral resolution of DIOS enables to identify the redshifts of several WHIMs at different emission energies, i.e., Oxygen emission line tomography of the WHIMs at different locations.

They concluded that within the exposure time of $`T_{\mathrm{exp}}=10^{56}`$ sec DIOS will be able to reliably identify Oviii emission lines (653eV) of WHIM with $`T=10^{67}`$ K and the overdensity of $`\delta =10^{0.52}`$, and Ovii emission lines (561, 568, 574, 665eV) of WHIM with $`T=10^{6.57}`$ K and $`\delta =10^{12}`$. The WHIM in these temperature and density ranges cannot be detected with the current X-ray observations except for the oxygen absorption features toward bright QSOs. DIOS is especially sensitive to the WHIM with gas temperature $`T=10^{67}`$K and overdensity $`\delta =10100`$ up to a redshift of $`0.3`$ without being significantly contaminated by the cosmic X-ray background and the Galactic emissions. Fang et al. (2003) also conducted a similar study and reached quite consistent conclusions. Thus such a mission, hopefully launched in several years, promises to provide a unique and important tool to trace the large-scale structure of the universe via dark baryons.

## 7. Conclusions

It turned out that I was able to review the cosmological simulations only for a time-scale of decades. Still the material may be heavily biased, which I have to apologize for the organizers and possible readers of these proceedings. A millennium is definitely too long for any scientist to make any reliable prediction for its eventual outcome. Thus the number of simulation particles that I predicted in Introduction might sound ridiculous, but in reality the progress in the new millennium may be even more drastic that whatever one can imagine. For instance, it is unlikely that one still continues to use currently popular particle- or mesh-based simulation techniques over next hundreds of years. In that case, the number of particles may turn out to be a totally useless measure of the progress or reliability of simulations. Nevertheless I believe that a historical lesson that I learned in preparing this talk will be still true even at the end of this millennium; good science favors the prepared mind, not the largest simulation at the time.

#### Acknowledgments.

I thank Stephan Colombi, Gus Evrard, Takashi Hamana, Y.P.Jing, Atsushi Taruya, Naoki Yoshida, and Kohji Yoshikawa for enjoyable collaborations. Naoki Yoshida also encouraged me to plot Figure 1 in order to illustrate the progress of cosmological N-body simulations. I am also grateful to Ed Turner for providing me a digitized version of the first movie of his cosmological N-body simulations that I was able to present in the talk. My presentation file for the symposium may be found in the PDF format at http://www-utap.phys.s.u-tokyo.ac.jp/~ suto/mypresentation\_2003e. This research was supported in part by the Grant-in-Aid for Scientific Research of JSPS. Numerical computations presented in this paper were carried out mainly at ADAC (the Astronomical Data Analysis Center) of the National Astronomical Observatory, Japan (project ID: yys08a, mky05a). DIOS (Diffuse Intergalactic Oxygen Surveyor) is a proposal by a group of scientists at Tokyo Metropolitan University, Institute of Astronautical Sciences, the University of Tokyo, and Nagoya University (P.I., Takaya Ohashi).

## References

Aarseth, S.J., Gott, J.R. & Turner, E.L. 1979, ApJ, 228, 664.

Bartelmann, M., Huss, A., Colberg, J. M., Jenkins, A., & Pearce, F. R. 1998, A&A, 330, 1.

Bode, P., Bahcall, N. A., Ford, E. B., Ostriker, J.P. 2001, ApJ, 551, 15

Bode, P., Ostriker, J. P., & Xu, G. 2000, ApJS, 128, 561.

Cen, R., & Ostriker, J.P. 1992, ApJ, 399, L113.

Cen, R., & Ostriker, J. 1999a, ApJ, 514, 1.

Cen, R., & Ostriker, J. 1999b, ApJ, 519, L109.

Davé, R., et al. 2001, ApJ, 552, 473.

Davis, M., Efstathiou, G., Frenk, C.S., & White, S.D.M. 1985, ApJ, 292, 371.

Efstathiou, G., & Eastwood, J.W. 1981, MNRAS, 194, 503.

Evrard, A.E., et al. 2002, ApJ, 573, 7.

Fang, T. et al. 2003, ApJ, submitted (astro-ph/0311141).

Fukugita, M., Hogan, C.J., & Peebles, P.J.E., 1998, ApJ, 503, 518.

Fukushige, T., & Makino, J. 1997, ApJ, 477, L9.

Fukushige, T., & Makino, J. 2001, ApJ, 557, 533.

Gelb, J. M., & Bertschinger, E. 1994, ApJ, 436, 467.

Groth,E.J., & Peebles, P.J.E. 1977, ApJ, 217, 385.

Hamana, T., Colombi, S., & Suto, Y. 2001a, A& A, 367, 18.

Hamana, T., Yoshida, N., Suto, Y., & Evrard, A.E. 2001b, ApJ, 561, L143.

Hamilton, A.J.S., Matthews, A., Kumar, P., & Lu, E. 1991, ApJL, 374, L1.

Jing, Y.P. 1998, ApJ, 503, L9.

Jing, Y.P., & Suto, Y. 1998, ApJL, 494, L5.

Jing, Y. P., & Suto, Y. 2000, ApJ, 529, L69.

Jing, Y. P., & Suto, Y. 2002, ApJ, 574, 538.

Kayo, I. et al. 2004, in preparation.

Lahav, O., & Suto, Y. 2003, Living Reviews in Relativity, in press (astro-ph/0310642).

Lee, J., & Suto, Y. 2003, ApJ, 585, 151.

Lee, J., & Suto, Y. 2004, ApJ, 601, February 1 issue, in press (astro-ph/0306217).

Katz, N., Hernquist, L., Weinberg, D. H. 1992, ApJ, 399, L109.

Meneghetti, M., Bolzonella, M., Bartelmann, M., Moscardini, L., & Tormen, G. 2000, MNRAS, 314, 338.

Meneghetti, M., Yoshida, N., Bartelmann, M., Moscardini, L., Springel, V., Tormen, G., & White S. D. M. 2001, MNRAS, 325, 435.

Miyoshi, K. & Kihara, T. 1975, Publ.Astron.Soc.Japan., 27, 333.

Mo, H.J., & White, S.D.M 1996,MNRAS, 282, 347.

Moore, B., Governato, F., Quinn, T., Stadel, J., & Lake, G. 1998, ApJ, 499, L5.

Navarro, J.F., Frenk, C.S., & White, S.D.M. 1995, MNRAS 275, 720.

Navarro, J.F., Frenk, C.S., & White, S.D.M. 1996, ApJ, 462, 563.

Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493.

Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, in press (astro-ph/0306102) .

Peacock, J.A., & Dodds, S.J. 1996, MNRAS, 280, L19.

Peebles,P.J.E. 1974, ApJ, 189, L51.

Sheth, R.K., & Tormen, G. 1999, MNRAS, 308, 119.

Smith, R. E. et al. 2003, MNRAS, 341, 1311.

Spergel, D.N. et al. 2003, ApJS, 148, 175.

Suto, Y., Magira, H., Jing, Y. P., Matsubara, T., & Yamamoto, K. 1999, Prog.Theor.Phys.Suppl., 133, 183.

Taruya, A. & Suto,Y. 2000, ApJ, 542, 559.

Totsuji, H. & Kihara, T. 1969, Publ.Astron.Soc.Japan., 21, 221.

Yamamoto, K., & Suto, Y. 1999, ApJ, 517, 1.

Yoshikawa, K., Taruya, A., Jing, Y.P., & Suto, Y. ApJ, 2001, 558, 520.

Yoshikawa, K., Yamasaki, N.Y., Suto, Y., Ohashi, T., Mitsuda, K., Tawara, Y. & Furuzawa, A. 2003, PASJ, 55, 879.


text:
# Field-induced magnetic order in La2-xSrxCuO4 (𝑥= 0.10, 0.115, 0.13) studied by in-plane thermal conductivity measurements

## I Introduction

The so-called 1/8 anomaly, namely, the anomalous suppression of superconductivity at $`p`$ (the hole concentration per Cu in the CuO2 plane) $``$ 1/8 in the La-based high-$`T_\mathrm{c}`$ superconductors has been understood in terms of the stripe model proposed by Tranquada et al.rf:Tranquada . From their neutron scattering experiment, a static stripe order of spins and holes was suggested to be formed at low temperatures in La1.48Nd0.4Sr0.12CuO4. Meanwhile, a spatially modulated dynamical spin correlation has been found to exist in a wide range of $`p`$ from the underdoped to overdoped region in La2-xSrxCuO4rf:Cheong ; rf:Mason ; rf:Thurston ; rf:Yamada . Since the dynamical spin correlation may be regarded as a spin part of the dynamical stripe correlations of spins and holes, it has been understood that the dynamical stripes tend to be statically stabilized around $`p=`$ 1/8 and are easily pinned, leading to the static stripe order and the 1/8 anomalyrf:Tranquada ; rf:Kumagai ; rf:Luke ; rf:Watanabe ; rf:Watanabe2 ; rf:Koike0 ; rf:Koike ; rf:Adachi . For the pinning, the buckling of the CuO2 plane in the tetragonal low-temperature (TLT) structure (space group: P42/ncm) is regarded as being effective in La1.6-xNd0.4SrxCuO4 or La2-xBaxCuO4. Moreover, a small amount of nonmagnetic impurities such as Zn has also been found to be effective for the pinning from the transport measurementsrf:Koike0 ; rf:Koike ; rf:Adachi  and also from the muon-spin-relaxation ($`\mu `$SR) measurementsrf:muSR0 ; rf:muSR1 ; rf:muSR2 ; rf:muSR3 ; rf:muSR3-2 ; rf:muSR4 .

Recently, Katano et al.rf:Katano  and Lake et al.rf:Lake010  have found from the neutron scattering experiments in magnetic fields that the intensity of the incommensurate magnetic Bragg peak corresponding to the long-range static stripe order is enhanced around $`p=`$ 1/8 in La2-xSrxCuO4 by the application of magnetic fields parallel to the $`c`$-axis. The enhancement of the magnetic Bragg peak intensity is marked in $`x=`$ 0.10rf:Lake010 , while it is observable but small in $`x=`$ 0.12rf:Katano . A similar enhancement has also been found in the excess-oxygen-doped La2CuO4+δrf:Khaykovich . Such a field-induced magnetic order may be interpreted as being due to the pinning of the dynamical stripes by induced vortex cores in the $`ab`$-plane, though the pinning effect of vortex cores is not certain. From the NMR measurements of the nearly optimally doped Tl2Ba2CuO6+δ, on the other hand, a magnetically ordered state has been detected only inside vortex cores, while the electronic state outside vortex cores remains superconductingrf:KumagaiNMR . This seems to be inconsistent with the above results of the neutron scattering experiments suggesting the development of the long-range magnetic order in magnetic fields. As for the theoretical study on the field-induced magnetic order, before the report of these experimental results, the SO(5) theory had already predicted that the electronic state in vortex cores was an antiferromagnetically ordered onerf:Zhang ; rf:Arovas . According to the recent theoretical study by Demler et al.rf:Demler , on the other hand, the results of the neutron scattering experiments in a magnetic field may be explained as being due to the approach of the superconducting phase to the coexisting phase of the superconducting state and the magnetically ordered state as a result of the increase of the spin fluctuations with low energy induced by the proximity effect. Accordingly, the origin of the field-induced magnetic order in the La-based high-$`T_\mathrm{c}`$ cuprates has not yet been settled.

The thermal conductivity measurement is a renewed technique for study of the spin and charge state in transition-metal oxides. For example, it has been found that the thermal conductivity due to phonons is markedly enhanced in the spin-gap state of CuGeO3rf:AndoCuGe  and SrCu2(BO3)2rf:KudoSr ; rf:Hof ; rf:KudoSr2  and also in the charge-ordered state of Sr1.5La0.5MnO4rf:Hess . The enhancement of the thermal conductivity due to phonons has also been reported in the static stripe-ordered state of La1.28Nd0.6Sr0.12CuO4rf:Baberski  and La2-xSrxNiO4rf:Hess , while the suppression of the thermal conductivity due to phonons has been reported in the possible short-range (dynamical) stripe-ordered state of YBa2Cu3O7-δ and HgBa2Can-1CunO2n+2+δ (n $`=`$ 1, 2, 3) around $`p=`$ 1/8rf:Cohn . In the antiferromagnetically ordered state of La2CuO4 rf:214-1 ; rf:214-2 ; rf:214-3  and YBa2Cu3O6rf:123 , large contribution of magnons to the thermal conductivity in the $`ab`$-plane has been observed.

In this paper, in order to investigate the field-induced magnetically ordered state of La2-xSrxCuO4 around $`p=`$ 1/8, we have performed the thermal conductivity measurements in the $`ab`$-plane of La2-xSrxCuO4 ($`x=`$ 0.10, 0.115, 0.13) and 1 % Zn-substituted La2-xSrxCu0.99Zn0.01O4 ($`x=`$ 0.115) single crystals in magnetic fields parallel to the $`c`$-axis and also parallel to the $`ab`$-plane. This is because the suppression of the superconductivity occurs around $`x=`$ 0.115 in La2-xSrxCuO4 and because the suppression becomes most marked at $`x=`$ 0.115 through the partial substitution of Zn for Curf:Koike115 .

## II experimental

Single crystals were grown by the Traveling-Solvent Floating-Zone method in a similar way to that described in Ref. rf:Kawamata . Thermal conductivity measurements were carried out by a conventional steady-state method. One side of a rectangular single-crystal, whose typical dimensions were $`3\times 1\times 1`$ mm3, was anchored on the copper heat sink with indium solder. A chip-resistance of 1 k$`\mathrm{\Omega }`$ (Alpha Electronics Corp., MP1K000) was attached as a heater to the opposite side of the single crystal with GE7031 varnish. The temperature difference across the crystal (0.02–1.0 K) was measured with two Cernox thermometers (LakeShore Cryotronics, Inc., CX-1050-SD). Magnetic fields up to 14 T were applied parallel to the $`c`$-axis and also parallel to the $`ab`$-plane, using a superconducting magnet.

## III Results and Discussion

Figures 1 (a)-(c) show the temperature dependence of the thermal conductivity in the $`ab`$-plane, $`\kappa _{\mathrm{ab}}`$, of La2-xSrxCuO4 ($`x=`$ 0.10, 0.115, 0.13) in magnetic fields parallel to the $`c`$-axis. As formerly reported by Nakamura et al.rf:214-1 , $`\kappa _{\mathrm{ab}}`$ in zero field decreases with decreasing temperature, and it exhibits a slight enhancement below the superconducting transition temperature $`T_\mathrm{c}`$, defined as the onset temperature of the Meissner diamagnetism. By the application of magnetic fields, $`\kappa _{\mathrm{ab}}`$ is suppressed at low temperatures. Strictly looking, it appears that the temperature $`T_\kappa `$, below which $`\kappa _{\mathrm{ab}}`$ is suppressed, is located above $`T_\mathrm{c}`$ and is almost independent of the magnitude of the magnetic field. The suppression by the application of magnetic fields is marked in $`x=`$ 0.10 and 0.13, while it is small in $`x=`$ 0.115. By the application of magnetic fields parallel to the $`ab`$-plane, on the other hand, the suppression is so small in all $`x`$ that $`T_\kappa `$ is hard to be determined, as shown in Figs. 1 (d)-(f).

Figures 2 (a)-(c) display the $`x`$ dependence of the relative reduction of $`\kappa _{\mathrm{ab}}`$, $`[\kappa _{\mathrm{ab}}(0\mathrm{T})\kappa _{\mathrm{ab}}(14\mathrm{T})]/\kappa _{\mathrm{ab}}(0\mathrm{T})`$, in a magnetic field of 14 T parallel to the $`c`$-axis and also parallel to the $`ab`$-plane. It is clearly seen that the relative reduction in a field parallel to the $`c`$-axis exhibits a dip at $`x=`$ 0.115 at low temperatures below $`T_\kappa `$, while that in the field parallel to the $`ab`$-plane is almost independent of $`x`$ in the whole temperature range. The former $`x`$ dependence with a dip-like shape is similar to the dip-like shape around $`p=`$ 1/8 in the $`T_\mathrm{c}`$ vs $`x`$ diagram, being associated with the 1/8 anomaly and the stripe order.

First, we try explaining the observed field-dependence of $`\kappa _{\mathrm{ab}}`$ in terms of conventional mechanisms. It is well known that the thermal conductivity in the high-$`T_\mathrm{c}`$ cuprates exhibits an increase just below $`T_\mathrm{c}`$ with decreasing temperature, which is explained as being due to the increase of the thermal conductivity due to both phonons and electrons on account of the decrease of the phonon-electron scattering rate and the increase of the life time of quasiparticles, respectivelyrf:Hagen-ph ; rf:Yu-ele . By the application of magnetic fields, the thermal conductivity due to quasiparticles is expected to be enhanced in such d-wave superconductors as the high-$`T_\mathrm{c}`$ cuprates owing to the Volovik effectrf:Volovik , namely, the increase of the density of states at the Fermi level, because the supercurrent surrounding vortices causes the Doppler shift of the energy of quasiparticles around the nodes of the superconducting gaprf:Sun1 ; rf:Izawa . In detail, recent theories have suggested two kinds of pictures for the state of quasiparticles in magnetic fields. One is a group of discrete states induced by the Andreev reflection which are composed of localized states within vortex coresrf:Kopnin1 ; rf:Kopnin2  and quantized states expanding outside vortex coresrf:Anderson ; rf:Gorkov ; rf:Janko ; rf:Kopnin3 . The other is an energy-band state due to the periodic vortex potentialrf:Melnikov ; rf:Franz . In any case, however, these are contrary to the present experimental result. In order to understand the suppression of the thermal conductivity due to quasiparticles in magnetic fields, we should consider that the increase of quasiparticles leads to the increase of the electron-electron scattering rate and/or the increase of the electron-vortex scattering rate. The latter has been proposed to understand the field dependence of the in-plane thermal conductivity in Bi2Sr2CaCu2O8rf:Krishana ; rf:Ando2212 . As for the thermal conductivity due to phonons, on the other hand, it is expected to be suppressed by the application of magnetic fields owing to the increase of the phonon-vortex scattering rate and/or the increase of the phonon-electron scattering rate, because the number of quasiparticles increases due to the Doppler shift. Hence, at a glance, the suppression of the thermal conductivity does not seem anomalous at all. However, these conventional mechanisms can not explain the fairly small field-dependence of $`\kappa _{\mathrm{ab}}`$ in $`x=`$ 0.115 compared with $`\kappa _{\mathrm{ab}}`$ in $`x=`$ 0.10 and 0.13. Therefore, a new concept is necessary to explain the suppression of the thermal conductivity by the application of magnetic fields in La2-xSrxCuO4.

Next, in order to explain the anomalous $`x`$ dependence of the suppression of the thermal conductivity in magnetic fields, we focus our mind on the development of the static stripe order in magnetic fields. In the following discussion, the field-induced magnetic order is regarded as the field-induced static stripe order. For comparison with the neutron scattering data, the temperature dependence of the relative reduction of the thermal conductivity, $`[\kappa _{\mathrm{ab}}(0\mathrm{T})\kappa _{\mathrm{ab}}(H)]/\kappa _{\mathrm{ab}}(0\mathrm{T})`$, in several magnetic fields $`H`$ parallel to the $`c`$-axis is plotted in Figs. 3(a)-(c). For $`x=`$ 0.10, it is found that the relative reduction rapidly increases at low temperatures below $`T_\kappa `$ with decreasing temperature and that $`T_\kappa `$ is almost independent of the magnitude of the magnetic field. The relative reduction is found to increase with increasing magnetic-field. These temperature and field dependences are quite similar to those of the intensity of the magnetic Bragg peak corresponding to the static stripe order in magnetic fields in $`x=`$ 0.10rf:Lake010 . Moreover, $`T_\kappa `$ of $`x=`$ 0.10 coincides with the temperature $`T_\mathrm{M}`$ of $`x=`$ 0.10 below which the magnetic Bragg peak in magnetic fields developsrf:Lake010 . Therefore, it is naturally understood that the suppression of the thermal conductivity in magnetic fields is related to the development of the static stripe order.

As mentioned in Sec. I, there are some previous reports on the relation between the thermal conductivity due to phonons and the static stripe order. It has been reported that the thermal conductivity due to phonons is enhanced by the formation of the static stripe order for La1.28Nd0.6Sr0.12CuO4rf:Baberski  and La2-xSrxNiO4rf:Hess  and that, on the contrary, it is suppressed in the possible short-range (dynamical) stripe-ordered state of YBa2Cu3O7-δ and HgBa2Can-1CunO2n+2+δ ($`n=`$ 1, 2, 3) around $`p=1/8`$rf:Cohn . The enhancement of the thermal conductivity due to phonons in La1.28Nd0.6Sr0.12CuO4 has been interpreted as being due to the disappearance of the scattering of phonons by the dynamical stripes which make the mean free path of phonons strongly limitedrf:Baberski . In this case, phonons are expected not to be scattered so strongly by the lattice distortion induced by the static stripe order in La1.28Nd0.6Sr0.12CuO4, because the correlation length of the static stripe order is more than $`170\mathrm{\AA }`$ so that the lattice distortion is regarded as being rather periodicrf:Tranquada2 . In the case of YBa2Cu3O7-δ, on the other hand, the suppression of the thermal conductivity due to phonons has been interpreted as being due to the possible formation of the short-range stripe order, namely, the dynamical stripe correlations, which was suggested from the inelastic neutron scattering experiment in YBa2Cu3O6.6 by Mook et al.rf:Mook . Accordingly, it is found to depend on the correlation length of the stripe order whether the thermal conductivity is enhanced or suppressed, but the present suppression of the thermal conductivity in magnetic fields in La2-xSrxCuO4 cannot be explained on these lines, because the correlation length of the static stripe order of $`x=`$ 0.10 in a magnetic field of 14.5 Trf:Lake010  is as long as that of La1.48Nd0.4Sr0.12CuO4 in zero fieldrf:Tranquada2 . According to the detailed investigation on the effect of the structural phase transition on the thermal conductivity for La2-x-yNdySrxCuO4 by Sera et al.rf:Sera , on the other hand, the enhancement of the thermal conductivity is explained as being due to the increase of the phonon velocity through the structural phase transition to the TLT structure. Very recently, moreover, it has been pointed out by Hess et al.rf:Hess3  that the enhancement in the TLT phase of La2-x-yREySrxCuO4 (RE: rare earth element) is probably not due to the formation of the static stripe order but due to the suppression of the lattice instability. At present, therefore, it appears that the formation of the static stripe order does not have a large effect on the thermal conductivity due to phonons and it is hard to clearly explain the present suppression of the thermal conductivity in magnetic fields on the basis of the thermal conductivity due to phonons.

Alternatively, we try to understand the suppression of the thermal conductivity on the basis of the change of the thermal conductivity due to quasiparticles through the field-induced static stripe order. In the CuO2 plane where the static stripe order is formed, charge carriers (holes) are confined in the one-dimensional path of the striperf:AndoR-T , so that they cannot move so easily as in the carrier-homogeneous CuO2 plane, leading to the decrease of the mobility of quasiparticles carrying heat. In fact, it has been pointed out from the electrical resistivityrf:Boe  and thermal conductivityrf:Sun1  measurements that quasiparticles tend to be localized by the application of magnetic fields in the underdoped region of La2-xSrxCuO4. Moreover, when both the static stripe-ordered phase and the superconducting phase coexist and form domains, the mean free path of quasiparticles carrying heat is expected to be reduced by the domain wall. Accordingly, the suppression of the thermal conductivity in magnetic fields in La2-xSrxCuO4 is able to be ascribed to the decrease of the thermal conductivity due to quasiparticles through the development of the static stripe order, and its anomalous $`x`$ dependence is also explained as follows. That is to say, the small field dependence in $`x=`$ 0.115 is attributed to the fact that the static stripe order is already developed in $`x=`$ 0.115 even in zero fieldrf:Suzuki ; rf:Kimura2 , while the large field dependence in $`x=`$ 0.10 and 0.13 is able to be regarded as being due to the marked development of the static stripe order by the application of magnetic fields. As shown in Figs. 3(a)-(c), in fact, the field dependence of the relative reduction of the thermal conductivity is very small in $`x=`$ 0.115, compared with the large field-dependence in $`x=`$ 0.10 and 0.13. Here, it is noted that the temperature dependence of the relative reduction of the thermal conductivity in $`x=`$ 0.115 and 0.13 roughly indicate the temperature dependence of the intensity of the magnetic Bragg peak corresponding to the static stripe order in magnetic fields, as in the case of $`x=`$ 0.10. Furthermore, it is worth while noting that Figs. 3(a)-(c) perhaps represent the development of the charge stripe order in $`x=`$ 0.10, 0.115 and 0.13. By the way, this model gives a possible answer for the well-known question why the enhancement of the thermal conductivity at low temperatures below $`T_\mathrm{c}`$ is relatively small in La2-xSrxCuO4rf:214-1  compared with the other high-$`T_\mathrm{c}`$ cupratesrf:Yu-ele ; rf:Hagen-ph . In the underdoped region, this is because at least a small amount of the static stripe order exists even in zero field to strongly suppress the thermal conductivity, though the magnetic Bragg peak detected in the neutron scattering experiments is very weakrf:Suzuki ; rf:Kimura2 ; rf:Matsushita .

The above model is confirmed by the following experimental result for the 1 % Zn-substituted La2-xSrxCu0.99Zn0.01O4 with $`x=`$ 0.115, as shown in Fig. 4. It is found that the thermal conductivity exhibits no field-dependence in magnetic fields parallel to the $`c`$-axis and also parallel to the $`ab`$-plane. According to the above model, no field-dependence is reasonably understood, because the static stripe order is fully developed even in zero field on account of the strong pinning effect of Znrf:Koike0 ; rf:Koike ; rf:Adachi ; rf:muSR0 ; rf:muSR1 ; rf:muSR2 ; rf:muSR3 ; rf:muSR3-2 ; rf:muSR4  so that further development of the static stripe order by the application of magnetic fields is not expected.

Now, there remains a question why $`T_\kappa `$ is almost independent of the magnitude of the magnetic field, though this question has already been pointed out from the neutron scattering experimentrf:Lake010 . Before answering the question, first we suppose the stripe pinning model by vortex cores. Taking into account the result that the dynamical stripes of spins and holes are pinned in such an inhomogeneous electronic background as the partially Zn-substituted CuO2 planerf:Koike0 ; rf:Koike ; rf:Adachi ; rf:muSR0 ; rf:muSR1 ; rf:muSR2 ; rf:muSR3 ; rf:muSR3-2 ; rf:muSR4 ; rf:Smith ; rf:Tohyama , vortex cores induced in the CuO2 plane are also expected to operate to pin the dynamical stripes, as mentioned in Sec. I. Surely, the existence of impurities such as Zn atoms or vortex cores induces the energy loss in the dynamical stripes, but the most important point is that the energy loss depends on whether an impurity such as Zn or a vortex core is located at a charge stripe or a spin stripe. The stripe order should be located in a way which make the energy loss the smallest, so that the dynamical stripes are pinned, resulting in the development of the static stripe order. The stripe pinning model by vortex cores is supported by the present experimental result that the suppression of the thermal conductivity is not observed by the application of magnetic fields parallel to the $`ab`$-plane. This is because, in magnetic fields parallel to the $`ab`$-plane, vortex cores penetrate the so-called blocking-layer preferably, so that no vortex core appears in the CuO2 plane, leading to neither pinning of the dynamical stripes nor development of the static stripe order. Recurring to the question, the distance between vortex cores is $``$ 130 $`\mathrm{\AA }`$ in a magnetic field of 14 T. Therefore, the pinning effect of vortex cores seems to be so local that $`T_\kappa `$ does not depend on the number of vortex cores at least in magnetic fields up to 14 T. Meanwhile, it is found that $`T_\kappa `$ is located above $`T_\mathrm{c}`$, but this is not inconsistent with the above discussion, because the superconducting fluctuation exists above $`T_\mathrm{c}`$. In fact, recent reports on the Nernst effect have suggested that the vortex state survives even above $`T_\mathrm{c}`$rf:Xu ; rf:Ong  . Thus, the field independence of $`T_\kappa `$ is explained in terms of the stripe pinning model by vortex coresrf:footnote .

Finally, in order to examine the stripe pinning model by vortex cores, we compare the present results with the neutron scattering resultsrf:Lake010 ; rf:Suzuki ; rf:Kimura2 ; rf:Matsushita  and also the recent $`\mu `$SR results revealing a peculiar $`x`$ dependence of the temperature $`T_\mu `$ below which the development of the dynamical spin correlations was detected in zero field from the $`\mu `$SR measurementsrf:highTG0 ; rf:highTG . Figure 5 shows $`x`$ dependences of $`T_\kappa `$, $`T_\mathrm{c}`$, $`T_\mu `$, $`T_\mathrm{M}`$(0) in zero field and $`T_\mathrm{M}`$(14.5 T) in a magnetic field of 14.5 T parallel to the $`c`$-axis. It is found that both $`T_\mathrm{M}`$(0) and $`T_\mu `$ exhibit the maximum at $`x`$ 0.115, while $`T_\kappa `$ exhibits the minimum at $`x=`$ 0.115. That is, the $`x`$ dependences of $`T_\mathrm{M}`$(0) and $`T_\mu `$ are contrary to that of $`T_\kappa `$, so that $`T_\kappa `$ is directly correlated with $`T_\mathrm{c}`$ rather than $`T_\mathrm{M}`$(0) and $`T_\mu `$. This suggests that $`T_\kappa `$ is regarded as the temperature below which vortex cores are formed in the superconducting fluctuation region in the CuO2 plane so as to pin the dynamical stripes and develop the static stripe order. In fact, the observed field dependence of the thermal conductivity is well explained as follows. For $`x=`$ 0.115, $`T_\kappa `$ is lower than $`T_\mathrm{M}`$(0), so that further development of the static stripe order is not large at low temperatures below $`T_\kappa `$, leading to the fairly small field-dependence of the thermal conductivity. For $`x=`$ 0.10 and 0.13, on the other hand, $`T_\kappa `$ is much higher than $`T_\mathrm{M}`$(0), so that the development of the static stripe order is marked at low temperatures below $`T_\kappa `$, leading to the large field-dependence of the thermal conductivity. The result that $`T_\kappa `$ is much higher than $`T_\mathrm{M}`$(0) for $`x=`$ 0.10 and 0.13 is analogous to the result obtained from the $`\mu `$SR measurements that the long-range magnetic ordering temperature $`T_\mathrm{N}`$ is much higher in the lightly Zn-substituted La2-xSrxCu1-yZnyO4 than in the Zn-free La2-xSrxCuO4 for $`x=`$ 0.10 and 0.13rf:muSR1 ; rf:muSR2 ; rf:muSR3 ; rf:muSR3-2 ; rf:muSR4 , suggesting that vortex cores operate to pin the dynamical stripes as Zn atoms. This is consistent with the fact that $`T_\mathrm{M}`$(14.5 T) deviates from $`T_\mathrm{M}`$(0) and roughly coincides with $`T_\kappa `$. Accordingly, the present results strongly suggest that the origin of the field-induced magnetic order is the pinning of the dynamical stripes by vortex cores.

## IV Conclusions

We have measured the thermal conductivity in the $`ab`$-plane of La2-xSrxCuO4 ($`x=`$ 0.10, 0.115, 0.13) in magnetic fields up to 14 T parallel to the $`c`$-axis and also parallel to the $`ab`$-plane. By the application of magnetic fields parallel to the $`c`$-axis, the thermal conductivity has been found to be suppressed at low temperatures below the temperature $`T_\kappa `$ which is located above $`T_\mathrm{c}`$ and almost independent of the magnitude of the magnetic field. The suppression is marked in $`x=`$ 0.10 and 0.13, while it is small in $`x=`$ 0.115. Furthermore, no suppression is observed in the 1 % Zn-substituted La2-xSrxCu0.99Zn0.01O4 with $`x=`$ 0.115. Taking into account the experimental results that the temperature dependence of the relative reduction of the thermal conductivity in $`x=`$ 0.10 is quite similar to the temperature dependence of the intensity of the magnetic Bragg peak corresponding to the static stripe order in magnetic fields and that the Zn substitution tends to stabilize the static order, it is concluded that the suppression of the thermal conductivity in magnetic fields is attributed to the development of the static stripe order. The temperature dependence of the relative reduction of the thermal conductivity in $`x=`$ 0.115 and 0.13 may indicate the temperature dependence of the intensity of the magnetic Bragg peak corresponding to the static stripe order in magnetic fields, as in the case of $`x=`$ 0.10. Moreover, it has been found that $`T_\kappa `$ is directly correlated with $`T_\mathrm{c}`$ rather than $`T_\mathrm{M}`$(0) and $`T_\mu `$, so that $`T_\kappa `$ is regarded as the temperature below which vortex cores pin the dynamical stripes of spins and holes. The present results suggest that the field-induced magnetic order in La2-xSrxCuO4 originates from the pinning of the dynamical stripes of spins and holes by vortex cores.

###### Acknowledgements.

We are grateful to Professor T. Tohyama and Doctor M. Fujita for the useful discussion. The thermal conductivity measurements were performed at the High Field Laboratory for Superconducting Materials, Institute for Materials Research, Tohoku University. This work was supported by a Grant-in Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan. One of the authors (K. K.) was supported by the Japan Society for the Promotion of Science.