text:
# The structure and stability of molecular cloud cores in external radiation fields

## 1 Introduction

It has long been presumed that during the formation of a star, there is an intermediate phase in which the “protostar” is at least approximately in a state of hydrostatic equilibrium or magneto–hydrostatic equilibrium (see e.g. Shu et al. 1987). While this idea in origin was merely based on plausibility arguments, it has received support from the discovery that in least some nearby molecular clouds, one finds embedded “cores” of higher density than the surroundings where the observed linewidths are thermally dominated. That is to say, while there may be local subsonic turbulence, there is no evidence for collapse onto a point mass and it appears that thermal pressure is capable of balancing gravity. On the other hand, comparisons of the gravitational, thermal, magnetic and turbulent energies of such cores show that all these quantities are equal to within the (considerable) uncertainties (Myers & Goodman 1988a, b; Myers et al. 1991; Crutcher 1999). These data suggest though they do not prove the existence of equilibrium structures which are an intermediate state in the evolution of a prestellar core.

This picture has received more observational support with the advent of high quality maps of the millimeter continuum emission of the dust grains within such cores as well as the possibility to study their extinction in the near and mid infrared (Johnstone et al. 2000, 2001; Bacmann et al. 2000; Lada et al. 1994). These have allowed a more unbiased view to be obtained of the density distribution in prestellar cores. In particular, they have shown that the early molecular line data was highly biased because above a critical density of $`5\times 10^4`$ cm-3 (see Tafalla et al. 2002; Caselli et al. 2002a, b; Kramer et al. 1999) most molecular tracers including CO condense out onto dust grain surfaces. This depletion has the consequence that in molecular lines one sees mainly a lower density outer shell, whereas the dust emission (or absorption) offers a more unbiased view of the density distribution.

The results from the millimeter continuum and infrared absorption studies have been compared with a variety of theoretical models of hydrostatic cores (Bacmann et al. 2000). One result of such studies has been evidence for a “flattening” in the density distribution for radii below a critical value of $`r_{\mathrm{cr}}2000`$-8000 AU implying a steeper fall off in radius above $`r_{\mathrm{cr}}`$. Thus for example, Bacmann et al. (2000) model L1544 with a roughly uniform density inside 1900 AU but a rapid fall off outside this radius. Structures whose support against gravitational collapse is mainly due to the magnetic field are plausible both because the observed dust continuum maps show large departures from spherical symmetry and because of time scale arguments. Collapse on a free fall time scale would produce a larger star formation rate than that observed.

There are however some cores where thermal pressure may dominate magnetic pressure and spherical symmetry may be a good assumption. Particularly worthy of note is B68 where Alves et al. (2001) have recently demonstrated that the density distribution derived from their NIR measurements is consistent with a purely thermally supported hydrostatic model. In fact, they find results consistent with the equilibrium structures discussed by Ebert (1955) and Bonnor (1956), which we will refer to in the following as Bonnor-Ebert spheres. To what extent B68 is exceptional is presently unclear but we note that a large number of cores in the Ophiuchus and Orion clouds (Johnstone et al. 2000, 2001) appear to have characteristics compatible with Bonnor-Ebert spheres at temperatures of 15–30 K.

This accord between theoretical expectation and observation suggests that a study of the theoretical assumptions may be worthwhile. One of these assumptions has been that of isothermality which is often based on the concept of “low optical depth for the cooling transitions”. In reality, the gas in prestellar cores is thought to be cooled mainly by optically thick CO transitions (Goldsmith 2001) although, as mentioned above, the CO seems to disappear at high densities. Moreover, the typical density for which dust grains and gas become thermally coupled is roughly of the same order as that observed in prestellar cores (see, e.g. Krügel & Walmsley 1984). Thus, we decided that a new look at the gas temperature distribution to be expected in such cores seemed warranted.

In this paper we replace the assumption of an isothermal gas with the more realistic condition of thermal balance in the gas, and we evaluate the consequences of an external heating source (the interstellar radiation field) on the structure of the cloud (density and temperature profiles) and its stability properties. Recent studies of the dust temperature distribution in such objects (Evans et al. 2001; Zucconi et al. 2001) have shown that the dust temperature typically falls by a factor of 2 from edge to center. It seems reasonable to ask how the gas temperature will react in such circumstances. One might also ask whether the density distribution in hydrostatic equilibrium will depart appreciably from that expected under the isothermal assumption. Will the expected density contrast differ for example from that expected for a Bonnor-Ebert sphere when one calculates the gas temperature in self–consistent fashion? This article represents an attempt to answer such questions.

The outline of this paper is as follows. In Sect. 2 we give a brief introduction to the theory of structures in hydrostatic equilibrium including the results for an isothermal equilibrium Bonnor-Ebert sphere. In Sect. 3 we discuss the input to our calculations and the simplifications which we have made. In Sect. 4, we present our results for the gas temperature distribution in two model cores whose density distribution has been assumed similar to that observed in L1544 and B68. In Sect. 5 we present our results for non-isothermal hydrostatic equilibria for a variety of conditions. Here we show among other things that the equilibria obtained depend relatively sensitively on the external radiation field. In Sect. 6 we discuss the observational consequences of our results, and compare the properties of our model clouds to isentropic polytropes. In Sect. 7 we summarize our conclusions.

## 2 Theoretical background

An isothermal gas sphere embedded in an external medium of given pressure $`p_{\mathrm{ext}}`$ has a critical Bonnor-Ebert mass above which no state of equilibrium can exist (Ebert 1955; Bonnor 1956; McCrea 1957),

$`M_{\mathrm{BE}}1.182{\displaystyle \frac{a^4}{\sqrt{G^3p_{\mathrm{ext}}}}}`$ (1)
$`2.6\left({\displaystyle \frac{T_\mathrm{g}}{10\text{K}}}\right)^2\left({\displaystyle \frac{p_{\mathrm{ext}}}{2\times 10^4\text{K}\text{cm}^3}}\right)^{1/2}M_{},`$
where $`a=(kT_\mathrm{g}/\mu m_\mathrm{H})^{1/2}`$ is the sound speed in the gas, $`T_\mathrm{g}`$ the gas temperature and $`\mu `$ the mean molecular weight (see Fig. 2).

Below this critical value, the equation of hydrostatic equilibrium admits single or multiple solutions, characterized by different degrees of density concentration. It has been shown by Bonnor (1956) and Ebert (1957) that equilibrium configurations with center-to-boundary density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}<13.98`$ are stable, whereas equilibrium solutions with $`\rho _\mathrm{c}/\rho _\mathrm{b}>13.98`$ are unstable to radial collapse. Stable isothermal equilibria can also be characterized by the condition

$$\xi _{\mathrm{max}}\frac{R}{a}\sqrt{4\pi G\rho _\mathrm{c}}<6.451,$$
(2)
where $`R`$ is the cloud’s radius.

Recently, Lombardi & Bertin (2001) have re-examined the classical Bonnor-Ebert problem relaxing the assumption of spherical shape of the cloud. The result is that for non-spherical isothermal clouds, the critical Bonnor-Ebert mass is larger than the value given by Eq. (1); a sphercal shape is more prone to gravitational instability than clouds of other shapes. The density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ for marginally stable non-spherical clouds is different from the Bonnor-Ebert value 13.98 for spherical clouds, but, remarkably, the condition for stability, when expressed in terms of the ratio of the average density $`\overline{\rho }`$ to the boundary density $`\rho _\mathrm{b}`$, i.e. $`\overline{\rho }/\rho _\mathrm{b}<2.465`$, is independent of cloud shape.

To express the condition for stability ($`M<M_{\mathrm{BE}}`$ and $`\rho _\mathrm{c}/\rho _\mathrm{b}<13.98`$) in terms of observable quantities one can eliminate the external pressure writing $`p_{\mathrm{ext}}=a^2\rho _\mathrm{b}>13.98^1a^2\rho _\mathrm{c}`$, obtaining, for a given cloud’s mass $`M`$,

$$\rho _\mathrm{c}<19.53\frac{a^6}{G^3M^2},$$
(3)
that, for a cloud of molecular hydrogen with mean molecular weight $`\mu =2.33`$ becomes

$$n_\mathrm{c}<1.6\times 10^5\left(\frac{T_\mathrm{g}}{10\mathrm{K}}\right)^3\left(\frac{M}{M_{}}\right)^2\text{cm}\text{-3}.$$
(4)
Thus, stable isothermal equilibria are only possible provided the central density is relatively low (in the absence of magnetic fields): this causes a conflict with observations if one attempts to model observed molecular cores with such equilibrium configurations. For example, if dense cores with $`M>4`$ $`M_{}`$ and $`T_\mathrm{g}=10`$ K are modeled as stable Bonnor-Ebert spheres, their density according to Eq. (4) will fall below the critical value of $`10^4`$ $`\text{cm}^3`$ reqired to excite the molecular lines observed in these regions. Moreover, the central densities inferred in many pre–stellar cores based on the observed mm–submm dust emission (e.g. Ward–Thompson et al. 1999) is often above $`10^5`$ $`\text{cm}^3`$ and hence for gas temperatures of $`10`$ K, only cores of mass below 1 $`M_{}`$ can be in stable equilibrium.

Many of the observed objects however have higher masses (see the discussion of Bacmann et al. (2000) and also column densities of up to $`10^{23}`$ $`\text{cm}^2`$ in molecular hydrogen or 100 visual magnitudes of extinction. One can express the stability condition for Bonnor-Ebert spheres as a condition on the cloud’s total (from edge to edge) extinction through the center $`A_{\mathrm{V},\mathrm{tot}}=N(\mathrm{H}_2)/10^{21}\mathrm{cm}^2`$, obtaining

$$A_{\mathrm{V},\mathrm{tot}}\mathrm{¡}\mathrm{}29\left(\frac{T_\mathrm{g}}{10\mathrm{K}}\right)^2\left(\frac{M}{M_{}}\right)^1\text{mag}.$$
(5)
This is clearly much less than often observed and provides a strong argument for the support of many observed structures being due to forces other than pure thermal pressure (though there are also cases of thermal support, e.g. Alves et al. 2001; Johnstone et al. 2000, 2001). In the following, we consider “observed cores” as having a minimum $`A_\mathrm{V}`$ of 10 mag., and thus exclude processes such as UV photoelectric effect from consideration. We do not moreover in this work compute models including magnetic fields and turbulence, but briefly mention in the next Sections how these agent of support modify the picture outlined above.

In summary, non-magnetic stable isothermal cores represent only moderate enhancements (a factor $`2.5`$ in the average density) of the ambient gas density, quite independent on cloud shape and the turbulent field. This limited density range represents a serious limitation when modeling molecular cloud cores as non-magnetized hydrostatic equilibria. For example Boland & de Jong (1984) were able to reproduce column density, brightness temperatures and abundances of several chemical species measured in dark clouds like L134, L183 and TMC-1 only with models characterized by density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}2000`$. Hasegawa (1988) however verified that these models were well into the unstable regime. The same criticism was made by Chièze & Pineau des Fôrets (1987) to the hydrostatic models computed by Falgarone & Puget (1985). Hasegawa (1988) concluded that stable spherical, hydrostatic models able to satisfy all observational constraints on extinction, external pressure, mass, and chemical composition of molecular cloud cores, could no be constructed. In this respect, it is interesting to recall that when molecular cloud cores are fitted by Bonnor-Ebert spheres, the best fit is obtained for models that are in the unstable regime: $`\xi _{\mathrm{max}}=6.9`$ for B68 (Alves et al. 2001), $`\xi _{\mathrm{max}}=12.5`$ for B335 (Harvey et al. 2001, although B335 is not a starless core).

### 2.1 Magnetic fields

The additional support provided by a large-scale magnetic fields can be evaluated by a virial-theorem analysis or with the help of detailed calculations (see e.g. McKee et al. 1987, for a review). The maximum value of a cloud mass for the existence of stable equilibria is approximately equal to

$$M_{\mathrm{cr}}\left[1\left(\frac{0.13}{\lambda }\right)^2\right]^{3/2}M_{\mathrm{BE}},$$
(6)
where $`\lambda `$ is the magnetic mass-to-flux ratio of the cloud in units of $`G^{1/2}`$ (Mouschovias & Spitzer 1976; Tomisaka et al. 1988). For magnetically subcritical clouds ($`\lambda \mathrm{¡}\mathrm{}0.13`$), no amount of external compression can induce the collapse of the cloud (as long as the magnetic field remains frozen in the matter). Magnetically supercritical clouds ($`\lambda \mathrm{¿}\mathrm{}0.13`$) can exist in stable equilibrium for masses larger than the Bonnor-Ebert critical mass, and can reach values of $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ much higher than $`13.98`$ (up to $`100`$, see Fig. 4 and 6 of Tomisaka et al. 1988).

### 2.2 Non-uniform turbulence and temperature gradients

If the turbulent motions of the gas are characterized by a uniform (three-dimensional) mean square velocity $`𝐯_{\mathrm{turb}}^2`$, the effects of a turbulent pressure, in addition to the thermal pressure, can be trivially incorporated in the classic Bonnor-Ebert picture summarized above. The condition for the existence of stable equilibria is the same expressed by Eq. (1) with the sound speed $`a^2`$ replaced by an “effective” sound speed $`a_{\mathrm{eff}}^2=a^2+𝐯_{\mathrm{turb}}^2/3`$, allowing the existence of stable equilibria for values of the central density (or peak extinction) larger than for clouds supported by thermal pressure alone. Notice, however, that (as in the purely thermal case) marginally stable clouds are characterized by a density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}=13.98`$ or $`\overline{\rho }/\rho _\mathrm{b}=2.465`$, where $`\overline{\rho }`$ is the average density, irrespectively on the cloud’s effective temperature.

One should remember also that the observed line widths (which in some cases are close to being thermal) place limits on possible values of $`𝐯_{\mathrm{turb}}`$. In dense cores, as well as in giant molecular clouds, turbulent linewidths are often observed to increase at larger scales (see discussion in Barranco & Goodman 1998; Goodman et al. 1998), suggesting that these objects can be modeled as negative-index polytropes (de Jong et al. 1980; Dickman & Clemens 1983; Maloney 1988). This is easily understood comparing the polytropic equation of state $`P\rho ^{1+1/n}`$ with the expression for the thermal pressure of a gas, $`P\rho T`$, implying $`T\rho ^{1/n}`$. Since the density is always a decreasing function of radius, negative index polytropes represent equilibrium configurations in which the effective temperature increases outward. However, from the point of view of the stability properties, our models show a behavior similar to polytropes with $`n1`$. We will return to this point in Sect. 6.2.

The logatropic equation of state analyzed by McLaughlin & Pudritz (1996) to incorporate the contribution of on non-thermal motions to the support of molecular clouds (and cloud cores) makes possible the existence of stable equilibria for large density constrasts ($`\rho _\mathrm{c}/\rho _\mathrm{b}100`$) but the resulting configurations are much more extended than Bonnor-Ebert spheres and are characterized by a maximum ratio $`\overline{\rho }/\rho _\mathrm{b}=3/2`$, even smaller than for critical Bonnor-Ebert spheres. Indeed, the observed column density distributions (Bacmann et al. 2000) of cores do not fit well to logatropes.

## 3 Model assumptions

Following the approach outlined by Goldsmith (2001), we compute the gas thermal equilibrium assuming that the dust temperature is unaffected by collisions with gas particles (see Evans et al. 2001, for estimates of the validity of this). Moreover, we have used the analytical fit of Zucconi et al. (2001) to determine dust temperature as a function of radius in our spherically symmetric model cores. The justification for this is simplicity. We are interested here in understanding qualitatively to what extent departures from isothermality may affect the properties of pre-stellar cores. High precision in determining the gas cooling rates and the dust temperature distribution is unlikely to fundamentally affect our results.

### 3.1 Dust and gas temperature

Thus we determine the gas temperature $`T_\mathrm{g}`$ by solving the equation of thermal balance of the gas in the presence of dust grains heated by the external radiation field,

$$\mathrm{\Gamma }_{\mathrm{cr}}\mathrm{\Lambda }_\mathrm{g}\mathrm{\Lambda }_{\mathrm{gd}}=0,$$
(7)
where $`\mathrm{\Gamma }_{\mathrm{cr}}`$ is 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. We follow Goldsmith (2001) and adopt a cosmic-ray heating rate

$$\mathrm{\Gamma }_{\mathrm{cr}}=10^{27}\left(\frac{\zeta }{3\times 10^{17}\text{s}\text{-1}}\right)\left[\frac{n(\mathrm{H}_2)}{\mathrm{cm}^3}\right]\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1.$$
(8)
The gas cooling rate $`\mathrm{\Lambda }_\mathrm{g}`$ for typical conditions of molecular cloud cores has been computed by Goldsmith (2001) in the LVG approximation (assuming a velocity gradient of 1 km s-1 pc-1) for various degrees of depletion of molecules from the gas phase onto grain surfaces. We adopt the parametrization of Goldsmith (2001)

$$\mathrm{\Lambda }_\mathrm{g}=\alpha \left(\frac{T_\mathrm{g}}{10\mathrm{K}}\right)^\beta \mathrm{erg}\mathrm{cm}^3\mathrm{s}^1,$$
(9)
where $`\alpha `$ and $`\beta `$ are parameters that depend on the H2 density and the molecular depletion factors. When not explicitly noted, we assume that molecular abundances in the gas have their standard (undepleted) values. As for the the gas-dust energy transfer rate we adopt the expression given by Goldsmith (2001) based on calculations of Burke & Hollenbach (1983),

$`\mathrm{\Lambda }_{\mathrm{gd}}=2\times 10^{33}\left({\displaystyle \frac{n}{\mathrm{cm}^3}}\right)^2\left({\displaystyle \frac{T_\mathrm{g}T_\mathrm{d}}{\mathrm{K}}}\right)`$ (10)
$`\times \left({\displaystyle \frac{T_\mathrm{g}}{10\mathrm{K}}}\right)^{1/2}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1.`$
To determine the dust temperature $`T_\mathrm{d}`$ we ignore the gas-dust coupling and we solve the equation of thermal balance of the dust

$$\mathrm{\Gamma }_{\mathrm{ext}}\mathrm{\Lambda }_\mathrm{d}=0,$$
(11)
where $`\mathrm{\Gamma }_{\mathrm{ext}}`$ is the dust heating rate by the external radiation field and $`\mathrm{\Lambda }_\mathrm{d}`$ the dust cooling rate by radiation. In practice, we use the analytical solution of Zucconi et al. (2001) to compute $`T_\mathrm{d}`$ as function of extinction from the cloud’s boundary.

As in Zucconi et al. (2001), we have adopted as a standard the interstellar radiation field (ISRF) in the solar neighborhood given by Mathis et al. (1983) and Black (1994). It consists of four components: the V-NIR component, peaking at $`\lambda 1`$ $`\mu \text{m}`$, the MIR and FIR components, peaking at $`\lambda 100`$ and 140 $`\mu \text{m}`$ respectively, and the cosmic background radiation, peaking at $`\lambda 1`$ mm. We allow for variations of the ISRF, assuming, for simplicity, that the spectral shapes of these components are the same as the local ISRF, and we scale the intensity of the V-NIR, MIR, and FIR components by a factor $`G_0`$ (with $`G_0=1`$ for the local ISRF).

The existence of a cloud radius leads to the necessity of an external pressure to keep the system in equilibrium. We embed our model cloud cores in a spherical envelope that represents the ambient molecular cloud. We do not need to specify the physical properties of this envelope; we just assume that it provides the needed external pressure and a shielding from the ISRF directly incident on the cloud cores. We set the latter at the value $`A_\mathrm{V}^{\mathrm{env}}=1`$ (McKee 1999), sufficient to absorb the UV and most of the V-NIR radiation. We define $`A_{\mathrm{V},\mathrm{c}}`$ to be the centre-to-edge extinction of our model cloud cores. The range of masses and densities of interest here corresponds to the range $`5\mathrm{¡}\mathrm{}A_{\mathrm{V},\mathrm{c}}\mathrm{¡}\mathrm{}100`$, where the central dust temperature is determined mostly by optically thick absorption of V-NIR radiation and optically thin absorbtion of MIR and FIR radiation.

The fundamental parameters of our models are the external pressure $`p_{\mathrm{ext}}`$ and the intensity of the external radiation field $`G_0`$. In this paper, we consider independent variations of these two quantities, and we analyse separately the dependence of the critical mass for increasing $`p_{\mathrm{ext}}`$ and fixed $`G_0`$ (Sect. 5, Fig. 2), and the effect of varying $`G_0`$ for a given value of $`p_{\mathrm{ext}}`$ (Sect. 6, Fig. 6).

### 3.2 Effects of molecular depletion

We include in our models the effect of molecular depletion onto dust grains, parametrized for simplicity by the single parameter $`f_\mathrm{d}`$ (depletion factor), defined as the ratio of the “standard” abundance of CO isotopes and CS to the actual abundances. The corresponding abundances of other molecular species are given in Table 3 of Goldsmith (2001). Gas-phase cooling rates for several values of $`f_\mathrm{d}`$ in the range 1–100 were computed and tabulated by Goldsmith (2001) for H2 density from $`10^3`$ cm-3 to $`10^6`$ cm-3. We assume the dependence of $`f_\mathrm{d}`$ on density given by Tafalla et al. (2002),

$$f_\mathrm{d}=\mathrm{exp}\left(\frac{n}{n_{\mathrm{dep}}}\right),$$
(12)
where the critical density for CO depletion is taken to be $`n_{\mathrm{dep}}=5.5\times 10^4`$ $`\text{cm}^3`$.

Although we are mainly interested in applications to molecular cloud cores of moderate to high extinction, where molecular depletion is effective only in a region of limited spatial extent, one should keep in mind that large fraction of the cloud’s mass is contained in the outer parts (roughly $`M(r)r`$), and therefore the overall stability of the cloud may be modified by the thermal properties of these external layers.

In our models we find in general that the effects of depletion on the temperature of the gas are significant only for densities $`n\mathrm{¡}\mathrm{}10^5`$ cm-3. At higher density, gas-dust coupling overwhelms the reduction of the cooling rate of the gas, and the gas temperature becomes insensitive to depletion. This is in qualitative agreement with the results of Goldsmith (2001). We conclude with Goldsmith (2001) that molecular depletion in the central, dense regions of molecular cloud cores ($`n\mathrm{¿}\mathrm{}10^{4.5}`$$`10^5`$ cm-3) is not associated with an increase in the gas temperature.

In the extreme case where CO (and other molecules containing C,N,O,S, etc.) depletes out onto grains causing Eq. (7) to reduce to $`\mathrm{\Gamma }_{\mathrm{cr}}=\mathrm{\Lambda }_{\mathrm{gd}}`$. Then using Eq. (8) and (10), one finds

$`{\displaystyle \frac{T_\mathrm{g}T_\mathrm{d}}{T_\mathrm{g}}}`$ (13)
$`5\left({\displaystyle \frac{\zeta }{3\times 10^{17}\text{s}^1}}\right)\left({\displaystyle \frac{n}{10^4\text{cm}^3}}\right)^1\left({\displaystyle \frac{T_\mathrm{g}}{10\text{K}}}\right)^{1.5},`$
This shows that in the extreme case where all heavy species deplete out, coupling between gas and dust temperatures will be assured for densities larger than $`5\times 10^4`$ $`\text{cm}^3`$ if the cosmic ray ionization rate is standard or lower. Since it is essentially at densities above this value that CO is observed to deplete out, we conclude that a rough equality between gas and dust temperatures in the depleted region is a reasonable assumption. In this context, it is interesting to note that Bergin et al. (2002) suggest that even such relatively volatile species as molecular nitrogen have depleted out in the central core of B68.

## 4 Temperature distributions for L1544 and B68

Before considering what one expects for theoretical equilibrium cores, it is useful to consider what one predicts for the temperature variation within an object whose density structure is similar to those observed. Specifically, we will consider the case of L1544 (Caselli et al. 2002a, b) and B68 Alves et al. (2001), two well studied prestellar cores.

To model L1544 we assume $`G_0=1`$ and the density structure of Tafalla et al. (2002), which, in turn, is based on the 1.3mm dust emission and absorption maps of Bacmann et al. (2000). Thus, we assume the H2 density distribution $`n(r)`$ to be given by

$$n(r)=\frac{n_\mathrm{c}}{1+(r/r_0)^{2.5}},$$
(14)
where the central $`\mathrm{H}_2`$ density is $`n_\mathrm{c}=1.4\times 10^6`$ $`\text{cm}^3`$ and $`r_0=0.014`$ pc. It is worth noting here that we have on purpose taken here a somewhat extreme case with large density contrast relative to the surroundings and large central column density ($`6.0\times 10^{22}`$ $`\text{cm}^2`$ corresponding to 60 mag. of extinction, from Bacmann et al. 2000). One should also be aware that the aspect ratio of the core based on the dust emission contours is $`0.6`$ and hence spherical symmetry is a crude assumption. Nevertheless, the case is of interest in that it is a model which approximates the actually observed density distribution and in which, one expects that the central temperature will be relatively low (due to the large extinction).

For B68 we assume the density structure of a Bonnor-Ebert sphere with the physical parameters given by Alves et al. (2001). The central density $`n_\mathrm{c}=2.4\times 10^5`$ cm-3 is a factor $`6`$ lower than in the case of L 1554, and consequently the amount of CO depletion at the center of B68 is about two orders of magnitude lower than at the center of L1544. The intensity of the IRSF incident on B68 can be estimated as follows. We first estimate the external UV radiation field based on the 90 $`\mu \text{m}`$ data of Ward-Thompson et al. (2002) and the conversion between FIR intensity and incident field found by Boulanger & Perault (1988). Another approach is to use the 7$`\mu \text{m}`$ intensity from Bacmann et al. (2000) and the results of Boulanger et al. (1996) to infer $`G_0`$. Combining these, we conclude that the radiation field incident on B68 is $`2.5`$ times larger that the standard ISRF.

In Fig. 1, we show the inferred gas temperature dependence for L1544 and B68, computed with cosmic ray ionization rate $`\zeta =1.3\times 10^{17}`$ s-1 and accomodation coefficient $`=0.3`$. The assumed CO depletion factor (cut off at a value of 100) is shown for comparison. In the inner part of both cores (densities above $`3\times 10^4`$ $`\text{cm}^3`$), the gas temperature is coupled to that of the dust and thus decreases gradually towards the core center. This coupling is more effective in the case of L1544, which has a higher central density. The gas is slightly hotter than the dust however due to cosmic ray heating and this difference increases as the density decreases. The CO depletion appears to have only small effects on the gas temperature profile because while at high density, gas–dust coupling dominates, at densities below $`10^5`$ $`\text{cm}^3`$, the depletion is not very large (and the optically thick nature of CO cooling causes the dependence of CO cooling on abundance to be minor). We draw the conclusion from this study that in the depleted core nucleus (densities above $`10^5`$ $`\text{cm}^3`$) where many molecular species may be unobservable because they are in solid form, the gas temperature may be expected to be close to that of the dust. In outer regions where the density is below $`10^5`$ $`\text{cm}^3`$ and where shielding from the external radiation field is much less, the gas and dust temperatures become uncoupled and one may find gas temperatures lower than that of the dust.

It is worth pointing out that the scenario depicted in Fig. 1 yields gas temperatures in conflict with observation. Tafalla et al. (2002) find based on their NH3 observations a temperature of 8.7 K throughout the L1544 core with no evidence for density gradients (40″ HPBW equivalent to 0.03 parsec), whereas Bourke et al. (1995) give $`T_\mathrm{g}=16`$ K for B68, also derived from NH3 data. These are in rough agreement with the values in Fig. 1 though the L1544 data do not show evidence for increasing temperature with radius as predicted.

Another conclusion that one can draw from Fig. 1 is that the gas temperature and hence pressure in the high density core nucleus depends somewhat on the external radiation field even for cores of extremely high visual extinction. Thus one can expect cores in Taurus for example (where the incident radiation field appears to be average) to have different characteristics than Ophiuchus cores where the radiation field appears in general to be an order of magnitude higher. Of course, there are other effects such as the external pressure which may differ in the two cases also but we stress that the radiation field alone may cause core temperatures and hence accretion rates to be higher in Ophiuchus than in Taurus.

## 5 Self consistent hydrostatic equilibrium results

The equation of hydrostatic equilibrium for a spherically symmetric cloud is

$$\frac{dp}{dr}+\rho \frac{dV}{dr}=0,$$
(15)
where $`p`$ is the gas pressure and $`V`$ the gravitational potential. The latter satisfies Poisson’s equation, which in spherical symmetry reads

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dV}{dr}\right)=4\pi G\rho .$$
(16)
The gas pressure $`p`$, neglecting the contribution of turbulent pressure, is

$$p=\frac{\rho }{\mu m_H}kT_\mathrm{g},$$
(17)
where $`T_\mathrm{g}`$ is the gas temperature, computed as described in Sect. 3, and $`\mu =2.33`$ the mean molecular weight.

The behavior of our model clouds with respect to an external compression is qualitatively similar to that of isothermal spheres. Fig. 2 shows the locus in the radius-external pressure diagram of cloud models with $`M=1`$ $`M_{}`$, for different values of the ISRF intensity $`G_0`$. For large radii the cloud follows closely the $`T_\mathrm{g}=10`$ K Bonnor-Ebert isotherm, independent of the radiation field (the effect of gas-dust coupling is negligible for these diffuse configurations). At higher values of the external pressure the cloud radius is smaller and the central (and average) density becomes larger. As a result, the central gas temperature increases and the sequence of equilibria departs from the Bonnor-Ebert isotherm. Eventually, a critical value of the pressure is reached above which no equilibria exist, and the curve turns down. Past this critical point, equilibria are characterized by an increasing central density and density contrast, and are unstable (see following Section). For clarity, in Fig. 2 we truncate the curves of equilibrium slightly after the critical point.

Fig. 3 shows the density of H2, the temperature of the gas and the temperature of the dust as function of radius for the marginally stable clouds of mass $`M=1`$ $`M_{}`$ with $`G_0=1`$ (central H2 density $`n_\mathrm{c}=2.7\times 10^5`$ cm-3). The density profile is very close to that of an isothermal sphere with the same values of the central density, shown for comparison in Fig. 3 for temperatures equal to the central and boundary temperature of the model. The gas temperature decreases from the outer regions to the center of the cloud. For marginally stable clouds with $`M=1`$ $`M_{}`$the central density is sufficient to ensure a modest gas-dust coupling, but the gas temperature is mainly determined by the balance of cosmic-ray heating and gas cooling.

Fig. 4 shows the H2 density, the gas and dust temperatures and the extinction profile of the marginally stable cloud in thermal equilibrium with an ISRF with $`G_0=100`$ and bounded by an external pressure $`p_{\mathrm{ext}}=1\times 10^6`$ K cm-3. These roughly correspond to the structures observed by Johnstone et al. (2001) in the $`\rho `$Oph and Orion clouds. For comparison, the figure also shows the density profiles of marginally stable Bonnor-Ebert spheres with the same central density and uniform temperature equal to the gas temperature at the center and at the boundary of the cloud.

In Tables 1 and 2 we list the values of central H2 density, gas and dust temperature, density contrast, and optical extinction at the center for marginally stable clouds with $`M=1`$ and 5 $`M_{}`$and various scaling factors for the ISRF. In these models we have varied $`p_{\mathrm{ext}}`$ until we find the marginally stable model.

A comparison of the results for the two different masses is instructive. Small mass clouds ($`M1`$ $`M_{}`$) exposed to the average ISRF behave like isothermal spheres with $`T_\mathrm{g}10`$ K, and are characterized by a maximum density contrast and critical central density very close to those of the corresponding Bonnor-Ebert values. Increasing the ISRF ($`G_0=10`$) has the effect of raising the dust temperature, and of increasing slightly the gas temperature in the central regions. This increase in the central gas temperature causes a higher central density in the marginally stable core and this, in turn, increases the efficiency of dust-gas coupling. A further increase of the ISRF ($`G_0=100`$) establishes the approximate equality of the gas and dust temperatures at the center. Thus, the fact that small mass clouds are dynamically stable for larger values of the central density than their more massive counterparts (see Eq. 4), has the consequence that they are “stabilized” by an external ISRF simply because they are able to become warmer if the dust does. Their gas temperature profiles, however, remain quite uniform, and therefore the maximum density contrast cannot deviate significantly from the Bonnor-Ebert value $`\rho _\mathrm{c}/\rho _\mathrm{b}13.98`$, as shown in Table 1.

A different behaviour characterizes more massive clouds ($`M5`$ $`M_{}`$), as they become dynamically unstable before an efficient gas-dust coupling can be established, for any value of the intensity of the ISRF (see Table 2). Their central extinction is low ($`A_{\mathrm{V},\mathrm{c}}3`$–5), and hence they do not correspond to “real observed cores”. Nevertheless, it is interesting that they can develop significant inward temperature gradients, because molecular cooling is less efficient than cosmic-ray heating in the central regions, characterized by densities $`n10^4`$$`10^5`$ cm-3 (where $`\mathrm{\Lambda }_\mathrm{g}n^{0.5}`$) than in the outer parts where $`n10^3`$ (and $`\mathrm{\Lambda }_\mathrm{g}n`$, as cosmic-ray heating). An inward temperature gradient stabilizes the cloud (see Sect. 6), and allows large deviations in the maximum density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}`$, up to 30–40, as shown in Table 2. Notice that any molecular depletion effect, ignored in these calculation, would make the central temperature even higher, and thus increase the slope of the temperature gradient.

In Fig. 5 we summarize our results for the maximum density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ allowed for clouds in the mass range $`M=1`$–10 $`M_{}`$. Again $`p_{\mathrm{ext}}`$ has been varied until the marginally stable model is found. Deviations of $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ from the Bonnor-Ebert value 13.98 indicate deviations from a condition of uniform gas temperature in the cloud. As anticipated in the above discussion, the largest deviations from isothermality are obtained for masses $`M4`$–6 $`M_{}`$. Marginally stable clouds of small mass ($`M\mathrm{¡}\mathrm{}2`$ $`M_{}`$) are characterized by relatively large values of the central density, and by a sufficiently good coupling of gas and dust: they are uniformly warmer or colder according to the intensity of the external ISRF, but because of this coupling they do not develop significant temperature gradients, and therefore behave essentially as isothermal spheres. Marginally stable high-mass clouds ($`M\mathrm{¿}\mathrm{}8`$ $`M_{}`$) also are similar to isothermal spheres, but for a different reason: they are characterized by relatively low values of the central density ($`n_\mathrm{c}10^3`$ cm-3) and also have low $`A_{\mathrm{V},\mathrm{c}}`$ (unlike observed cores), which implies an uniform value of the gas temperature since both cosmic-ray heating and radiative gas cooling depend approximately linearly on density (and gas-dust coupling plays no role at these densities). In the intermediate mass range, both because of the different dependence of heating and cooling rates from density and partial gas-dust coupling, the deviations from uniform gas temperature are larger, and the maximum density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ has a peak for $`M4`$–5 $`M_{}`$. Fig. 5 also shows, for comparison, the values of the density contrast for marginally stable polytropic spheres of index $`n`$, discussed in the next Section.

## 6 The effect of the external radiation field

The presence of an external ISRF affects the stability of clouds bounded by an external pressure. In the first place, it makes the dust (and therefore the gas) warmer on average, thus increasing the critical mass for gravitational instability (proportional to $`T_\mathrm{g}^2`$); in the second place, it enforces a temperature gradient in the cloud which modifies the value of the classic Bonnor-Ebert critical mass.

The latter point is best investigated by the so-called static method originally devised by Zeldovich (1963) for white dwarfs and neutron stars (see Tassoul 1978). This method allows one to perform a stability analysis for axisymmetric motions on the basis of the properties of the equilibrium models only. The static criterion asserts that a turnover in the $`M(\rho _\mathrm{c})`$ curves at fixed $`P_{\mathrm{ext}}`$, occurring at say $`\rho _c=\rho _c^{\mathrm{cr}}(p_{\mathrm{ext}})`$, marks the onset of dynamical instability to radial motions. The condition $`(M/\rho _\mathrm{c})_{p_{\mathrm{ext}}}=0`$ is both a sufficient and necessary condition for a stability transition: equilibria with $`\rho _\mathrm{c}<\rho _\mathrm{c}^{\mathrm{cr}}(p_{\mathrm{ext}})`$ are stable, while equilibria with $`\rho _\mathrm{c}>\rho _\mathrm{c}^{\mathrm{cr}}(p_{\mathrm{ext}})`$ are unstable to at least one normal mode.

In Fig. 6 we show how the intensity of the external ISRF may affect the stability of a cloud bounded by a given external pressure. For this particular example, we have chosen $`p_{\mathrm{ext}}=2\times 10^4`$ K cm-3 (McKee 1999). The three curves labelled with $`G_0=1`$, 10 and 100 show the mass versus central H2 density for externally heated clouds, compared with the similar quantities for isothermal spheres. According to the static method, continuous (dashed) lines represent stable (unstable) equilibria. The dotted lines indicates the boundaries of the stability domain for isothermal spheres with arbitrary temperatures (only the case $`T_{\mathrm{gas}}=10`$ K is shown) and non-isothermal spheres in thermal equilibrium with the external ISRF. The figure shows that for a given external pressure the condition of thermal equilibrium allows the existence of stable cloud configurations for values of mass and central density not permitted under the condition of uniform gas temperature.

Although not very large, this effect may resolve in part the problem posed by objects like B68 which are well fitted by slightly unstable Bonnor-Ebert spheres (magnetic fields are of course another possibility). For example, for the three marginally stable non-isothermal clouds of Fig. 6 the central-to-boundary density ratio are $`(\rho _\mathrm{c}/\rho _\mathrm{b})=19`$, 25 and 39, for $`G_0=1`$, 10 and 100, compared with the corresponding ratio 13.98 for isothermal spheres, and 16.5 inferred for B68 (Alves et al. 2001).

### 6.1 Comparison with polytropic models

From the point of view of the stability properties and the maximum density contrast, the models shown in this paper have a behavior similar to polytropes with $`n1`$ (see Fig. 3). Increasing the intensity of the ISRF produces an increase in the dust temperature which is felt by the gas only in the central parts of the cloud, where the coupling of dust and gas is more efficient. As a result, the central gas temperature $`T_{\mathrm{g},\mathrm{c}}`$ increases slightly, whereas the gas temperature at the boundary $`T_{\mathrm{g},\mathrm{b}}`$ remains unchanged ($`T_{\mathrm{g},\mathrm{b}}10`$ K), since for densities $`n10^3`$ cm-3 both the gas cooling rate and the the cosmic-ray heating rate are roughly proportional to $`n`$. Thus, one can extend Eq. (1) to polytropic clouds bounded by the same external pressure and characterized by the same thermal sound speed $`a_\mathrm{b}`$ at the boundary, and write

$$M_{\mathrm{cr}}\alpha _n\frac{a_\mathrm{b}^4}{\sqrt{G^3p_{\mathrm{ext}}}}.$$
(18)
Table 3 shows the value of $`\alpha _n`$ and the density contrast $`\rho _\mathrm{c}/\rho _\mathrm{b}`$ for clouds with different polytropic index $`n`$. For decreasing positive $`n`$, the coefficient $`\alpha _n`$ increases, indicating a larger critical mass than an isothermal sphere bounded by the same external pressure and with the same temperature at the boundary. For $`n=6`$, for example, the critical mass is expected to be a factor $`2.06/1.181.75`$ larger than critical Bonnor-Ebert mass for the same gas temperature at the boundary ($`M_{\mathrm{BE}}2.6`$ $`M_{}`$ for $`T_\mathrm{g}=10`$ K, see Eq. 1).

## 7 Discussion and Conclusions

We have examined the gas temperature distribution to be expected in interstellar pre–protostellar cores heated by the external ISRF. We find that when (as in observed cores), the central density exceeds $`3\times 10^4`$ $`\text{cm}^3`$, there is coupling between the gas and dust temperatures and hence the gas temperature (like the dust) decreases with decreasing radius. At larger radii and smaller densities, the dust and gas decouple and the gas temperature may (for low external pressures as in the Taurus cloud) decrease towards the values expected for heating by galactic cosmic rays of around 10 K. The region where gas and dust temperatures are coupled is somewhat interior to the region where CO is highly depleted. However, CO depletion does not seem greatly to affect the temperature distribution mainly because cooling by gas–grain collisions becomes dominant.

The observed values of the temperature of around 10 K in many pre-stellar cores allow limits to be placed on the cosmic ray ionization rate similar to the standard value of order $`10^{17}`$ s-1 based on the measured cosmic ray flux. The fact that measured gas temperatures in cores show so little spatial variation (Benson & Myers 1989; Tafalla et al. 2002) suggests to us that the observed thermometers mainly trace layers of moderate depletion where gas–grain coupling is not playing an important role.

We have also examined the consequences of such a temperature distribution for the density distribution in hydrostatic equilibrium cores. The changes caused by the “real temperature distribution” are minor and the characteristics of a “marginally stable Bonnor-Ebert sphere” are similar to those in the isothermal case.

An interesting point which emerges from these calculations is that the temperature in the core nucleus is sensitive to the external radiation field (because the core nucleus as a rule is heated by small particle MIR emission from the borders of the surrounding cloud). This is in particular the case for cores with high external pressure (and hence high density) such as those studied by Johnstone et al. (2001, 2000) in the Orion B and $`\rho `$Oph clouds where $`p_{\mathrm{ext}}`$ appears to be above $`10^6`$ $`\text{cm}^3`$K. In these cases, we expect the dust temperature to follow the gas temperature and thus the behavior should be roughly like a negative index polytrope. The external radiation fields are also higher in Orion and Ophiuchus than in Taurus and thus the core masses may also rise somewhat. For constant external pressure, the mass of the marginally stable Bonnor-Ebert sphere rises with about the 0.35 power of the external radiation field. It would be useful to have direct temperature estimates to confirm these expectations.

Another result of this study is that high mass thermally supported cores are incompatible with high column density. Thus it is difficult to imagine observed high column density (greater than $`10^{22}`$ $`\text{cm}^2`$) high mass (greater than 10 $`M_{}`$) cores going through a series of quasi-equilibrium states en route to collapse. Hence we suspect that higher mass cores are either not stable or have other (magnetic) means of support. These other means of support become apparent observationally both because observed cores show large departures from spherical symmetry and also because the observed line widths in most cores can only be explained in terms of turbulence and often of supersonic turbulence.

Thus thermally supported cores may be the exception rather than the rule. It is significant nonetheless that cores in regions such as Taurus where predominantly low mass star formation is taking place have in general line profiles showing a large component of thermal broadening. This suggests that where (as in Taurus), the pivotal state may be a core with predominantly thermal support, only low mass stars are likely to form. In clouds such as Ophiuchus and Orion with higher radiation fields and pressures and with cores having predominantly non–thermal support, higher mass stars may become possible and star formation may take another course.

It is also significant that occasionally, one finds cases like B68 where the data are consistent with hydrostatic equilibrium (marginally stable) and pure thermal support. In fact, we find that when one takes the gas temperature dependence into account, B68 is marginally stable. But irrespective of whether this is true or not, B68 gives every sign of being close to the pivotal state from which protostellar collapse will commence.

###### Acknowledgements.

It is a pleasure to thank A. Natta for constructive criticism and Jens Kauffman for useful comments on the manuscript. We acknowledge financial support from the EC Research Training Network “The Formation and Evolution of Young Stellar Clusters” and grant COFIN-2000. JG acknowledges support from scolarship SFRH/BD/6108/2001 awarded by the Portuguese Ministry of Education.


format: text

text:
# Relaxing in foam

## Abstract

We investigate the linear mechanical response of an aqueous foam, and its relation to the microscopic rearrangement dynamics of the bubble-packing structure. At rest, even though the foam is coarsening, the rheology is demonstrated to be linear. Under flow, shear-induced rearrangements compete with coarsening-induced rearrangements. The macroscopic consequences are captured by a novel rheological method in which a step-strain is superposed on an otherwise steady flow.

Aqueous foams are tightly packed collections of gas bubbles separated by a continuous liquid phase Prud’homme and Khan (1996); Weaire and Hutzler (1999). Like elastic solids, bulk foams resist shear, completely unlike the gases and liquids from which the they are comprised. The origin of this striking behavior is that the bubbles are jammed, unable to flow around one another and explore configuration space under thermal energy. Thus bubbles distort, rather than rearrange, when subjected to small shear deformations. The resulting extra internal gas-liquid surface area costs energy in proportion to the surface tension, and this provides a restoring force. The shear modulus is roughly surface tension divided by bubble radius, depending precisely on the volume fraction, $`\epsilon `$, of the continuous liquid phase Princen (1983); Mason et al. (1995); Saint-Jalmes and Durian (1999). As a foam is made wetter, the bubbles become progressively rounder, and the shear modulus decreases. The elasticity completely vanishes at the point where the bubbles are close-packed spheres. This is one example of unjamming Liu and Nagel (2001).

In this paper we explore other ways to unjam the bubbles in a foam. Each could correspond to a different trajectory in a global jamming phase diagram Liu and Nagel (1998). For example, one could imagine raising the temperature so that $`k_BT`$ is greater than $`\sigma R^2`$, which would allow the bubbles to rearrange thermally like Brownian particles. However, a typical value is $`\sigma R^2/k_B10^{12}`$ K; therefore, foams are athermal far-from-equilibrium systems and raising the temperature is not feasible. But there are at least two other ways to drive the system so that bubble rearrangements occur. One is through coarsening, the diffusion of gas from smaller to larger bubbles. As this proceeds, local stress inhomogeneities repeatedly build up to some theshold and relax by sudden avalanche-like local rearrangements. Such microscopic dynamics have been observed previously via diffusing-wave spectroscopy (DWS) Durian et al. (1991a); Gopal and Durian (1997); Cohen-Addad and Höhler (2001). This raises a string of interesting questions. How do localized coarsening-induced rearrangements affect the linear shear rheology of bulk foams? Is the effect similar to thermally excited dynamic heterogeneities in glassy systems? An entirely different way to induce rearrangments is by shear. Here local stresses also accrue and relax, but in a more correlated manner. These rearrangment dynamics have similarly been studied by DWS Earnshaw and Jafaar (1994); Gopal and Durian (1995, 1999). How does the elastic character of foam vanish as the shear rate is increased? Can we think of the two driving mechanisms as “heating” and ultimately unjamming the sample? To make progress, we study the linear mechanical response of an aqueous foam at very long times and low frequencies. In a novel twist, we also study linear response during uniform shear flow.

Our samples are a commercial aqueous foam, consisting of nearly-spherical polydisperse gas bubbles, 92% by volume, that are tightly packed in an aqueous solution of stearic acid and triethanolamine (Gillette Company, Boston MA). Samples are surrounded by a water bath held at $`25.0^{}`$C and are measured after the foam has aged for approximately 100 minutes. By this time, the average bubble size is approximately 60 microns and is growing reproducibly via coarsening Durian et al. (1991b); stresses due to the loading process have also relaxed. Test durations are sufficiently short that gravitational liquid drainage and bubble coalescence are negligible.

Our measurements are performed with a Paar Physica UDS 200 rheometer, controlling the rotational speed and angular displacement of a solid cylinder whose axis is vertical and concentric with a fixed surrounding cup. The sample cell has an inner radius of 20.0 mm, and a 4.1 mm gap; these dimensions ensure that the foam can be treated as a bulk material with uniform stress. To minimize end-flow effects, the 98 mm long inner cylinder is much shorter than the depth of the surrounding cup yet much longer than the gap width. Wall slip is precluded by coating both the cylinder and cup with a fine-grade sandpaper. Samples are loaded through a 5 mm hole at the bottom of the cup, after lowering the cylinder to the test position. DWS measurements indicate the absence of shear-banding and other secondary flows Gopal and Durian (1999). To re-confirm, we compared with a similarly coated cone-and-plate cell, with a 10 cm radius and $`10^{}`$ cone angle. Identical rheology results were obtained for both cell configurations, implying that the imposed shear deformation is uniform. Diffuse light transmission indicates that bubbles do not burst, and voids do not form, even under high shear Gopal and Durian (1995, 1999).

Effect of coarsening on mechanical response. To quantify linear response, we measure both the complex shear modulus, $`G^{}(\omega )=G^{}(\omega )+iG^{\prime \prime }(\omega )`$ Ferry (1980); Macosko (1994), and the stress relaxation modulus, $`G(t)`$ Ferry (1980); Macosko (1994). Our data are displayed in Fig. 1. Note that $`G^{}(\omega )>G^{\prime \prime }(\omega )`$, meaning that the response is primarily elastic rather than dissipative. The frequency range for our $`G^{}(\omega )`$ data spans almost six decades, from $`2\pi `$ divided by sample age (an absolute minimum below which measurement is not possible) up to a maximum set by limitations of the rheometer. Note that different amplitude strains give the same result, demonstrating absence of wall-slip and other geometry-dependent artifacts and hinting at linearity of response. A more stringent test of linearity is comparison with $`G(t)`$, which should be related to $`G^{}(\omega )`$ by Fourier transform Ferry (1980); Macosko (1994). The time range for our $`G(t)`$ data spans over five decades, from the time needed to achieve the step strain up to the time beyond which stress is zero to within instrumental limits. An empirical fit to the $`G(t)`$ data is shown in Fig. 1(b) by a solid curve; this fit is transformed and plotted in Fig. 1(a) over a frequency range corresponding to the time range of the fit. The agreement is very good, demonstrating conclusively that the sample is linear. This is further supported by an empirical fit to $`G^{}(\omega )`$ data at high-$`\omega `$ and the comparison of its transform with $`G(t)`$ data at short-$`t`$.

Let us now consider the frequency and time dependence of the moduli in Fig. 1. The fit for $`\omega >5`$ rad/s is to the form $`G^{}(\omega )=G_{}(1+\sqrt{i\omega /\omega _n})`$ with free parameters $`G_{}=2300`$ dyne/cm2 and $`\omega _n=156`$ rad/s. The former represents the static shear modulus, roughly surface tension divided by bubble size Princen (1983); Mason et al. (1995); Saint-Jalmes and Durian (1999). The latter represents the effect of nonaffine deformation of the bubbles under shear due to local packing configurations that are strong or weak with respect to the shear direction Liu et al. (1996). This fit, including the parameters, is consistent with the $`G^{}(\omega )`$ data in Ref. Cohen-Addad et al. (1998), where the frequency range ($`0.320`$ rad/s) was too small to fully demonstrate the functional form. According to the theory of Liu et al. (1996), the characteristic frequency is $`\omega _nG_{}/\eta _{\mathrm{}}`$, where the very-high frequency response is $`G^{}(\omega )=i\eta _{\mathrm{}}\omega `$. The numerical prefactor was not predicted; experimentally, it was found to depend on $`\epsilon `$ and was not of order 1. To compare, the value is $`\omega _n600`$ rad/s for a $`\epsilon =0.38`$ emulsion of 0.5 $`\mu `$m oil droplets in water Mason and Weitz (1995); Liu et al. (1996) (NB: by our definition, $`\omega _n`$ may be easily read off a plot by locating where $`G^{}(2\omega _n)=2G_{}`$ and $`G^{\prime \prime }(2\omega _n)=G_{}`$).

The fit to $`G_{}(1+\sqrt{i\omega /\omega _n})`$ fails for $`\omega <5`$ rad/s; the corresponding transform fails for $`t>20`$ s. At longer times, the $`G(t)`$ data decay slowly below $`G_{}`$, almost logarithmically, over a few decades before relaxing more rapidly at around 1000 s. The transform of this final decay corresponds to the peak in $`G^{\prime \prime }(\omega )`$ at $`10^3`$ rad/s. At lower frequencies, $`G^{}(\omega )`$ is unmeasurable; but since the integral of $`G(t)`$ over all time is finite, the very-low frequency behavior is formally $`G^{\prime \prime }(\omega )\omega `$ as required by causality Buzza et al. (1995). Thus the full frequency-dependence of $`G^{}(\omega )`$ for our foam is truly known and well-behaved. This resolves a long-standing puzzle Buzza et al. (1995) raised by earlier measurements Khan et al. (1988); Mason and Weitz (1995); Cohen-Addad et al. (1998) where $`G^{}(\omega )`$ was roughly constant down to the lowest measured frequencies.

All that remains is to understand the origin of the low-$`\omega `$ / long-$`t`$ behavior. We contend that evolution of the foam structure by coarsening is responsible. One clue is that the onset of deviation from the high-$`\omega `$ fits corresponds to the time $`\tau _{oq}=20`$ s given by DWS for the time between coarsening-induced rearrangements at each site. Another clue is that the final decay of $`G(t)`$ and, equivalently, the peak in $`G^{\prime \prime }(\omega )`$ correspond to the sample age. Since coarsening gives power-law growth, it takes of order the sample age for the structure to completely change. It is interesting that coarsening-induced rearrangments relax microscopic coarsening-induced stress inhomogeneities far more quickly than the relaxation of macroscopically-imposed stress. Rather, the cumulative effect of many rearrangements and a change in bubble size is needed to relax global stress. The net result is a rheology that obeys linear response. This is remarkable given that the microscopic relaxation mechanism isn’t thermal motion, but rather evolution. In effect, coarsening unjams the foam, so that at low frequencies the rheology is $`G^{}(\omega )G^{\prime \prime }(\omega )=\eta \omega `$ like an ordinary equilibrium liquid.

Effect of shear on mechanical response. Next we investigate how the application of shear causes a similar loss of elasticity. For this we superimpose a small-amplitude step-strain, $`\mathrm{\Delta }\gamma `$, on top of otherwise steady shear at rate $`\dot{\gamma }`$. The resulting temporary increase in stress defines the transient stress relaxation modulus, $`G(t,\dot{\gamma })=[\sigma (t)\sigma (0)]/\mathrm{\Delta }\gamma `$. A frequency-domain version of this technique was developed for polymers Booij (1966). Example data are shown in Fig. 2(a). Note that as the strain rate increases, the elastic character of the foam melts away. There are two signs of this. First, the progressively smaller intercept, $`G(0,\dot{\gamma })`$, means less elastic energy storage. Second, the progressively shorter decay time means more liquid-like dissipation. At high enough strain rates, where there is no transient storage and only dissipation, the foam behaves like an equilibrium liquid. Similar behavior has now been observed in simulations of a model foam Ono (2002).

To further quantify the unjamming behavior apparent in Fig. 2(a), we first deduce the relaxation spectrum, $`H(t,\dot{\gamma })dG(t,\dot{\gamma })/d\mathrm{ln}t`$ Ferry (1980); Macosko (1994). Results are shown in Fig. 2(b). For low and zero strain rates, there are two competing relaxation processes, reflected by a broad peak at very late times and long tail of short-time modes of nearly equal weight. The former reflects the coarsening process, and the latter the nonaffine bubble motion. As the shear rate increases, the coarsening peak gradually falls and another peak gradually rises over the nonaffine tail at short times. Presumably this is due to shear-induced rearrangements.

The salient features of the stress relaxation are shown in Fig. 3 vs strain rate. The first plot is of elastic storage, $`G(0,\dot{\gamma })`$, the value when the superimposed step-strain is achieved (below about 0.1s). This decreases very slowly, and is nearly constant, for strain rates less than about 0.05/s; for higher strain rates it abruptly vanishes. The second plot is of stress relaxation times. One such measure is $`t_e`$, when stress falls to $`1/e`$ of the initial value. Another measure is $`t_p`$, where $`H(t,\dot{\gamma })`$ reaches a global maximum. At zero and very low strain rates, these times are different since there are two competing relaxation mechanisms (evolution and shear). For strain rates higher than about $`3\times 10^4`$/s, the stress relaxation is essentially exponential and the two relaxation times are hence indistinguishable. In this regime, shear completely dominates the relaxation. It is puzzling that the relaxation time decreases with increasing strain rate as $`\dot{\gamma }^{1/2}`$, since on dimensional grounds one would have expected $`\dot{\gamma }^1`$.

Now we may compare the macroscopic rheology with the nature of the microscopic bubble dynamics. Previously we used DWS to measure the strain rate dependence of two microscopic time scales: $`\tau _o`$, the time between localized discrete rearrangements, and $`\tau _s`$, the time for adjacent scattering sites to convect apart by one wavelength of light Gopal and Durian (1999). The observed DWS data are reproduced by the dashed curves in Fig. 3(b). For very low strain rates, below about $`3\times 10^4`$/s, the rearrangements are discrete and the time between events, $`\tau _o`$, is not noticeably different from the quiescent value. The stress relaxation time is much much longer than $`\tau _o`$ \- indicating that the cooperative effect of many coarsening-induced rearrangements is required. For slightly higher rates, the time between events decreases, roughly as $`\dot{\gamma }^{1/2}`$, just like the stress relaxation time. Since the stress relaxation time is now shorter than $`\tau _o`$, not every site in the foam must rearrange in order to relax the overall stress. Shear-induced events begin to dominate over coarsening-induced events at strain rates near $`\dot{\gamma }_c=\gamma _y/\tau _{oq}=0.025`$/s, where $`\gamma _y=0.05`$ is the yield strain Gopal and Durian (1999). At still higher rates, events merge together and the flow becomes progressively more homogeneous and smooth. The crossover strain rate is $`\dot{\gamma }_m=\gamma _y/\tau _d=0.5`$/s, where $`\tau _d=0.1`$ s is the duration of rearrangements Gopal and Durian (1999). Many physical quantities display a change in character above and below $`\dot{\gamma }_m`$ Ono et al. (2002a). Indeed, just below $`\dot{\gamma }_m`$, the new rheological measures of elasticity vanish at precisely where we no longer can detect dicrete rearrangments (i.e. at the endpoint of the dashed curve for $`\tau _o`$). At higher strain rates, the bubble motion is dominated entirely by uniform shear, as characterized by the DWS time scale $`\tau _s`$. The onset of this correlated shearing motion of bubbles (i.e. the starting point of the $`\tau _s`$ curve) corresponds quite well to where $`G(0,\dot{\gamma })`$ begins to drop.

Finally, we note that the dramatic changes in the elastic character of the foam with strain rate are not reflected very strongly in the viscosity of the foam. As seen in Fig. 3(c), the viscosity decreases across the whole strain rate range, though not quite as fast as $`\dot{\gamma }^1`$ (which would have indicated constant stress). There is at most a slight change in exponent at $`\dot{\gamma }_m`$. Similar behavior is found in simulations Langer and Liu (2000). Altogether this emphasizes the importance of the superposition method to measure $`G(t,\dot{\gamma })`$, since it provides a clear dramatic signature of the unjamming transition. With this tool, and by comparison with DWS data, we have succeeded in connecting the nature of microscopic bubble dynamics with the resulting macroscopic rheological behavior. We have thereby shown that the unjamming of foam can be accomplished both by time and by application of shear. The unjammed liquid-like state is very similar to what would be achieved by raising the temperature for a thermal system. An important next step would be to deduce an effective temperature, recently shown in simulation of a sheared model foam to have many of the attributes of a true statistical mechanical temperature Ono et al. (2002b).

We thank A.J. Liu and P.T. Mather for helpful discussions. This work was supported by NASA Microgravity Fluid Physics grant NAG3-2481.


text:
# Young open clusters in the Carina region

## 1. Introduction

We have obtained deep UBVRI CCD photometry for Trumpler 14, Trumpler 16 and Collinder 232, the three young open clusters closest to $`\eta `$ Carinaæ . In this paper we briefly summarize our analysis methods and results. First of all we have tried to reconstruct the extinction pattern in the direction of the three clusters. By combining optical photometry with near infrared photometry and spectral classification for a sample of stars in the three clusters, we find by using a variety of methods that the total to selective absorption ratio varies from cluster to cluster, being $`R_V=3.48\pm 0.11,4.16\pm 0.07`$, and $`3.73\pm 0.01`$ for Trumpler 16, Trumpler 14 and Collinder 232, respectively. Then, by using the method devised by Romaniello et al. (2002), we derived individual star color excesses and built up reddening corrected Colour-Magnitude Diagrams (CMDs), from which we constrain clusters’s distances, finding that Trumpler 14, Trumpler 16 and Collinder 232 are located $`2.9\pm 0.3`$, $`4.4\pm 0.3`$, and $`2.5\pm 0.3`$ kpc from the Sun, respectively. The method devised by Romaniello et al. (2001) provides also stars luminosities and effective temperatures, that we used to construct the HR diagram (HRD) for the clusters, which, together with CMDs, have been used to derive clusters’ ages and age spreads. To this aim pre-Main Sequence (MS) isochrones have been adopted from Ventura et al. (1998). We find that Trumpler 16 is older than Trumpler 14, and that all the clusters have a prominent pre-MS population. As an example, in Fig. 1 we present our results for Collinder 232, that we conclude to be a real open cluster.

## References

Carraro, G., et al. (2002), A&A, submitted
Romaniello M., et al. (2002), AJ 123, 915
Ventura P., et al.1998, A&A 334, 953


text:
# “Vortex-melting front” in thin superconductors with pinning

## I Introduction

One of the most fascinating features of the vortex structure in layered superconductors is the so-called “melting” of the vortex ensemble. It is clearly established by now that this transition is a first-order phase transition above which the long-range correlation along the vortex lines is lost. The most thoroughly investigated material, in which the “vortex melting” occurs well below the second critical field, is the high-$`T_c`$ superconductor Bi2Sr2CaCu2Ox (BSCCO). There are attempts to add some other phase transition lines to the $`B`$-$`T`$ phase diagram of this material ($`B`$ is the local induction and $`T`$ the temperature). For example, the field of abrupt increase of the screening current with increasing $`B`$, the so-called second peak of the magnetization, was ascribed to a second order phase transition connected with the “melting” line at a tricritical point. However inspection of the available data suggests that the line of the second peak in the phase diagram is the continuation of the “melting” line into the low-temperature area where bulk currents do not completely relax during the experiment.

It is natural to assume that, in spite of the macroscopic irreversibility, the vortex phases are in equilibrium locally. Then the unified “melting line” including the second peak-line, which shifts as one curve under the introduction of additional disorder and/or change of anisotropy by oxygen doping or electron irradiation,  is the only first-order phase transition line. Recently the absence of the tricritical point at the border of observable pinning was elegantly demonstrated by unmasking the equilibrium magnetization step at the vortex melting point from the hysteretic magnetization using ac demagnetization.  This picture implies different pinning properties of these two vortex phases with correspondingly different current–voltage characteristics, which were extracted from local Hall probe measurements of the magnetic relaxation in the vicinity of the second magnetization peak at different temperatures and starting fields. The resulting jump of the screening current at the first-order phase transition is then observed as a second front of the penetrating magnetic flux.

Magneto-optical observations of this second magnetization front initially did not bring much new information and did not attract enough attention.  Such experiments were initially interpreted straightforwardly as the direct observation of the phase boundary which was already well established by local Hall probe measurements. However, after more careful consideration it became clear that the suggested “naive” interpretation is not compatible with the (macroscopic) electrodynamics of thin samples. Namely, in thin superconductors a sharp boundary between two different screening currents should produce the same logarithmic peak in the induction as at the sample edges. However, this moving peak is not observed in experiments. The correct theory indeed yields that a transition region with finite width and smooth current distribution should be expected, as it is known from the exact solutions for the penetration of perpendicular flux into thin samples, see e.g. Ref. . We propose here a model which predicts that this transition area in thin samples should inevitably contain a mixture of both phases.

A comprehensive description of the peculiarities of magnetic flux dynamics became particularly important when new observations of metastable vortex phases under fast remagnetization were reported.  In order to separate new dynamic effects from the (already complicated) geometry-dependent “regular” flux dynamics occurring when the vortex phases are in local equilibrium, we make a step backward in our investigations and describe the equilibrium case more carefully. The nature of flux flow at the phase transition is simulated qualitatively using an algorithm developed for the flux dynamics in thin conductors with nonlinear current-voltage law.

## II Magneto-optical observations

The penetration of a magnetic field into BSCCO single crystals was visualized by means of the magneto-optical imaging technique using a ferrimagnetic garnet indicator film with in-plane anisotropy.  In general, higher light intensity in the image means higher local magnetic induction normal to the sample surface. The images were recorded by a CCD camera during the continuous ramp of the perpendicular magnetic field with $`dB_a/dt=200`$ G s-1 and digitalized every 250 ms by a video-frame grabber card. The camera was working at the standard TV frame frequency of 25 Hz, so each image of the obtained movie is an average over about 40 ms. The increase of $`B_a`$ during the acquisition of a single image is comparable to the measurement noise. The relatively fast ramp of the field was used in order to observe the critical slope of the initial penetration, which otherwise relaxes. The intensity profiles measured across the digital images were calibrated in terms of field units using an additional recording under equivalent conditions but above $`T_c`$. A separate calibration curve was obtained from the additional “empty” images for each point on the flux profiles. In this way we reduce the effect of inhomogeneities in the illumination and in the sensitivity of the indicator across the image. For demonstration, an underdoped BSCCO crystal grown at the FOM-ALMOS center (the Netherlands) by the traveling-solvent floating zone technique under 130 mbar O2 partial pressure was specially selected for the uniformity of flux penetration and for the low field $`B_{\mathrm{sp}}180`$ G at which the second magnetization peak occurs. The crystal was cut into a rectangle with dimensions $`1450\times 700\times 43\mu `$m3.

The initial penetration into the zero-field-cooled sample resembles the numerous magneto-optical observations of the Bean-like critical states in type-II superconductors. The crystal edges are well observed as a sharp change from the very bright field outside the sample and the grey interior. The flux penetrating from the crystal edges is visible as a grey belt (higher brightness represents higher

local value of the perpendicular magnetic induction) which squeezes the central (dark grey) Meissner area in which $`B`$ is zero (Fig. 1 a). The gradual decrease of the brightness of the penetrated perimeter corresponds to the well-known Bean gradient of the penetrating flux. When the field at the edges reaches $`B_{\mathrm{sp}}`$, a second sharp flux front starts to penetrate from the crystal edges (Fig. 1 b).  When the second flux front has penetrated deep into the sample, the resulting image is becoming nearly identical to the initial penetration (compare Fig. 1 a and b). However, one should remember that the induction in the central dark area in Fig. 1 c approximately corresponds to that in the light grey field outside the sample in Fig. 1 a, where $`B=0`$ in the center.

The whole course of flux penetration during field ramping from zero to 500 G is visualized in Fig. 2 by $`B(x)`$ profiles measured across the sample width. After reaching a well-determined field level the increase of the induction $`B`$ pauses and $`B`$ is stabilized until a new flux front approaches which increases $`B`$ further. The plateau of constant induction is associated with the phase transition at $`B_{\mathrm{sp}}`$. The slight inclination of this plateau, which is marked by two dashed lines, is due to inhomogeneity of the crystal. To our knowledge, such inhomogeneity is found in all BSCCO crystals and is inevitable (see e.g. Refs. ). During the subsequent decrease of the field a similar feature appears at the same line $`B_{\mathrm{sp}}`$ without notable hysteresis in the transition field. The second flux front looks very similar to what the ordinary Bean-like flux profile would have looked like if some flux relaxation had occurred during a temporal halt of the field ramp. Another close analogy exists with respect to flux penetration into a superconductor after field cooling. That analogy, in particular, shows that the initial interpretation of the second front as a phase boundary separating regions with very different critical currents, is not compatible with the perpendicular geometry of the experiment. Such a sharp border between two currents at a definite field is possible only in the longitudinal geometry, where it occurs at the junction of two different field gradients. In perpendicular geometry, however, any flux front penetrating into a thin superconductor is accompanied by strong sheet currents inside the screened area.

## III Model calculations and discussion

In order to understand the experimental observations, let us first consider the simplest possible situation of a superconducting strip occupying $`|x|<a`$, $`|z|<d/2`$ ($`da`$) in which the vortex matter undergoes a first order transition at $`B_z=B_{\mathrm{sp}}`$. In the vortex solid below $`B_{\mathrm{sp}}`$, the critical current is zero, and in the vortex liquid above $`B_{\mathrm{sp}}`$, it has some non-zero value $`J_c`$. In other words, the critical current discontinuously jumps from zero to $`J_c`$ as the vortex matter melts. Furthermore, let us neglect surface and geometrical barriers. The flux penetration into this strip is illustrated in Fig. 3. As long as the applied field $`B_az`$ does not exceed $`B_{\mathrm{sp}}`$, the strip is perfectly transparent to the field (flat profiles 1 and 2 in Fig. 3). Once the applied field exceeds $`B_{\mathrm{sp}}`$, the internal induction is stabilized at the value $`B_{\mathrm{sp}}`$, and further flux penetration is screened by a non-zero $`J(x)`$. The distribution of this $`J(x)`$ and of the internal induction are described by the well–known solution of Ref. :

$`\begin{array}{ccc}J(x)& =& \{\begin{array}{cc}\left(2J_c/\pi \right)\mathrm{arctan}\frac{cx}{\sqrt{b^2x^2}},\hfill & |x|<b\hfill \\ J_cx/|x|,\hfill & b<|x|<a\hfill \end{array}\\ & & \\ B(x)& =& \{\begin{array}{cc}B_{\mathrm{sp}},\hfill & |x|<b\hfill \\ B_{\mathrm{sp}}+B_c\mathrm{arctanh}\frac{\sqrt{x^2b^2}}{c|x|},\hfill & b<|x|a\hfill \end{array}\end{array}`$

where $`B_c=\pi \mu _0J_c/d`$ is the characteristic field, $`\pm b=a/\mathrm{cosh}[(B_aB_{\mathrm{sp}})/B_c]`$ is the position of the flux front, and $`c=\sqrt{a^2b^2}`$. For Fig. 3, we chose $`B_c=\frac{1}{2}B_{\mathrm{sp}}`$. The situation is illustrated by the profiles 3 and 4. In the region in which $`B>B_{\mathrm{sp}}`$, the current $`J=J_c`$, while in the region with $`B=B_{\mathrm{sp}}`$, the strip geometry imposes the smooth decrease of $`J`$ to zero in the interval between the “second” flux front and the sample center. While it is evident that the former region is occupied by the high field vortex liquid, the latter region must also consist of vortex liquid, since the vortex solid (with zero critical current) cannot sustain the non-zero $`J(x)`$. Therefore, the observed flux front does not correspond to the front of phase transformation.

The above consideration, based on the simple Bean model, does not account for a non-zero current in the solid phase, nor for the possibility of flux creep. We include these effects using the algorithm previously developed for studying the flux penetration into thin superconductors with a nonlinear current–voltage law $`E(J,J_c)=\rho _0|J/J_c|^\sigma J`$ where $`J`$ is the sheet current (the current density integrated over the strip thickness).  We model the sharp jump of the current-voltage law from $`E^{}(J)=E(J,J_c^{})`$ in the low-field low-pinning phase ($`J_c=J_c^{}`$) to $`E^+(J)=E(J,J_c^+)`$ in the high-field high-pinning phase ($`J_c=J_c^+`$) by introducing a step in $`J_c(B)`$ at the transition field $`B_{\mathrm{sp}}`$:

$$J_c(B)=J_{c0}+\frac{\mathrm{\Delta }J_c}{2}\mathrm{tanh}\frac{BB_{\mathrm{sp}}}{\mathrm{\Delta }B},$$
(1)
with jump width $`\mathrm{\Delta }BB_{\mathrm{sp}}`$ and jump height $`\mathrm{\Delta }J_c=J_c^+J_c^{}`$. Other models of steps in the current–voltage characteristics, e.g. a jump in the creep exponent $`\sigma (B)`$, or different values of the parameters, give qualitatively similar results.

The following graphs (Fig. 4 a–c) show profiles (across the strip) of the magnetic field $`B(x)`$, screening current $`J(x)`$, and the electric field $`E(x)`$ under a constant ramp rate of the external field $`dB_a/dt=1`$ for $`\sigma =20`$, $`B_{\mathrm{sp}}=1`$, $`\mathrm{\Delta }J_c=0.5`$ and $`\mathrm{\Delta }B=0.04`$, in dimensionless units where $`J_c^+=1`$ (critical sheet current), $`a=1`$ (half width of the strip), $`B=H`$, $`dE/dx=dB/dt`$ (left-hand coordinate system), and $`\rho _0=1`$. Our computed profiles of the induction $`B(x)`$ (Fig. 4a) well resemble the experimental profiles (Fig. 2). The initial penetration at $`B<B_{\mathrm{sp}}`$ is the same as was calculated analytically  and numerically . As usual, the corresponding $`J(x)`$ in the flux penetrated area saturates at the “spatial projection” $`J_{sat}^{}`$ of the low field current–voltage characteristic. The reason for this is very simple: after the complete penetration of the field, $`B`$ increases everywhere with the same constant rate. From the Maxwell law $`dE(x)/dx=dB(x)/dt`$ it follows that $`E^{}(J_{\mathrm{sat}}^{}(x))=xdB_a/dtx`$ everywhere. Thus $`J_{\mathrm{sat}}^{}(x)=`$ $`J^{}(xdB_a/dt)`$, where $`J^{}(E)`$ is the inverse function of

$`E^{}(J)`$ defined below $`B_{\mathrm{sp}}`$. The constant critical current of the classical Bean model corresponds to the nearly constant flat part of $`J_{\mathrm{sat}}^{}(x)J_{\mathrm{sat}}^{}`$. Unlike the Bean model, in our more realistic model the abrupt change of direction of the critical current at $`x=0`$, the so-called discontinuity line (d-line), has a finite thickness. The transition from $`+0.8J_{\mathrm{sat}}^{}`$ to $`+0.8J_{\mathrm{sat}}^{}`$ extends over $`\mathrm{\Delta }x_d2E^{}(0.8J_c^{})/(dB_a/dt)`$. For the particular $`E(J)`$ used in our simulations (Fig. 4) $`\mathrm{\Delta }x_d0.02a`$.

When $`B`$ increases further and approaches the value $`B_{\mathrm{sp}}`$, the field gradient decreases, and further flux penetration occurs via the penetration of a new flux front moving on top of the line $`BB_{\mathrm{sp}}`$. It is important to note that, as in the simplified case outlined above, the screening current $`J(x)`$ does not jump sharply at the flux front but changes smoothly from $`J_{\mathrm{sat}}^{}(x)`$ to the new saturation $`J_{\mathrm{sat}}^+(x)J_{\mathrm{sat}}^+`$ before the second front has arrived (Fig. 4b). Here again $`E^+(J_{\mathrm{sat}}^+(x))=xdB/dtx`$ yields the “spatial projection” $`J_{\mathrm{sat}}^+(x)`$ of the high-field current-voltage law. This behavior is valid for any pair of functions $`E^{}(J)`$ and $`E^+(J)`$, which play the role of saturation profiles of the sheet current $`J(x)`$ after the complete penetration of the field, below and above $`B_{\mathrm{sp}}`$, correspondingly. Preceding the second front there is a region where the flux density is nearly constant in space and time, $`B(x,t)B_{\mathrm{sp}}`$ (Fig. 4 a); in this region the slope of the electric field is reduced since $`dE/dx=dB/dt`$ is small. In order to reach its final saturation shape $`E(x)=xdB_a/dt`$, $`E(x)`$ exhibits a sharp upward bend at the second flux front, which is clearly seen in Fig. 4c. We should mention that a non-zero $`\mathrm{\Delta }B`$ is included here only for the stability of the numerical algorithm; its value only affects the rounding of sharp features in the flux profile for $`BB_{\mathrm{sp}}`$, but not, as we checked, the wide smooth transition of the current. The simplified solution outlined above shows that the smooth transition in $`J`$ exists for $`\mathrm{\Delta }B=0`$ as well.

Within the framework of our assumption that only two phases with fixed current–voltage laws $`E^{}(J)`$ and $`E^+(J)`$ can exist in the superconductor, the above outlined smooth transition of the resulting $`J(x)`$ between the corresponding saturation profiles can be explained only by a phase mixture. One can imagine various configurations of the phase domains corresponding to different conditions, varying from the equality of $`E`$ in the two phases to the homogeneity of $`J`$. In general the dynamics of the corresponding phase mixture should be solved, but this complication is outside the scope of our paper. Here, we present the simplest and most natural model in which the intermediate region consists of narrow lamellar domains of alternating phases oriented parallel to the flux flow, i.e.perpendicular to the current, across the width of the strip. This model corresponds to homogeneous $`J`$ across the phase boundaries and a modulated $`E`$. In the mixed region, we then have $`J_{\mathrm{liquid}}=`$ $`J_{\mathrm{solid}}=`$ $`J(x)`$ and $`E(x)=nE^+(J(x))+(1n)E^{}(J(x))`$. Here, $`n=n(x)`$ is the relative concentration of one phase; $`n=0`$ corresponds to the “vortex solid” and $`n=1`$ to the “vortex liquid”. From the calculated $`J(x)`$ and $`E(x)`$ we obtain $`n(x)`$ in the form

$$n(x)=\frac{E(x)E^{}(J(x))}{E^+(J(x))E^{}(J(x))}$$
(2)
(see Fig. 4d). In case of the simplest domain structure of the phase mixture, equal parallel lamellae along the flux flow, the relative width of the “liquid” domains is equal to $`n(x)`$. An example of such domains at $`B_a=1.4`$ is shown in Fig. 5, which is obtained by flipping the corresponding profile $`n(x)`$ (outlined by the bold line in Fig. 4c) several times and cutting it at the right at the position of the second flux front, and at the left, e.g. at a value $`n(x)=0.01`$ ($`0.01`$ is an arbitrary choice, the final result depends on it logarithmically).

Interestingly, “melting” (i.e. $`n(x)>0`$) starts well before the second flux front has arrived, and the second flux front corresponds not to the midpoint $`n=\frac{1}{2}`$, but to the point where $`n(x)`$ has already well saturated to $`n=1`$. At this point the “melting” transition is complete. The smooth transition of the current is characterized by the fast flux flow through the weaker “solid” lamella (light grey), the width of which plays the role of a “flux valve”, while the vortices in the “liquid” lamella (white) are nearly immobile until the second flux front arrives. At $`xx_{\mathrm{sp}}`$, the average electric field equals

$$E(x)E_{\mathrm{sp}}=_0^{x_{\mathrm{sp}}}\frac{dB(x,t)}{dt}𝑑x\mathrm{const}$$

because $`B(x,t)`$ is fixed to a value $`B_{\mathrm{sp}}`$. The vortex flow through the region of the phase mixture is concentrated in the “solid” phase. The local electric field (vortex velocity) in the solid is very high, $`E_{\mathrm{solid}}=E^{}(J(x))=E_{\mathrm{sp}}`$, because the electric field in the “liquid” phase $`E_{\mathrm{liquid}}=E^+(J(x))`$ is nearly zero. Consequently, the width of the “solid” valve adapts itself so that the necessary average flux flow is conserved, i.e. $`1n_{\mathrm{sp}}E_{\mathrm{sp}}/E_{\mathrm{solid}}`$.

Particularly impressive is the moment when the second front approaches and the high second critical current $`J_{\mathrm{sat}}^+`$ passes through the “weak” “solid” phase. Here all necessary flux flow squeezes into the very narrow “solid neck” with a huge speed corresponding to $`E^{}(J_{\mathrm{sat}}^+)`$. The order of the introduced values $`E_{\mathrm{liquid}}`$ and $`E_{\mathrm{solid}}`$ can be read from the “spatial projections” of $`E^+(J)`$ and $`E^{}(J)`$ in Fig. 4b (see the corresponding lower and upper asymptotics $`J_{\mathrm{sat}}^+(x)`$ and $`J_{\mathrm{sat}}^{}(x)`$). For this one has to take a horizontal projection of, for example, the cross-point of the bold line $`J(x)`$ (at $`B_a=1.4`$) and the vertical dashed line, indicating the second flux front, onto $`J_{\mathrm{sat}}^+(x)`$ and $`J_{\mathrm{sat}}^{}(x)`$. Then the first one will give $`x^+=E_{\mathrm{liquid}}0.05`$ and the second one will be well outside the graph at the right $`x^{}=E_{\mathrm{solid}}10^3`$. The relative width of the “solid neck” is here $`1n_{\mathrm{sp}}2\times 10^4`$. In reality, even for such domains as thick as the sample the “neck” width given by our qualitative theory is well below the vortex spacing of 0.3 $`\mu `$m at $`B_{\mathrm{sp}}`$. This means that probably only a single line of “solid” vortices will move before the second front arrives.

The proposed dynamical intermediate structure of the coexisting phases should not be confused with the ordinary intermediate phase that exists at any first order magnetic phase transition in the narrow field range $`\mathrm{\Delta }B=B^+B^{}`$ over which the jump occurs. The relative size of the equilibrium domains $`n_{\mathrm{eq}}`$ in that well known static case is fixed by the condition that the average local field equals the $`B(x)`$ which is determined by the entire current distribution, $`n_{\mathrm{eq}}B^++(1n_{\mathrm{eq}})B^{}=B(x)`$.

As opposed to this, in our case the relative phase ratio $`n(x)`$ is determined by the different dynamics of the flux flow in the two different phases and by the global electrodynamics of the system, see the profiles in Fig. 4d and Eq. (2). In the ordinary static case, the region of the coexisting phases may be very narrow because $`\mathrm{\Delta }B`$ is much smaller than the typical variation of $`B(x)`$ over the sample; its width equals $`(dB/dx)^1\mathrm{\Delta }B`$ at the phase transition. In our case the width of the transition region containing the phase mixture is a sizable fraction of the sample size, while our dynamic $`B(x)`$ varies even more due to the presence of high bulk currents.

The question arises whether the intermediate structure of dynamically coexisting solid and liquid vortex phases also occurs in thick samples (thicker than the penetration depth $`\lambda `$). In the case of the simplified model with zero current for $`B<B_{\mathrm{sp}}`$ and a critical current $`J_c`$ for $`B>B_{\mathrm{sp}}`$, we refer to Ref. , in which it is shown that the shape of the flux front in the interior of a thick superconductor is very well described by the shape of the current profile in an infinitely thin superconductor. Using Fig. 3, we can thus immediately draw the current front inside a thick sample. The region between this front and the sample surface carries a current $`J_c`$ and must therefore be in the liquid state. From the Bean model only, there is no objection to the inner core remaining in the solid state with $`J_c=0`$. In that case, one would also have a coexistence of solid and liquid phases, but now across the sample thickness. Such a scenario would depend on the stability of the curved phase boundary — any phase mixture in the inner core up to its complete transformation to the solid phase also does not contradict to electrodynamics.

Recently, a real-time differential treatment of magneto-optical images allowed to visualize the propagating vortex “melting” front through a BSCCO crystal as a contour of the related equilibrium magnetization jump.  The melting point appeared to be very sensitive to tiny inhomogeneities of the sample defect disorder and a very irregular “vortex melting” with no sign of equilibrium domains is observed even in the best crystals. Our lamella should be more robust because they are determined by the critical currents, which are much higher than the tiny equilibrium current drop at the “vortex melting” point. But unfortunately, the elaborate procedure of such observations cannot be applied to our case of irreversible magnetization. Very recently, after preparation of this paper, a mixture of phases with different $`E(J)`$ has been observed in the peak effect region in NbSe2 using a scanning Hall probe magnetometer and local ac excitation.  We have to mention that the phenomenon observed there, is static and the suggested model of the “vortex capacitor” differs from our dynamic model. Furthermore, we deal with an equilibrium transition, contrary to the metastable state considered in Refs.  and . The new observation technique, nevertheless, may be applicable to our case; its magneto-optical imaging counterpart can be realized as well. Note however that in our case the picture is evolving rapidly (on the time scale of 10-2 s) and the irreversible currents which would screen the small ac excitation are much higher.

## IV Conclusions

We have obtained experimental and numerical flux profiles that render the major features of flux flow at “vortex-lattice melting” in the presence of pinning and magnetic hysteresis. The flux profiles show a second front corresponding to the well-known second peak in the magnetization curve. In thin superconducting plates, “melting” occurs well before the arrival of the second flux front. At that moment, the phase transition is complete and the low-field phase (“vortex solid”) with lower $`J_c`$ completely disappears there. The appearance of the high field phase (“vortex liquid”) does not produce any sharp feature on the observable field profiles. Rather, there is a large region preceding the second flux front in which domains of “solid” and “liquid” coexist dynamically: the relative width of these domains is controlled by the electrical field induced by the vortices moving under the force exerted by the current, which flows continuously across the domains; the current in its turn is controlled by the overall electrodynamics of the system. This intricate feature should be taken into account in a correct interpretation of experiments.

Evidently this effect occurs for any flux dynamics in the presence of a first order phase transition, irrespective of whether it is driven by ramping the applied field as in the presented experiment, or it happens during magnetic relaxation as in Ref. . Moreover, when the process is slow and the current in the “weak” “solid vortex-lattice” phase is allowed to relax to zero, the “melting” proceeds nearly immediately everywhere because the flux flow required for this is very low. In this case the remaining “solid” lamella should be very narrow and the second front moves on top of the nearly homogeneous “liquid” without any relation to the phase transition itself.

The authors are grateful to Ming Li and P.H. Kes of the FOM-ALMOS center (the Netherlands) for the samples. M.V.I. ***Contact author:
Mikhail Indenbom, e-mail: indenbom@issp.ac.ru acknowledges financial support by the Max-Planck-Gesellschaft.