text:
# Dynamic cores in hydrostatic disguise

## 1 Introduction

Stars and planets form in dense cores within molecular clouds. In the case of low-mass star forming regions, it has traditionally been thought that these cores are quasi-static equilibrium configurations supported against gravitational collapse by a combination of magnetic and thermal pressures (see, e.g. Shu, Adams & Lizano 1987). When considering purely thermal support, stable hydrostatic solutions of the self-gravitating fluid equations with finite central densities and sizes exist, provided the sphere is confined by an external pressure $`P_{\mathrm{ext}}`$ acting on the surface. These solutions are known as Bonnor-Ebert (BE) spheres (Ebert, 1955; Bonnor, 1956). Recently, Alves, Lada & Lada(2001) have determined with small error bars the column density profile of B68, a core embedded in an HII region. They have also successfully fitted a BE profile. Similar studies (Johnstone et al., 2000; Harvey et al., 2001; Evans et al., 2001) have attempted to fit BE profiles to the column density profiles of selected molecular cloud cores.

However, the picture of isothermal cores in hydrostatic equilibrium may be in conflict with the fact that molecular clouds are turbulent. It appears difficult that quasistatic equilibrium structures may appear and survive in isothermal, supersonic, highly compressible turbulent flows, where density fluctuations are in general transient, and have further substructure in a self-similar hierarchy (Scalo, 1985; Houlahan & Scalo, 1990; Falgarone, Phillips & Walker, 1991; Houlahan & Scalo, 1992; Vázquez-Semadeni, 1994; Vázquez-Semadeni et al., 2000; Williams, Blitz & McKee, 2000; Chappell & Scalo, 2001; Klessen, 2001; Mac Low & Klessen, 2003), which may possibly end at small enough scales that the turbulent velocity fluctuations are no longer supersonic, and cannot produce further turbulent fragmentation111Note that by small enough scale we just mean a typical average scale at which the velocity dispersion equals the sound speed, but around which there may be a large scatter in the sizes of individual clumps satisfying this condition. (Vázquez-Semadeni, Ballesteros-Paredes & Klessen 2003a). Nevertheless, in the turbulent case, there is no reason why these structureless density fluctuations should be hydrostatic, as they have formed from larger-scale turbulent compressions within a uniform-temperature medium. It is more likely that they should either re-expand or proceed to gravitational fragmentation and collapse, since the probability of them reaching precise balance between self-gravity and their internal pressure is vanishingly small (Tohline, Bodenheimer & Christodoulou 1987; Taylor, Morata & Williams 1996; Ballesteros-Paredes, Vázquez-Semadeni & Scalo 1999; Vázquez-Semadeni, Shadmehri & Ballesteros-Paredes 2003b). Note that the condition for constructing a BE sphere, namely the availability of a hotter, more tenuous confining medium is in general not realized between molecular clouds and their embedded cores, as they are both at roughly the same temperature.

The transient character of most turbulent fluctuations and the eventual induction of collapse on one of them in a gravitationally unbound medium was shown by Vázquez-Semadeni, Passot & Pouquet (1996). Moreover, simulations of turbulent, globally gravitationally bound clouds consistently show the collapse of local peaks, but never the formation of hydrostatic cores, either in the purely hydrodynamic case (Léorat, Passot & Pouquet 1990; Klessen, Heitsch & Mac Low 2000) or in the magneto-hydrodynamic (MHD) case (Heitsch, Mac Low & Klessen, 2001). One exception is the case of MHD simulations in boxes with subcritical mass-to-magnetic flux ratios (Ostriker, Stone & Gammie 1999), which evolve towards flattened sheets that cannot collapse in the direction perpendicular to the field. However, being performed in closed boxes, these simulations cannot represent the fact that more mass can continue to accrete along field lines until a core becomes supercritical (Hartmann, Ballesteros-Paredes & Bergin 2001). There also exist both observational evidence (Crutcher, 1999; Andre, Ward-Thompson & Barsony, 2000; Bourke et al., 2001; Hartmann et al., 2001) and theoretical arguments (Nakano 1998; Hartmann et al. 2001; as summarized by Mac Low & Klessen 2003) showing that the overwhelming majority of cores must be supercritical. Moreover, the long lifetimes that quasi-static cores would have are difficult to reconcile with observational statistics of cloud cores (Taylor et al., 1996; Lee & Myers, 1999; Visser, Richer & Chandler, 2002), and with the suggestion of short molecular cloud formation time scales ($``$ a few Myr), based on the observed lack of post-T-Tauri stars in Taurus (Herbig, 1978; Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999; Hartmann 2002).

It is important to mention that the preceding discussion is not in contradiction with the fact that stars are objects in hydrostatic equilibrium within a turbulent medium, since they do not correspond to quasi-isothermal flows. In stars, energy is trapped since the opacity has increased, and the cooling time is about 1010 times the free-fall time. This is also consistent with the fact that the formation of a stable hydrostatic structure requires a polytropic exponent $`\gamma >4/3`$, where $`P\rho ^\gamma `$ (Chandrasekhar, 1939). On the other hand, molecular clouds and their cores are not able to trap energy effectively because they are approximately isothermal (Goldsmith & Langer, 1978; Goldsmith, 1988; Scalo et al., 1998; Ballesteros-Paredes et al., 1999b; Spaans & Silk, 2000). Thus, they are not expected to be able to reach a hydrostatic configuration if they are formed by a turbulent compression (Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999; Vázquez-Semadeni, et. al 2003b). This conclusion may hold even when magnetic fields are considered because, as mentioned above, the existence of subcritical cores is being questioned, and supercritical configurations are qualitatively equivalent to non-magnetic ones. Finally, even is magnetically subcritical cores do exist in clouds, their flattened configurations are very different from those of the BE sphere, and no reason exists for them to be well fitted by a BE profile.

In view of the above, it is clear that the recent attempts to fit hydrostatic BE profiles to observed molecular cloud cores appear to be in serious contradiction with the theoretical results suggesting that non-magnetic hydrostatic structures are unlikely to form in a turbulent medium and that the only hydrostatic magnetic configurations ought to be flattened. In this paper, we show that the apparent discrepancy can be resolved because fitting a BE-type profile to an observed cloud core is not an unambiguous test of it being in hydrostatic equilibrium. To show this, we take some of the cores (which are not in hydrostatic equilibrium) in the numerical simulations of Klessen, Burkert & Bate (1998), Klessen & Burkert (2000), Klessen & Burkert (2001), and Klessen et al. (2000), and apply to them a fitting procedure similar to that used by Alves et al. (2001), showing that reasonable BE column-density profile fits can be made on these non hydrostatic structures as well.

The plan of the paper is as follows. In §2 we summarize the main characteristics of the simulations analyzed. In §3 we explain the numerical technique in fitting BE column density profiles, and in §4 we present the results of the fits. In §5 we discuss the results, and present the main conclusions.

## 2 Numerical Simulations

A prerequisite for adequately describing the inner structure of cores in models of turbulent molecular clouds is to use a numerical technique that is able to resolve high density contrasts, at arbitrary locations within the cloud. The method of choice is smoothed particle hydrodynamics (SPH), which is a particle-based Lagrangian scheme to solve the equations of hydrodynamics. Flow properties are obtained by averaging over an appropriate subset of SPH particles, and high resolution is achieved where needed by increasing the particle concentration and reducing the averaging volume (Benz, 1990; Monaghan, 1992).

For the current investigation we re-analyze numerical models of turbulent molecular cloud evolution first presented by Klessen et al. (1998); Klessen & Burkert (2000, 2001) and Klessen et al. (2000) using 205,379 particles. These calculations solve the equations for a self-gravitating isothermal ideal gas in cubic volumes with periodic boundary conditions, and are imagined to be located within larger, roughly self-similar clouds. For the density range of interest ($`10^2\mathrm{cm}^3n(\mathrm{H}_2)10^7\mathrm{cm}^3`$) interstellar molecular gas is seen at temperatures close to $`10`$K and the isothermal equation of state is approximately valid (Goldsmith & Langer, 1978; Goldsmith, 1988; Scalo et al., 1998; Ballesteros-Paredes et al., 1999b; Spaans & Silk, 2000). Under these conditions the dynamical behavior of the gas is scale-free and depends only on the ratio of internal and turbulent kinetic energy to gravitational energy. The numerical simulations are thus performed in dimensionless units and need to be rescaled to obtain physical quantities. In the simulations, once the density in a region exceeds a density contrast of $`10^4`$, a sink particle is created, with the same mass of the region where it formed, but with a fixed radius of the order of the Jeans length at the threshold density (Bate, Bonnell & Price, 1995). The internal structure of sink particles is not resolved.

To facilitate a direct comparison with the observations we decompose the gas distribution in each numerical model into individual clumps using a 3-dimensional clump-finding algorithm (see Appendix A in Klessen & Burkert 2000) and select the first 21 gas clumps without sink particles (as their analysis would be meaningless) determined by the clump-finding algorithm.

Once we have the individual clumps, we define the centers of control subregions (with volume $`(1/10)^3`$ of the full simulation cube) as the location of the density peak. The SPH density distribution within a control subregion is assigned onto a $`128^3`$ grid, and integration along the principal axes produces three column-density maps – one for each direction of projection. Note that the location of the density peak needs not correspond to the maximum of the projected column density.

To compare our maps with observed stellar extinction maps of dark globules, we adopt a physical scaling such that, for a typical protostellar core in the simulation, the resulting column density map and pixel sizes roughly correspond to the observations of the Bok globule B68 by Alves et al. (2001). Our maps then have a physical size of $`0.18`$pc$`\times 0.18`$pc and column densities $`N(\mathrm{H}_2)`$ ranging from about $`10^{20}`$cm-2 for the background value to a few times $`10^{22}`$cm-2 at the core centers. For details on the scaling procedure the reader is referred to Klessen & Burkert (2000) and Klessen et al. (2000).

We analyze three different numerical models. These constitute the extreme ends of the possible range of molecular cloud dynamics. In the first two models, the molecular gas is subject to strong supersonic turbulence. In the first (hereafter LSD), turbulence is driven on large scales, with wavelength $`\lambda 1/2`$ of the computational box. In the second model (denoted SSD), energy is injected on small scales, with $`\lambda 1/8`$ of the computational box. With a box size of 1.54 pc per side, these correspond to driving scales of $`0.77`$pc and $`0.19`$pc, respectively. For each model we consider two different evolutionary stages. First, a stage ($`=t_0`$) of fully developed supersonic turbulence without self-gravity, and second, one snapshot after self-gravity has been switched on and gravitational contraction has led to the formation of collapsed cores accumulating roughly 5% of the total mass ($`=t_1`$). These two models are complemented by one simulation without driving, and where the gas is allowed to collapse freely from a field of random Gaussian density fluctuations (hereafter GC) under the influence of self-gravity. Again, two times are considered, $`t=t_0`$, now corresponding to the evolutionary phase just before the first collapsed core occurs, and $`t=t_1`$, when $`5`$% of the mass is in collapsed objects. For a more detailed discussion of the dynamical evolution the models and the implications for star formation, see Klessen (2001). The properties relevant for the current investigation are summarized in Table 1.

## 3 Analysis of the Simulations

### 3.1 Bonnor-Ebert Models

In the classical analysis of Ebert (1955) and Bonnor (1956), the equation that describes a self-gravitating isothermal sphere in hydrostatic equilibrium is the modified Lane-Emden equation:

$$\frac{d}{d\xi }\left(\xi ^2\frac{d\mathrm{\Psi }}{d\xi }\right)=\xi ^2e^\mathrm{\Psi }$$
(1)
where $`\xi =(r/c_s)\sqrt{4\pi G\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}}`$ is a non-dimensional radial variable, $`r`$ is the radial coordinate, $`c_s=\sqrt{kT_{\mathrm{core}}/m}`$ is the thermal velocity at temperature $`T_{\mathrm{core}}`$, $`k`$ is the Boltzmann constant, $`G`$ is the gravitational constant, $`\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}=m_pn_{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}`$ is the central mass density, $`m_p`$ is the mean mass per particle, $`n_{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}`$ is the central number density of particles, and $`\mathrm{\Psi }=\mathrm{ln}(\rho /\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}})`$. The solutions of eq. (1) are required to be finite at the center of the core, which implies that they must satisfy the boundary conditions $`\mathrm{\Psi }(0)=0`$ and $`d\mathrm{\Psi }/d\xi (0)=0.`$

Equation (1) can be solved by the change of variables $`y_0=\xi ^2\eta `$, $`y_1=\mathrm{\Psi }`$ (with $`\eta =d\mathrm{\Psi }/d\xi `$), which splits up the equation into two first-order equations:

$`{\displaystyle \frac{dy_1}{d\xi }}`$ $`=`$ $`{\displaystyle \frac{y_0}{\xi ^2}}`$ (2)
$`{\displaystyle \frac{dy_0}{d\xi }}`$ $`=`$ $`\xi ^2\mathrm{exp}(y_1),`$ (3)
with the initial conditions $`y_0(0)=0`$, and $`y_1(0)=0`$. This system has a unique solution out to the nondimensional parameter

$$\xi _{\mathrm{max}}=\frac{R_{\mathrm{core}}}{c_s}\sqrt{4\pi G\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}},$$
(4)
where $`\xi _{\mathrm{max}}`$ is the value of $`\xi `$ at the outer boundary $`R_{\mathrm{core}}`$, beyond which the density is negligible, but the temperature is assumed to be sufficiently high as to maintain pressure equilibrium at the boundary.

Here we stress two points. First, a particular solution is defined by the radius, temperature and central density of the core, $`R_{\mathrm{core}}`$, $`T_{\mathrm{core}}`$, and $`\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}`$; and second, if the solution has $`\xi _{\mathrm{max}}>`$ 6.5 (i.e., if the density contrast between the center and the edge of the BE sphere is larger than 14.3), the configuration is in an unstable equilibrium. We will get back to these points in next section.

### 3.2 Numerical Technique

#### 3.2.1 Construction of column density maps and radial profiles

In order to mimic the observational procedure, we assume that the center of the BE sphere is located at the position of the column density maximum, $`(x_0,y_0)`$. Then, we slice the column density map into 36 angular sectors of width 10, centered on $`(x_0,y_0)`$. For each sector we construct the radial column density profile. If in a particular sector more than one measurement falls into the same radial bin, we take the arithmetic average. This procedure gives a column density profile in which the values of the abscissas (the radii) are unevenly spaced. We thus redefine the profile as an interpolation of the latter at homogeneously distributed radii. Finally, we define the angle-averaged radial column density as the average profile over the 36 slices. We also define “dispersion” profiles as the mean profile $`\pm `$ the standard deviation of the 36 profiles. We generically refer to both of them as the “dispersion profile”.

Our cubic subregions typically contain several thousand particles. For the inner parts (i.e. the central regions of the considered clumps), the subregion resolution of $`128^3`$ cells roughly reflects the effective spatial resolution of the SPH method in these high-density regions. Towards the outer parts, where the density is lower, the smoothing volume of an SPH particle may exceed the cell size which we adopt, so that in the outer low-density wings of the cores we do not fully resolve the small-scale structure. However, this does not influence the BE fitting procedure significantly, because we use logarithmic bins in the radial direction for constructing the radial profiles, implying that we have larger stepsizes in the low density regions corresponding to the lower spatial resolution there. We expect that the additional smoothing in the wings leads to a slight overestimate of the “true” density (i.e. if we had infinite mass and spatial resolution), because we have a positive gradient inwards. However, there being one or more decades in column density below the peak, we expect the effect on the overall result of the BE fit to be negligible.

#### 3.2.2 Fitting procedure

We calculate the solutions of the modified Lane-Emden equation by solving equations (2) and (3) using an ordinary fourth order Runge-Kutta method. Similarly to the criterion used by Alves et al. (2001), who fitted the observational data only out to the radius where the column density profile started to deviate significantly from a BE-like profile (roughly 100 arcsec in their fig. 2; J. Alves, private communication), in this paper we choose the value of outer radius of the core ($`R_{\mathrm{core}}`$) by inspection of the column density maps, selecting either the radius at which noticeable kinks start to appear in the average column density profile, or else that at which the profile begins to clearly deviate from a BE one. Subsequently, we calculate a set of BE spheres by varying both the central density and the temperature (see below), and we define the BE profile that “best matches” the core as the one that minimizes the function:

$$\chi ^2=\frac{1}{m}\underset{i}{}\left[\mathrm{log}(N_{i,\mathrm{BE}})\mathrm{log}(N_{i,\mathrm{core}})\right]^2,$$
(5)
where $`N_{\mathrm{BE},i}`$ and $`N_{\mathrm{core},i}`$ are, respectively, the column density of the fitted Bonnor-Ebert sphere, and of the actual column density of the core at the $`i`$-th radius, and $`m`$ is the number of points used to construct the radial profile. The sum extends over all points at which $`N_{\mathrm{BE}}`$ and $`N_{\mathrm{core}}`$ are evaluated (i.e., $`i=1,\mathrm{},m`$). We have chosen to minimize the difference of the logarithms in order to give the same weight to the central and outer parts of the cores, even though they have very disparate typical column densities.

The ranges of central density and temperature values used to find the best BE fit are, for $`\rho _{0\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}}`$ in eq (4), from 0.1 to 10 times the actual value of the central density, with linear increments of size 0.1 below this density, and of size 1 above it. For the temperature, we actually vary the value of $`\xi _{\mathrm{max}}`$, from 1 to 20 with increments of 0.1, which is equivalent to varying the sound speed, and thus the temperature (see eq. ).

A few comments regarding the choice of method are in order. Our method follows as closely as possible those used in recent observational works (e.g., Johnstone et al. 2000; Evans et al. 2001; Alves et al. 2001; Harvey et al. 2001), although in general not much detail is given in those papers. Other choices for the fitting procedure are clearly possible. In fact, when analyzing the simulations there is more freedom in choosing the procedure, because we have much more information about the physical fields than is available observationally. Observational data generally cover only the core itself, and extend only out to where the data become noisy. In the simulations, one has information on the density field out to distances up to two orders of magnitude larger than the core, and moreover, there is no instrument noise. This means that there is no clear boundary to the cores. Thus, the choice of the truncation radius becomes more arbitrary. Our criterion, based mostly on the shapes of the profiles, is probably as arbitrary as one based on where the data become too noisy. We could have instead chosen to truncate our fits where the column density contrast reaches a certain threshold, but this would probably bias our results towards BE fits with a certain $`\xi _{\mathrm{max}}`$ (as is probably the case in some of the observational papers). Alternatively, one could fix the central density and temperature, and then use $`R_{\mathrm{core}}`$ and the external pressure as the fitting parameters, but this is not what is done observationally. Keeping in mind our goal of showing that those works do not provide unambiguous evidence of the hydrostaticity of the cores, we consider that our conventions are adequate.

## 4 Results

### 4.1 Statistical Analysis

We have analyzed the first 21 cores found by the clump finding algorithm, at two different timesteps of each one of the three simulations. Since we have projected each core on its 3 directions, we have analyzed a total of 378 column density profiles. We stress again that none of these cores are in hydrostatic equilibrium. Instead, they either get dispersed or collapse (see §4.3).

Table 2 lists the results of the analysis for one core in each physical situation, chosen to illustrate the variety of situations that can occur. Column 1 gives the name of the run (see §2 for its physical parameters.) Columns 2, 3 and 4 denote the time at which the run was analyzed, the order number of the core, and the projection plane, respectively. Columns 5 and 6 denote the central density of the core ($`n_{0\mathrm{core}}`$) and the central density of the fitted BE sphere ($`n_{0\mathrm{fit}}`$), respectively. Columns 7 and 8 respectively give the temperature ($`T_{\mathrm{fit}}`$) and the value of the nondimensional radius $`\xi _{\mathrm{max}}`$ for the best core fit. These values have been obtained using a radius ($`R_{\mathrm{core}}`$) listed in column 9. Column 10 gives the value of the rms error, given by eq. (5), between the fit and the actual angle-averaged column density profile, $`\chi _1^2^{1/2}`$. Similarly, in column 11 we give the rms separation between the “dispersion” profile and the angle-averaged one, $`\chi _2^2^{1/2}`$.

From the 378 projections, we reject those that have fits whose rms distance to the column density profile ($`N_{0,\mathrm{core}}`$) is larger than the rms distance of the dispersion curves to $`N_{0,\mathrm{core}}`$, i.e., if $`\chi _1^2^{1/2}>\chi _2^2^{1/2}`$. Thus, only 65.08% of the 378 projections are well described by a BE profile, for which the values of the ratio $`\chi _1^2^{1/2}/\chi _2^2^{1/2}`$ are plotted against $`\chi _1^2^{1/2}`$ in Fig. 1. From this figure, we note that the bulk of the fitted profiles are roughly a factor of 2 or 3 closer to the mean column density profile than the dispersion curves.

In Fig. 2 we show the number of projections that have good BE fits, for all run types and times considered. We note that, for runs GC and LSD, there are fewer cores with BE-like fits at $`t_1`$ than at $`t_0`$. In contrast, for run SSD, there are substantially more BE-like cores at $`t_1`$ than at $`t_0`$. These results suggest that the structure of the cores depends sensitively on the global parameters and initial conditions of the flow. Indeed, at $`t=`$$`t_0`$ the cores in runs LSD and SSD are due exclusively to advection, the ones in run SSD having more substructure (Mac Low & Ossenkopf, 2000; Ballesteros-Paredes & Mac Low, 2002), because the smaller-scale driving implies that there is more turbulent energy at small scales in this run than in LSD. Run GC, on the other hand, was started with an already relatively smooth initial density distribution and furthermore has been subject to the action of self-gravity from the start, so it too, as run LSD, has little substructure within the cores. Thus, a BE profile (with $`\xi _{\mathrm{max}}20`$) is a better fit for the smoother cores of runs GC and LSD than for the more irregular cores in run SSD at $`t=`$$`t_0`$.

At $`t_1`$, instead, gravity has already been acting for some time on all runs, and has produced both smoothing and collapse of the cores. In this case, the originally smoother cores in runs LSD and GC have already reached advanced stages of collapse, and fitting them would require larger values of $`\xi _{\mathrm{max}}`$ than we have considered here, while the cores in run SSD, which started out more turbulent, had to first overcome that turbulence (either by dissipating it or by accreting more mass) and only then could start collapsing. Thus, they appear to be in an earlier collapse phase with more of them in the range $`1\xi _{\mathrm{max}}20`$. Note that, however, the above results do not imply that less mass has been accreted into collapsed objects in run SSD compared to the other two, since $`t_1`$ is defined as the time at which this mass is 5% of the total. Instead, it is simply a reflection of the fact that collapse proceeds at a higher rate in runs LSD and GC (Klessen et al., 2000).

We now turn to the analysis of the physical parameters as obtained from the fitting procedure. Fig. 3 shows the ratio of fitted-to-actual central density, $`n_{0\mathrm{fit}}`$/$`n_{0\mathrm{core}}`$, as a function of $`n_{0\mathrm{core}}`$ at times $`t_0`$ (left panel) and $`t_1`$ (right panel). We find that the central density obtained with the BE fit tends to be smaller than the actual density of the core.

In Fig.’s 4a and b, we show the distribution of the fit parameters $`\xi _{\mathrm{max}}`$, and $`T_{\mathrm{fit}}`$. The contribution of each one of the six physical situations is represented with a different gray-scale tone. In the first histogram, the vertical dotted line at $`\xi _{\mathrm{max}}`$ $`=`$ 6.5 denotes the critical value above which a BE density profile corresponds to an unstable equilibrium. We note that roughly half (47.2%) of the fits have $`\xi _{\mathrm{max}}`$$`6.5`$, making the cores to appear as stable configurations. In Fig. 4b the histogram exhibits a broad distribution of fit temperatures, ranging between 5 K and 60 K, although most of them have temperatures between 5 and 30 K. Recall that the numerical simulations are isothermal, with the scaling taken such that $`T=11.3`$ K. Thus, we see that the temperature derived from the BE fit procedure does not recover the actual temperature of the core very well. Finally, in Fig. 4c (right), we show that most of the fits extend to radii $`R_{\mathrm{core}}`$ between 0.04 and 0.1 pc, but there are several cores whose BE density distributions extend up to $`0.13`$ pc. These values are similar to the values of the BE fits shown in the literature (Shirley et al., 2000; Johnstone et al., 2000; Alves et al., 2001; Evans et al., 2001; Langer & Willacy, 2001; Harvey et al., 2001).

An important remark is that, amongst the 248 accepted fits, 195 come from 65 cores that can be fit in all three directions ($`x`$-$`y`$, $`x`$-$`z`$, and $`y`$-$`z`$) simultaneously. Out of the remaining fits, 40 come from 20 cores that are fitted in only two of their three projections, and 11 cores resemble BE spheres in only one projection. Furthermore, cores that are accepted in more than one projection do not generally yield the same values for $`T_{\mathrm{fit}}`$, and/or $`\xi _{\mathrm{max}}`$. This fact is shown in Fig. 5, where we give the values obtained for the temperatures (upper 3 panels) and dimensionless radii of the BE sphere, $`\xi _{\mathrm{max}}`$ (lower 3 panels). Under the light of these results, the usefulness of BE-type fits to molecular cloud cores is seen to be suspect. We discuss this issue further in §5.

### 4.2 Column density profiles and projection effects

We now turn to the column density structure of the cores and to the confusion that may arise due to projection effects. In Figures 68 we show column density maps and radial profiles of the six cores listed in Table 2. Upper panels show the logarithmic column density profiles, and lower panels show the column density maps in logarithmic gray-scale. In the first, the dashed line corresponds to the actual averaged density profile, dotted lines are the average $`\pm `$ the standard deviation column density profiles (the “dispersion” profiles), and the solid line denotes the column density of the fitted BE profile. The left, middle, and right panels respectively depict the $`x`$-$`y`$, $`x`$-$`z`$, and $`y`$-$`z`$ projections respectively. The white circle in the lower panels shows the size of the BE sphere fitted. The first point to notice is that the column density distribution within these circles is far from circularly symmetric. Instead, it is highly irregular, often elongated, and sometimes contains more than one local maximum. Nevertheless, radial profiles frequently appear soft and monotonically decreasing.

In Figs. 6a and 6b we show the maps and profiles for clump 0 in SSD at $`t_0`$ and clump 13 in SSD at $`t_1`$, respectively. For each of the three projections of the first case we fit a BE sphere with relatively good confidence at $`R_{\mathrm{core}}`$$`=`$ 0.03, 0.08 and 0.04 pc. We notice that this core has very different values of $`\xi _{\mathrm{max}}`$($`=`$ 4.1, 18.7, and 9), and of the estimated temperature $`T_{\mathrm{core}}`$($`=`$ 20.59, 11.73, and 17.73 K) for each of the projections. In the case of SSD at $`t_1`$, the values obtained using $`R_{\mathrm{core}}`$$`=0.08`$, 0.025, and 0.025 pc are $`\xi _{\mathrm{max}}`$$`=`$ 20, 3.4, and 7.5; and $`T=`$ 13.68, 15.41, and 14.25 K. In particular, note that in the case of SSD at $`t_1`$, the three values for $`\xi _{\mathrm{max}}`$ are such that the BE fit will imply unstable, stable, and unstable (but closer to critical) configurations for the $`x`$-$`y`$, $`x`$-$`z`$, and $`y`$-$`z`$ projections, respectively.

Regarding run LSD, we analyze clump 5 at $`t=t_0`$, and clump 19 at $`t=t_1`$. For the former (Fig. 7a), we find again different BE configurations for each projection; in the $`x`$-$`y`$ and $`x`$-$`z`$ projections we obtain very small values of $`\xi _{\mathrm{max}}`$(3.3 and 4.1 respectively, implying a core in stable equilibrium), and the core is fitted up to only 0.02 and 0.04 pc respectively, while in the $`y`$-$`z`$ projection we find $`\xi _{\mathrm{max}}`$ = 8.4, and the fit can be made out to 0.08 pc. The fitted temperatures for this core are 24.47, 21.13, and 25.17 K. The profiles and the maps for clump 19 in LSD at $`t_1`$ are shown in Fig. 7b. This core exhibits an elongated structure in each projection. The fitted values of $`\xi _{\mathrm{max}}`$ are 8.4, 15.1, and 9.6, respectively. The temperatures in each projection are more scattered than in the previous examples, being $`T=25.95`$, 10.71, and 33.12. It is convenient to mention that this core appears in this figure not because of the great quality of its fit, but to stress the problem of column density maps. For instance, while the maximum of the column density is located close to the maximum of the volumetric density for the $`x`$-$`y`$ and $`y`$-$`z`$ projections, the maximum column density of the $`x`$-$`z`$ projection is located almost at the edge of the box. In fact, we have found that frequently, the position of the maximum column density is shifted by several pixels from the position of the maximum volumetric density (1 pixel $`=`$ 0.0012 pc for LSD and SSD, and 0.0014 pc for GC), the $`x`$-$`z`$ projection being the most critical case.

Finally, we show cores 4 and 26 for run GC at $`t_0`$ and $`t_1`$ respectively. In this run, cores tend in general to be more roundish than in the other cases and thus have smaller error. Figure 8a illustrates the effects of projection. The $`x`$-$`z`$ projection (middle panel) shows a well-defined, round core, and is very well fitted by a BE sphere with $`\xi _{\mathrm{max}}`$ $`=`$ 5.1. Nevertheless, this core is just the overlap of two cores seen in projection, as shown by the $`x`$-$`y`$ and $`y`$-$`z`$ projections. A similar case is presented for GC at $`t=t_1`$ in Fig. 8b, where the $`x`$-$`y`$ projection gives an unstable fit ($`\xi _{\mathrm{max}}`$=8.6), while the other two projections give stable configurations (5.7 and 2.7 for the $`x`$-$`z`$ and $`y`$-$`z`$, respectively).

Before ending this section, there are some points worth noting. First of all, our criterion for a “good” fit is that $`\chi _2^2^{1/2}>\chi _1^2^{1/2}`$ (see §4.1). The algorithm was constructed to find the maximum column density in the box, which some times is shifted from the center of the box. This is because, even if the maximum volumetric density is well-centered, there is no guarantee that the column density may will be located at the same place. In particular, there are some cores that fall at the edge of the box (as mentioned before, Fig. 7b, projection $`x`$-$`z`$), and others for which $`R_{\mathrm{core}}`$ goes out of the box boundaries (Fig.’s 6a projection $`x`$-$`z`$, 6b projection $`x`$-$`y`$, and 7a projection $`y`$-$`z`$). In these cases, the radial profile was computed by considering only the density structure inside the box.

### 4.3 The physical conditions in the clumps

Throughout this paper, we have repeatedly stated that the cores and clumps in the simulations are not hydrostatic, but are instead dynamic, transient entities. This is obviously the case of the cores selected at $`t_0`$ in runs LSD and SSD, since self gravity has not been turned on at this stage yet. So, the fact that we have many acceptable fits at this time is a clear proof that BE fits to the cores’ column density profiles do not provide unambiguous evidence of the cores being in hydrostatic equilibrium.

It is moreover instructive to analyze the cores’ physical structure, and compare it to the observational data. A detailed study of the clump evolution and a comparison between the kinematics and line profiles of the simulations and the observations will be presented in a future contribution. Nevertheless, here we wish to give just a brief discussion of the the density and velocity profiles along the three coordinate axes for each one of the six cores presented in the previous section, in order to show their similarity with observed cores, and further justify our claim of non-hydrostatic conditions within them.

Figure 9 shows density (solid lines) and velocity (dotted lines) cuts for the clumps presented in Fig. 6 (clumps 0 at $`t_0`$ and 13 at $`t_1`$); Fig. 10 shows similar cuts for the clumps presented in Fig. 7 (clumps 5 at $`t_0`$ and 19 at $`t_1`$); and Fig. 11 shows cuts for clumps 4 at $`t_0`$ and 26 at $`t_1`$, presented in Fig. 8. Thin lines denote $`x`$ axis cuts of the density (solid) and the $`x`$-component of the velocity field ($`v_x`$, dotted). Similarly, intermediate bold lines denote cuts along the $`y`$ axis of the density and $`v_y`$, while thick lines denote $`z`$-axis cuts of the density and $`v_z`$. From these figures, various points are worth noting. First, the density profiles are asymmetrical at least in one of the projections. Second, the velocity profile across the clump exhibits, at least in one of the three cuts, a difference larger than the sound speed (=0.1, in code units). Third, the velocity gradients exhibited by the cores show that, (a) in the SSD cores (Fig. 9), there are some directions of contraction (negative gradients of the velocity), and some directions of expansion (positive gradients). (b) The same occurs for clump 19 at $`t_1`$ in LSD, but in the case of clump 5 at $`t_0`$ for LSD, the three directions show positive gradients of the velocity field, suggesting that this clump is actually re-expanding. (c) The cores in GC always show negative gradients, a natural result because these runs are not turbulent, and gravitational contraction is the only possible mechanism to form condensations.

The shapes and amplitudes of the velocity profiles across the clumps show that they are qualitatively similar to observed cores, in the following senses: a) The simulated clumps are transonic, exhibiting both sub- and super-sonic velocity differences across them, similarly to the reported velocity dispersions for observed cores (see, e.g., Jijina, Myers, & Adams, 1999). b) The velocity patterns include both inflow and outflow, which also occurs in real cores when observed at high enough resolution (Myers, Evans & Ohashi 2000, and references therein). Thus, our dynamic cores are not dissimilar to observed molecular cloud cores. Nevertheless, observation of the evolution of the simulations clearly shows that the cores are not hydrostatic structures, having typical lifetimes of the order of their crossing times (Kleesen & Lin 2003; see also Vázquez-Semadeni et al. 1996).

## 5 Summary and Discussion

### 5.1 Bonnor-Ebert fits

We have analyzed cores in SPH simulations of molecular clouds from Klessen et al. (1998), Klessen & Burkert (2000, 2001) and Klessen et al. (2000). The advantage of using SPH simulations over, for example, a regular fixed-grid Eulerian code, is that in the SPH case, density enhancements are better resolved spatially, allowing us to make a detailed study of cores at scales between 3$`\times `$10-3 pc and 0.3 pc, the relevant scales for molecular cloud cores.

We have found cores whose angle-averaged column density profiles are well fitted by BE profiles in a variety of physical situations: with and without self-gravity (respectively, times $`t_1`$ and $`t_0`$), with turbulent driving at large and small scales (respectively LSD and SSD), and even cases with merely random, Gaussian initial density fluctuations (case GC). Several results are found. First, we have shown that 65% of the column density profiles studied may resemble BE profiles, in spite of the fact that the cores are not in hydrostatic equilibrium, and the fact that at $`t=t_0`$ the self-gravity has not been included yet for the turbulent runs LSD and SSD.

Second, we have found that BE fits give temperature values in a range between 5 and 60 K, with most of them being between 5 and 30 K. The fitted central densities range between $``$ 104 cm-3 and $``$ 107 cm-3, with most of them between 2$`\times `$105 cm-3 and 106 cm-3. These values are similar to the values found in the literature (Shirley et al., 2000; Johnstone et al., 2000; Alves et al., 2001; Evans et al., 2001; Langer & Willacy, 2001; Harvey et al., 2001). Nevertheless, the fitted values of the temperature and density in general do not represent the values for the actual cores. Specifically, the actual temperature of the cores is always 11.3 K, while the fitted value of the density is typically smaller than the actual density. The values of the non-dimensional radius $`\xi _{\mathrm{max}}`$ range between 2 and 20, with 47% of them below 6.5, and would seem to suggest that an important fraction of the fitted cores are in hydrostatic equilibrium, while in reality they are transient objects.

### 5.2 Goodness of fit, and the stability of B68

A natural question refers to the goodness of our fits. We have proposed that two conditions must be satisfied in order for us to accept a fit. First, the core under consideration should not contain a collapsed object, i.e. it should not yet have formed a sink particle in our numerical scheme, as we are not able to resolve such object. Second, the rms separation between the fitted and the angle-averaged (“mean”) column density profiles must be smaller than the rms separation between the dispersion profile and the mean profile. These two conditions are met in about 1/2 of the analyzed cores. The typical dispersion curves of the fits for the accepted cores are not significantly different from the error bars of several observational studies (Shirley et al., 2000; Johnstone et al., 2000; Alves et al., 2001; Evans et al., 2001; Langer & Willacy, 2001; Harvey et al., 2001) and, like in those studies, our fits generally span one to two orders of magnitude in column density. One notable exception is the fit of Bok globule B68 by Alves et al. (2001), which has remarkably small error bars. In that work, the column density map is defined over $`1000`$ positions (the positions of the background stars), which corresponds to an equivalent resolution of $`\sqrt{1000}30`$ pixels per dimension on the plane of the sky. The reported error bar at each radius in that paper is calculated as the standard deviation of the observational uncertainties of all positions at such radius (J. Alves, private communication). In fact, those authors comment that their column density profile is the highest signal-to-noise radial column density profile ever obtained for Barnard 68. For comparison, in our simulations we know the value of the column density without uncertainty over more than $`128^2=16,384`$ positions on the projection plane, but our rms errors are nevertheless larger, except for some cases in the GC runs. This appears to suggest that B68 may actually be an especially smooth and roundish structure, at least as seen from Earth. Although this is consistent with the fact that B68 is actually located within an H II region, and so the BE paradigm of thermal pressure confinement by a hotter, more tenuous medium, is applicable in this case, there are still some pieces of inconsistency concerning this core: First, the core is not round, or elliptical. This makes it implausible that it can be in precise hydrostatic equilibrium. Second, their recently reported observations show motions of 0.25$``$0.5 the sound speed, suggesting that the cloud is near but not precisely in hydrostatic equilibrium. And third, the fitted Bonnor-Ebert profile by Alves et al. (2001) implies a temperature of 16 K, and, even at this temperature, marginal instability. The actual temperature is $`10`$ K, i.e., a factor of 30-50% lower that that (Hotzel, Harju, & Juvela, 2002), and thus the thermal support is even lower. Thus, even within the context of BE spheres, this core is unstable by a wide margin, and therefore the hydrostatic model does not appear to provide a satisfactory explanation of its physical state. In the next section we suggest a possible alternative for B68.

### 5.3 Projection effects

We have shown that the cores in the simulations are in general far from spherical, similarly to the situation for actual molecular cloud cores, and that they may exhibit substantially different column density profiles, depending on the projection direction. Cores often exhibit different morphologies in different directions, and the radius at which they seem to merge with their surroundings is not unique, also depending on the direction of the projection. The projection effects play an important role here, and a core apparently well-defined in one projection may appear as a highly structured column density profile in the others. An extreme case is presented by Boss & Hartmann (2001), who show that a reasonable Bonnor-Ebert fit can be obtained even to a disk-like structure seen edge-on.

In this regard it is important to mention that, in order to minimize confusion of cores that are in the same line of sight (LOS) but whose separation is large enough to be considered dynamically disconnected from each other (as was often the case in the structures analyzed by Ballesteros-Paredes et al. 1999; Ballesteros-Paredes & Mac Low 2002), we have chosen sub-boxes of size one-tenth of the full computational domain. Thus, the sub-boxes have sizes of 0.154 pc in the SSD and LSD cases, and of 0.18 pc in the GC case, around the center of the three-dimensional core, and the LOS integration for producing the column density maps was performed only over the length of these sub-boxes. Since, even within these small length scales there is substantial substructure in the simulations, it is reasonable to expect that observed regions of comparable size in molecular clouds should also contain significant amounts of substructure, the apparent observational smoothness of some cores possibly being the result of LOS crowding, as occurs for example in Fig. 8.

The possibility that there is substructure at such small scales is consistent with the fact that most newborn stars belong to multiple systems and thus they must be born in structured, non-uniform cores. Observational evidence that molecular cores at these scales can have substantial clumpy structure is given by Velusamy, Kuiper & Langer (1995) and Wilner et al. (2000), who show that B335, recognized earlier as one of the best candidates of an isolated, round globule in a collapsing phase (see Myers, et al. 2000, and references therein), exhibits a complex, asymmetrical structure, and where the physical conditions of the infalling gas suggest that the standard model of protostellar collapse fails.

This possibility, together with inspection of Fig. 8, suggests an interesting alternative for the nature of B68. In this figure, we show one sub-box that has an exceptionally good BE fit on the $`x`$-$`z`$ projection plane, while on the other two projection planes it is seen to consist of two density peaks. Thus, the goodness of the BE fit to B68 by Alves et al. (2001) is not an unambiguous proof of the core’s BE-like nature, even though, for this particular core, this is not unlikely either.

We conclude that the evidence based on the BE fitting procedure that cores in molecular clouds are in hydrostatic equilibrium is inconclusive. In order to discriminate between the standard picture of low-mass star formation that proposes that cores are quiescent and the turbulent picture that states that they are dynamical, transient entities, we need more detailed observations and theoretical work, although we emphasize that, to date, there is no numerical model which allows for a turbulent molecular cloud that has produced quiescent cores222Note that some of the simulations by Ostriker et al. (1999) do contain magnetostatic flattened structures of sizes comparable to the whole computational domain. These are the result of super-Jeans but sub-critical initial conditions, and periodic boundary conditions that do not allow for further mass accretion along field lines. However, we expect that in actual molecular clouds, cores can always accrete matter, eventually becoming supercritical, and proceeding to collapse., while the cores analyzed here do show similar physical conditions to those typically reported for real cores.

We thank L. Hartmann for careful reading of the manuscript. EVS and JBP acknowledge support from CONACYT’s grants 27752-E and I39318-E respectively. RSK acknowledges support by the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (DFG: KL1358/1) and funding by a NASA astrophysics theory program supporting the joint Center for Star Formation Studies at NASA-Ames Research Center, UC Berkeley, and UC Santa Cruz. This research has made use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Phase separation kinetics in compressible polymer solutions: Computer simulation of the early stages

## 1 Introduction

The phase behavior of polymer-solvent systems has important application in the industry for the production and processing of many kinds of plastic materials . An example is the formation of solid polystyrene foams . In addition, both the understanding of the equilibrium phase diagram of these systems and the kinetic mechanism of phase separation are challenging problems of statistical mechanics. On the one hand, variation of the molecular weight of the polymer (i.e. the “chain length” of the flexible linear macromolecule, N) offers a control parameter that leaves intermolecular forces invariant, and hence allows a more stringent test of the theory than would be possible for small molecule systems. On the other hand, it is crucial to take the compressibility of the system fully into account, since for practical reasons one normally wishes to use supercritical fluids (e.g., CO2) as a solvent . Then small changes of the pressure result in a large variation of the solvent density, and this fact obviously facilitates applications. While for an incompressible polymer solution (as it is modelled, for instance, by the well-known Flory-Huggins lattice theory ) the early stages of phase separation can only be studied if one realizes a rapid temperature quench from a state in the one phase region into the miscibility gap (cf. Fig. 1a), for compressible polymer solutions it is possible to induce phase separation by an (experimentally much more convenient) pressure jump (cf. Fig. 1b). Of course, in Fig. 1 we have made the tacit assumption that including the pressure does not change the qualitative form of the phase diagram, i.e., the point ($`T_c,x_c`$) in $`(T,x,p)`$-space becomes a line $`(T_c(p),x_c(p))`$ or, expressed in the appropriate inverse functions, $`(p_c(T),x_c(T))`$, respectively (Fig. 1b). In reality, the situation is more complicated, since one must treat both the densities of CO2 and hexadecane as two coupled order parameters, which both play a crucial role in the possible phase transitions (remember that for the pure solvent at its critical point gas-liquid phase separation starts). As a result, it can happen that both gas-liquid and liquid-liquid phase separation are present in the phase diagram, leading to the occurrence of a triple line (Fig. 2) . As will be discussed below, the presence of the triple line has a profound influence on the interface properties and hence the nucleation behavior in the system.

The present paper is devoted to a study of the early stages of phase separation kinetics in such systems by computer simulation. We use a coarse-grained model for mixtures of hexadecane and carbondioxide as an archetypical example. The reason for the choice of this particular system is the fact that recently extensive experiments on nucleation in this system have been performed . In Sec. 2, we shall briefly describe the model and the simulation techniques. Sec. 3 presents the relevant static equilibrium properties, in particular, the phase diagram and the properties of planar interfaces separating coexisting phases. Sec. 4 describes simulations of quenching experiments under different conditions, that lead either to nucleation and growth or spinodal decomposition. Sec. 5 summarizes our conclusions.

## 2 Model and Simulation Techniques

A fully atomistic simulation of mixtures of hexadecane (C16 H34) and carbondioxide (CO2), which aims at both establishing phase diagrams as a function of the three variables $`T,p`$ and $`x`$, and the study of phase transition kinetics under various conditions, would be an extremely formidable problem. The potentials for the length of the covalent bonds in these molecules (as well as the potentials for bond angles) are very stiff. Hence, an extremely short time step (of the order of about 1 fs) would be required in a Molecular Dynamics (MD) simulation. In a corresponding Monte Carlo (MC) simulation only very small random displacements of atoms would be admissible. Thus, for liquid alkanes it is a quite common and well-established practice to integrate CH2-groups (as well as the CH3 end groups) into “united atoms” , and also to work with a fixed C-C bond length. These simplifications already reduce the necessary computer time by a factor of 100. However, we have estimated that even for such a simplified model computer resources as they are available today still are not yet sufficient. Thus it was chosen to simplify the model even further, representing CO2 by a single pseudo-atom, and representing C16 H34 by a flexible chain of 5 subsequent effective segments, each of which then contains roughly 3 successive C-C bonds. All such effective monomeric units interact with a truncated and shifted Lennard Jones (LJ) potential,

$$V_{LJ}(r)=4\epsilon _{\text{hh}}[(\frac{\sigma _{\text{hh}}}{r})^{12}(\frac{\sigma _{\text{hh}}}{r})^6+\frac{127}{16384}],rr_c=22^{1/6}\sigma ,$$
(1)
while $`V_{LJ}(r>r_c)=0`$. Subsequent effective monomers along a chain are in addition exposed to a finitely extensible nonlinear elastic (FENE) potential

$$V_{\text{FENE}}(r)=33.75ϵ_{hh}\mathrm{ln}[1(\frac{r}{1.5\sigma _{\text{hh}}})^2]$$
(2)
The constants in this potential are chosen such that the most favorable distance between bonded neighbors is 0.96 $`\sigma _{\text{hh}}`$, while the preferred distance between non-bonded effective monomers is $`2^{1/6}\sigma _{\text{hh}}1.12\sigma _{\text{hh}}`$. This mismatch is desirable to prevent crystallization at high densities , which is appropriate for glass-forming polymers .

The parameters $`\sigma _{\text{hh}}`$, $`\epsilon _{\text{hh}}`$ set the scales for energy and length of our hexadecane model. Since we wish that our model fits the thermodynamic properties of this material in the fluid phase as faithfully as possible, we have adjusted them such that the critical temperature $`T_c=723`$ K and the critical density $`\rho _c=0.21`$ g/cm3 of hexadecane are correctly reproduced. In simulations, a critical point can be computed by a finite size scaling study of the liquid-gas transition along the lines of the techniques proposed by Wilding et al. . A comparison of critical temperature and density in Lennard-Jones units with the experimental values yields $`\epsilon _{\text{hh}}=5.78710^{21}`$ J, $`\sigma _{\text{hh}}=4.52310^{10}`$m .

Since the description of hexadecane is thus already reduced to a crudely coarse-grained model, it would not make sense to keep all atomistic detail for carbon dioxide. Thus the CO2 molecule is also coarse-grained into a point particle, and we require an interaction potential between CO2 molecules of exactly the same LJ form as in Eq. (1), but with parameters $`\epsilon _{\text{cc}}`$, $`\sigma _{\text{cc}}`$. Requiring once more that the critical temperature and density of real CO2 are correctly reproduced we obtained  $`\epsilon _{\text{cc}}=4.20110^{21}`$ J, and $`\sigma _{\text{cc}}=3.69310^{10}`$ m. It is clear that this procedure ignores some physical effects, e.g., the fact that the CO2 molecules carry electrical multipole moments is completely neglected. Nevertheless our model not only describes the gas-liquid coexistence curve of hexadecane and CO2 over a reasonable temperature range , but also other data (like the critical pressure, or the temperature dependence of the surface tension near $`T_c`$ ).

Of course, the choice of interaction parameters between CO2 and hexadecane is more subtle. We model the interaction between the pseudoatoms representing CO2 and the pseudoatoms representing three subsequent CH2 groups again by a LJ potential of the form Eq. (1), but now with parameters $`\epsilon _{\text{hc}}`$, $`\sigma _{\text{hc}}`$. It then remains to find an optimal choice for these parameters. The simplest choice would be the well-known Lorentz-Berthelot mixing rule

$$\epsilon _{\text{hc}}=\sqrt{\epsilon _{\text{hh}}\epsilon _{\text{cc}}},\sigma _{\text{hc}}=(\sigma _{\text{hh}}+\sigma _{\text{cc}})/2.$$
(3)
However, when one tries the choice Eq. (3) the resulting phase diagram of our model system would be of type I in the classification scheme of Van Konynenburg and Scott  (i.e., the critical points of pure hexadecane and CO2 are connected by a critical line of the mixture system. Liquid-liquid phase seperation does not exist and thus there is no three phase coexistence). However, experimentally it is known  that alkane – CO2 mixtures exhibit a type I – phase behavior only for very short alkanes, while for hexadecane – CO2 mixtures the phase diagram is of type III. Instead of a connecting line we rather observe a topology as shown in Fig. 3 . In the (p,T) projection of the phase diagram, a critical line emerges from the critical point of hexadecane and does not end at the critical point of CO2. In addition we observe liquid-liquid immisciblity and a three-phase line. When we want to obtain this behavior from our coarse-grained model we must allow for a deviation from the first equation of Eq. (3), assuming rather

$$\epsilon _{\text{hc}}=\xi \sqrt{\epsilon _{\text{hh}},\epsilon _{\text{cc}}},\xi <1.$$
(4)
It turns out that a suitable choice for the parameter $`\xi `$ which characterizes the deviation from the Lorentz-Berthelot mixing rule is  $`\xi =0.886`$. Eq. (4) with this choice of $`\xi `$ was in fact used to compute the phase diagram shown in Fig. 3. However, experimental measurements of the critical line vary considerably and a small modification in $`\xi `$ rises or lowers the critical line. Details about how such phase diagrams are in fact estimated from the simulation are given elsewhere .

For a study of the phase diagram (Fig. 3) as well as for a study of the interfacial free energy between coexisting gas and liquid phases  the grandcanonical ensemble is used, where volume $`V`$, temperature $`T`$, and the chemical potentials $`\mu _\text{c}`$, $`\mu _\text{h}`$ are fixed. Both the particle numbers $`N_\text{c}`$ and $`N_\text{h}`$ of CO2 and hexadecane and the pressure $`p`$ are then (fluctuating) observables “measured” in the simulation . Fig. 4 shows an isothermal slice through the phase diagram at $`T=486`$ K. This temperature is higher than the temperature of the critical end point where the critical line ending at the critical point of CO2 and the line of triple points ($`p_{\text{trip}}(T)`$) in the ($`p,T`$) phase diagram (Fig. 3) meet. Therefore unlike Fig. 2 there is no triple point in the phase diagram, and the qualitative features are the same as in Fig. 1b. At this temperature we shall examine the kinetics of the phase separation in Sec.4.

In order to observe the kinetics of bubble nucleation in the vicinity of the binodal we control the undersaturation of the system by fixing the chemical potential of both species. As a starting configuration we use a homogeneous state equilibrated at a higher temperature. This situation corresponds to a simulation cell in contact with a much larger reservoir which is held at constant undersaturation. As the number of particles in the simulation cell is not conserved, and the particles do not move according to a realistic dynamics, we do not obtain information about the time scale of bubble nucleation. The kinetics of phase separation in the vicinity of the binodal is, however, chiefly determined by the free energy barrier the system encounters on its path towards the equilibrium state. Hence, we expect the relaxation path to be similar to a simulation with realistic dynamics. The most natural choice would be to apply Molecular Dynamics methods . However, if the number of particles was conserved we would need prohibitively large simulation cells for the undersaturation not to decrease substantially even in the very early stages of nucleation. The simulation of spinodal decomposition took place in the canonical ensemble. Again, a homogeneous starting configuration was equilibrated at a higher temperature and quenched to a state below a critical point of the mixture system. We only allowed for local Monte Carlo displacements which should yield kinetics comparable to those obtained in Molecular Dynamics.

## 3 Static equilibrium properties

First we return to the phase diagrams in Figs. 34 and discuss their properties more closely. Note that at the temperature chosen in Fig. 4 ($`T=486`$ K) we are very far below the critical point of pure hexadecane (which is $`T_c=723`$ K). Therefore the hexadecane melt (with some dissolved carbon dioxide) coexists with almost pure CO2 the composition of which gradually changes from gas- to liquid-like with increasing pressure. In both cases there is almost no hexadecane present. This is reflected in the coexistence curve which rapidly approaches the CO2 molar fraction $`x=1`$.

We have also included the results of an analytical approximation for the phase diagram of our model, obtained from a thermodynamic perturbation theory  called TPT1, developed along the lines of Wertheim et al. . It is seen that this theory describes the coexistence curve in Fig.4 at pressures that are much smaller than the critical pressure $`p_{\text{crit}}(T)`$ almost quantitatively, while the critical pressure itself is clearly overestimated. Of course, such an overestimation of the critical parameters $`p_{\text{crit}}(T)`$, $`x_{\text{crit}}(T)`$ is quite typical for all theories that invoke a mean field description of the critical behavior, as TPT1 does. A similar analytical study for more alcane-CO2 systems has been performed recently . Note also that the TPT1 theory readily yields a spinodal curve, which has the standard meaning of separating “metastable states” from “unstable states” in the phase diagram . The naive (mean-field) description implies that in the metastable region of the phase diagram the initial stages of phase separation kinetics are described by nucleation and growth , while in the unstable region the decay mechanism is spinodal decomposition . In systems with short range forces there is no well-defined sharp spinodal line, the actual transition between both decay mechanisms occurs rather gradually in a broad transition region, and this region is not centered around the mean field spindodal but occurs close to the coexistence curve . However, for the quenching experiment performed at $`T=486`$ K and $`x=0.60`$, where the final state corresponds to a pressure of about $`p`$130 bar, we are so close to the coexistence curve and so far from the mean field spinodal, that classical nucleation-type behavior  should be observable. In order to demonstrate the mechanism of spinodal decomposition, we quenched into the unstable region as defined by the mean-field spinodal (compare with Fig. 4).

Classical nucleation theory assumes (in our case) the formation of spherical gas bubbles of essentially pure CO2 within a hexadecane matrix . For a quantitative understanding of the free energy barrier against homogeneous nucleation, clearly the understanding of the interface between macroscopic coexisting phases is a prerequisite. As discussed elsewhere , such interfaces are most conveniently studied in the grandcanonical ensemble in conjunction with a multicanonical preweighting scheme that generates mixed-phase configurations. Choosing a simulation box of rectangular shape $`L\times L\times cL`$ with a constant $`c=2`$ or larger and periodic boundary conditions, one generates system configurations where the two coexisting phases are separated by two parallel $`L\times L`$ interfaces (compare with Fig. 5a). “Measuring” the canonical probability of these two-phase configurations relative to the probability of the pure phases is a standard method for the estimation of interfacial free energies . The results for T=486 K are shown in the inlet of Fig. 4. In addition, this method can be used to generate well-equilibrated configurations of interfaces allowing a study of the interfacial structure. Fig. 5 shows, as an example, two snapshot pictures of such states containing interfaces, one at $`T`$=486 K and the other at $`T`$=243 K. The size of the particles is enlarged (radii shown are $`\sigma _{\text{cc}}`$, $`\sigma _{\text{hh}}`$ rather than $`\sigma _{\text{cc}}/2`$, $`\sigma _{\text{hh}}/2`$), for the sake of a clearer view. The temperature $`T`$=243 K is slightly lower than the temperature of the critical point of pure CO2, and thus the system possesses a triple point where a polymer-rich phase coexists with two CO2-rich phases, one being liquid CO2 and the other the gas. Consequently, when we study the interface between the polymer-rich phase and gaseous CO2, a layer of the third phase (liquid CO2) intrudes at this interface. The thickness of this layer is expected to diverge when the pressure approaches the value of the triple point pressure. For higher pressures, we then have a coexistence between the polymer rich-phase and the liquid-like dense CO2, similar to the situation that occurs at the higher temperature (Fig. 5a). At low enough temperatures, where a triple point occurs, one expects a strong decrease of the interfacial tension already when the triple point pressure is approached : the interfacial tension between the polymer-rich phase and liquid-like CO2 is much smaller than the interfacial tension between the polymer-rich phase and the gas. The barrier against nucleation also decreases when one approaches the spinodal from the metastable region. This fact facilitates nucleation of relatively large droplets observable on the time-scales of our Monte Carlo simulation.

## 4 Monte Carlo simulation of bubble nucleation

The animation in Fig. 6 visualizes how phase separation in polymer solutions proceeds via nucleation. In the beginning the system fluctuates around a metastable free energy minimum. One observes the “birth and death” of small density fluctuations in the hexadecane matrix (displayed as red spheres) which become visible whenever the blue background shines through the slice. These irregular voids are too small, however, to lead to an immediate decay of the initial metastable state. Usually the voids also contain a few CO2 molecules (displayed as small yellow spheres). Only after some time lag a void manages to grow to critical size and overcome the free energy barrier which seperates the metastable from the homogeneous equilibrium state. From now on the bubble will grow until it fills the whole simulation box. As expected, the critical void gets also filled with CO2 molecules, thus decreasing the supersaturation of the remaining polymer-rich phase. Remarkably, this filling of the bubble does not occur homogeneously: at the interface between the gas and the polymer, i.e., at the surface of the bubble, the CO2 density is clearly higher than it is in the center of the bubble. This surface enrichment effect may be interpreted as a precursor to interfacial wetting of the fluid CO2 predicted when one approaches the triple point . A quantitative analysis of this phenomena is in preparation .

Of course, the observation of a nucleation event displayed in Fig. 6 is by no means the first observation of nucleation phenomena by computer simulations: e.g., already Ref.  contains a series of snapshot pictures illustrating nucleation in the two-dimensional lattice gas model. Since then many more elaborate studies have appeared. However, the present work deals with nucleation in a system of industrial relevance, namely bubble nucleation in metastable polymer melts supersaturated with supercritical CO2. In addition, the nucleation mechanism here is quite nontrivial, since the “critical bubble” (which just has the size to be at the nucleation barrier, where there is a 50% chance of growth or decay) is characterized by two variables, its size and the number of CO2 molecules it contains. Similarly the interface is described by nontrivial profiles of two variables, the total density and the relative concentration of CO2. Such a nontrivial coupling between two order parameters has also been predicted by recent self-consistent field calculations based on a qualitatively similar model . Since these mean field descriptions do not take into account the effects of fluctuations, a check by computer simulation methods is clearly warranted.

When we quench to a position below the critical point (Fig. 7), the behavior is rather different: fluctuations in the relative concentration of CO2 get gradually more pronounced everywhere. These fluctuations are not localized, however, in the form of identifiable bubbles, but rather form an irregular percolating network. This is the hallmark of spinodal decomposition . Since the linear theory of spinodal decomposition only holds for systems with long range forces in the very early stages , we do not attempt a more quantitative analysis of this data. Note also that the build-up of these concentration fluctuations is relatively faster than the nucleation of bubbles. Hence the separation of time scales between the structural relaxation and the rate constant of spinodal decomposition is less well established. In the intermediate stages of spinodal decomposition, we also expect hydrodynamic mechanisms to prevail , which are not captured by Monte Carlo moves.

Returning now to the simulation of nucleation phenomena, we ask how one can quantify these observations. A straightforward but tedious method would be to repeat dynamic simulations as shown in Fig. 6 many times in order to extract an estimate of the nucleation rate from the average time it takes to form a bubble or a droplet. Of course, since the nucleation rate is expected to decrease dramatically if the quench depth is decreased slightly, such a direct method works hardly in practice. An indirect method to obtain information on the surface free energy of clusters as a function of cluster size simply uses the final equilibrium states of simulations. An estimate for the free energy as a function of density can be obtained directly from the probability distribution. In this context special care has to be taken of finite-size effects to identify regions of densities which correspond to a single cluster . If we assume that bubbles (or droplets) are spherical, filled with gas and surrounded by a homogeneous liquid of coexistence density (or vice versa), we can assign a radius to each density and hence a free energy estimate to each radius. Fig. 8 indicates that the surface free energy obtained this way is somewhat smaller than expected from the so-called “capillarity approximation”. A simple estimate of the free energy is $`F_s=4\pi R^2\gamma _{\mathrm{int}}`$, where $`R`$ is the radius of the cluster and $`\gamma _{\mathrm{int}}`$ the interfacial tension of flat interfaces as shown in Fig. 5 and derived in Fig. 4. As a consequence, the nucleation barrier is lower than the classical theory of homogeneous nucleation  would predict. Differences in $`\gamma _{\mathrm{int}}`$ also decrease with increasing cluster size and distance to the critical point of the mixture. At this point, however, the findings are still somewhat preliminary, and more work will be necessary to obtain sufficiently accurate results .

## 5 Concluding Remarks

In summary, we have presented a simple model and computational tools which allow us to examine equilibrium properties and kinetics of phase separation in systems of industrial relevance. We were able to produce a two-component mixture phase diagram with a complete critical line which is in qualitative agreement with experiments. For a system at T=486 K we determined the isotherm and interfacial tension of a flat surface as a function of pressure. In a small animations and a series of snapshots we visualized bubble nucleation and spinodal decomposition for this specific system. For the first case we found a precursor of interfacial wetting which was also observed on a flat surface in the vicinity of the three-phase coexistence line. Finally, we found evidence that for small bubbles the interfacial tension is smaller than expected by classical theory of homogeneous nucleation.

## 6 Acknowledgements

Financial support from the BASF AG (P.V.) and from the Deutsche Forschungsgemeinschaft (DFG) via a Heisenberg fellowship (M.M.) and by grant No BT 314/17-3 (L.G.M) is gratefully acknowledged. LGM would like to thank Ministerio de Ciencia y Tecnologia (MCYT) and Universidad Complutense (UCM) for the award of a Ramon y Cajal fellowship and for financial support under contract BFM-2001-1420-C02-01. We are grateful to E. Hädicke, B. Rathke, R. Strey, and D.W. Heermann for stimulating discussions, and to NIC Jülich, HLR Stuttgart, and ZDV Mainz for generous grants of computer time.


text:
# Spin Current Through a Magnetic-Oscillating Quantum Dot

## Abstract

Non-equilibrium spin transport through an interacting quantum dot is analyzed. The coherent spin oscillations in the dot provide a generating source for spin current. In the interacting regime, the Kondo effect is influenced in a significant way by the presence of the precessing magnetic field. In particular, when the precession frequency is tuned to resonance between spin up and spin down states of the dot, Kondo singularity for each spin splits into a superposition of two resonance peaks. The Kondo-type cotunneling contribution is manifested by a large enhancement of the pumped spin current in the strong coupling, low temperature regime.
PACS numbers: 72.25.-b, 85.35.Be, 73.23.-b, 03.67.Lx
Key words: Spintronics, quantum dot, Kondo effect

Introduction.—The rapid progress of nano-electronics and information technologies has prompted intense interest in exploiting the interplay of electron charge and spin properties, which results in the emergence of semiconductor spintronicsDatta ; Prinz . Since it is a manageable task to manipulate and measure electron and nuclear spins via optical or electric fieldKane , the solid-state implementation of a quantum information processing architecture is intensely studiedSamarth , aiming at the possibility of large-scalability of the quantum computer. One of the most important spin-based electronic devices is a mesoscopic quantum dot (QD) system, in which the spin coherence time for electrons or nuclei is relatively longKikkawa ; Kou . Spin-polarized transport through a QD has been extensively investigated recently. It has been shown theoreticallyMucciolo ; Wang  and demonstrated experimentallyWatson  that a QD system will function as a phase-coherent spin pump in the presence of sizable Zeeman splitting. A QD-based spin battery device has been proposed by use of a combination of a vertical magnetic field and an oscillating electric fieldSun1 . Very recently, spin-polarized current has been detected from a quantum point contact (QPC)Potok1  and from Coulomb-blockaded QDsPotok2 . For the latter, due to the large charging energy at low temperatures, a more subtle effect—the creation of new states of many-body character at the Fermi level by the Kondo effect—is expected to occurHewson . Although the effect of Kondo resonance on the charge current has been well studied in a QDKouwenhoven , its influence on the spin current and spin detection might be equally significant, since it provides a spin-flip cotunneling mechanism.

In this paper we explore spin current response to a magnetic-oscillating quantum dot, which is schematically shown in Fig. 1. The setup is similar to the one used in recent ESR-STM experimentsWach . The single-electron levels in the dot are split by an external magnetic field $`B_0`$, $`\epsilon _d\epsilon _d=g\mu _BB_0`$, where $`g`$ is the effective electron gyromagnetic factor and $`\mu _B`$ is the Bohr magneton. The two spin levels are coupled by a rotating magnetic field ($`B_1\mathrm{cos}\omega t`$, $`B_1\mathrm{sin}\omega t`$), applied perpendicularly to the field $`B_0`$. The possibility of the ESR setup based on GaAs QD has been discussed in Ref.. The QD is in some sense similar to an oscillating magnetic dipole moment of a two-level atom. However, we show that the spin transport properties are highly non-trivial for the QD system, due to the many-body excitations and coupling with an external electrode.

Theoretical formalism.—The model Hamiltonian under consideration can be written as

$`H={\displaystyle \underset{\sigma }{}}\epsilon _{d\sigma }d_\sigma ^+d_\sigma +Ud_{}^+d_{}d_{}^+d_{}`$ $`g\mu _BB_1(d_{}^+d_{}e^{i\omega t}+d_{}^+d_{}e^{i\omega t})`$ (1)
$`+{\displaystyle \underset{k\sigma }{}}ϵ_ka_{k\sigma }^+a_{k\sigma }`$ $`+{\displaystyle \underset{k\sigma }{}}\left[V_ka_{k\alpha \sigma }^+d_\sigma +\text{H.c.}\right].`$
Here $`d_\sigma ^+`$ and $`a_{k\sigma }^+`$ create an electron of spin $`\sigma `$ in the QD and in the lead, respectively. The first two terms in (1) represent correlated spin levels of the QD, the time-dependent terms describe the coupling between two QD spin states $`|`$ and $`|`$ by the rotating magnetic field, and the last term is the hybridization of the QD to the lead. A similar model has been proposed for electrical detection of single-electron spin resonanceEngel ; Jiang .

The time dependence of the Hamiltonian (1) may be eliminated by introducing a unitary transformation $`U=e^{i\frac{\omega t}{2}[_k(a_k^+a_ka_k^+a_k)+(d_{}^+d_{}d_{}^+d_{})]}`$ ($`\mathrm{}`$ is set to be unity) and thus redefining the Hamiltonian in the rotating reference (RF) as follows

$`H_{RF}`$ $`=U^1HU+i{\displaystyle \frac{dU^1}{dt}}U`$ (2)
$`={\displaystyle \underset{\sigma }{}}\stackrel{~}{\epsilon }_{d\sigma }d_\sigma ^+d_\sigma +Ud_{}^+d_{}d_{}^+d_{}g\mu _BB_1(d_{}^+d_{}+d_{}^+d_{})`$
$`+{\displaystyle \underset{k\sigma }{}}ϵ_{k\sigma }a_{k\sigma }^+a_{k\sigma }+{\displaystyle \underset{k\sigma }{}}\left[V_ka_{k\alpha \sigma }^+d_\sigma +\text{H.c.}\right],`$
where $`\stackrel{~}{\epsilon }_{d\sigma }=\epsilon _{d\sigma }\pm \frac{\omega }{2}`$ and $`ϵ_{k\sigma }=ϵ_k\pm \frac{\omega }{2}`$ are the dressed dot and lead levels for up and down spins, respectively. One can see that the rotating magnetic field, in effect, shifts lead electron energy to the opposite directions for up and down spins. This spin-oriented energy splitting originates from the hybridization term in (2), which propagates energy from the QD to the lead. Thus the spin-up chemical potential is now $`\mu _{}=\omega /2`$ (The reference energy is set $`E_F=0`$), while the spin-down chemical potential is $`\mu _{}=\omega /2`$. The charge chemical potential is still $`\mu _e=(\mu _{}+\mu _{})/2=0`$. The distinction between $`\mu _\sigma `$ and $`\mu _e`$ was recently discussed in Ref..

The physical picture of spin current generation is as follows. Due to the one-photon difference between the spin-down and spin-up chemical potentials of the lead, an electron in the spin-down channel can tunnel into the $`\stackrel{~}{\epsilon }_d`$ level in QD. After a spin-flip process given by $`B_1`$ term in (2), the electron can tunnel out of the QD into the spin-up channel. This process repeats, so a steady spin current builds up in the lead. Note that the net charge current is zero because there is only one electrode in the system. This zero charge current feature is retained when more electrodes are taken into account with no bias voltage applied between them. We emphasize that the present situation is not a simple realization of magnetization of the electrode by the rotating magnetic field. Here the generation of spin current is due to real charge transfer processes between the two spin channels of the lead. Therefore, the spin current can be conveniently controlled in experiments. To calculate the spin current we make use of the non-equilibrium Green function formalism.

In the absence of hybridization and Coulomb interaction, the QD is reduced to a simple driven two-level system (TLS). Its dynamics is featured by a coherent weight transfer (Rabi oscillations) between the two spin states, which is complete when the rotating frequency is tuned to the resonant condition $`\omega _R=g\mu _BB_0`$. The spin oscillation period is $`T=\pi /\mathrm{\Omega }`$, where $`\mathrm{\Omega }=\sqrt{\mathrm{\Delta }^2+4(g\mu _BB_1)^2}`$ is the Rabi frequency and $`\mathrm{\Delta }=\omega _R\omega `$ denotes the detuning from the resonance.

To proceed we introduce a canonical transformation

$$\left(\begin{array}{c}d_{}\hfill \\ d_{}\hfill \end{array}\right)=𝐮\left(\begin{array}{c}c_{}\hfill \\ c_{}\hfill \end{array}\right)\text{ with }𝐮=\left(\begin{array}{cc}\mathrm{cos}\varphi \hfill & \mathrm{sin}\varphi \hfill \\ \mathrm{sin}\varphi \hfill & \mathrm{cos}\varphi \hfill \end{array}\right),$$
(3)
where $`\varphi =\mathrm{tan}^1(\frac{2g\mu _BB_1}{\mathrm{\Omega }\mathrm{\Delta }})`$. In terms of new fermion operators, the dot Hamiltonian in Eq. (2) is has a diagonal form $`_\sigma \epsilon _{c\sigma }c_\sigma ^+c_\sigma +Uc_{}^+c_{}c_{}^+c_{}`$ with $`\epsilon _{c\sigma }=(\epsilon _d+\epsilon _d\pm \mathrm{\Omega })/2`$ for the up and down spins.

The spin-$`\sigma `$, say spin-$``$, current can be calculated from the time evolution of the occupation number $`N_{}(t)=_kc_k^+c_k`$ for electrons in the electrode using the non-equilibrium Green function:

$$J_{}=e\frac{d}{dt}N_{}(t)=\frac{2e}{\mathrm{}}\mathrm{Re}\left\{\underset{k}{}V_k\frac{dϵ}{2\pi }[\mathrm{𝐮𝐆}_{k,c}^<(ϵ)]_{}\right\},$$
(4)
where we have defined the lesser Green function matrix $`[𝐆_{k,c}^<(t)]_{\sigma \sigma ^{}}=ic_\sigma ^{}^+(0)a_{k\sigma }(t)`$. Next, we use Dyson’s equation to calculate $`𝐆_{k,c}^<`$ and express $`J_{}`$ by the Green functions of the dot as follows

$$J_{}=\frac{ie}{\mathrm{}}\frac{dϵ}{2\pi }\mathrm{\Gamma }(ϵ)\left(𝐮\{𝐆_c^<(ϵ)+f_{}(ϵ)[𝐆_c^r(ϵ)𝐆_c^a(ϵ)]\}𝐮^\text{T}\right)_{},$$
(5)
where $`\mathrm{\Gamma }(ϵ)=2\pi _k|V_k|^2\delta (ϵϵ_k)`$ and $`[𝐆_c^<(t)]_{\sigma \sigma ^{}}=ic_\sigma ^{}^+(0)c_\sigma (t)`$ are the lesser Green functions of the dot, $`𝐀=𝐆_c^r𝐆_c^a`$ is the spectral function matrix, and $`f_\sigma (ϵ)=f(ϵ\mu _\sigma )`$ are the Fermi distribution functions for the spin-up and spin-down channels.

The non-interacting spin current.—In the non-interacting ($`U=0`$) case or within the mean field (MF) treatment of the Coulomb interaction for finite $`U`$, the retarded Green function matrix is given by

$$[𝐆_c^r(ϵ)]_{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}\frac{1}{[𝐠_c^{r1}(ϵ)]_{\sigma \sigma }[𝚺_0^r(ϵ)]_{\sigma \sigma }},$$
(6)
where $`𝚺_0^r=i\mathrm{\Gamma }𝐈`$ is the tunneling contribution to the retarded self-energy, and $`𝐠_c^r(ϵ)`$ is the unperturbed Green function matrix of the dot, which in the mean-field approximation has a simple Hartree-Fock type $`[𝐠^R(ϵ)]_{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}[\frac{n_{c\overline{\sigma }}}{ϵϵ_{c\sigma }U}+\frac{1n_{c\overline{\sigma }}}{ϵϵ_{c\overline{\sigma }}}]`$ with $`n_{c\sigma }=c_\sigma ^+c_\sigma `$. $`n_{c\sigma }`$ must be calculated self-consistently via the relation $`n_{c\sigma }=\mathrm{Im}\frac{dϵ}{2\pi }[𝐆_c^<(ϵ)]_{\sigma \sigma }`$. The lesser Green functions in (5) are obtained via the Keldysh equation $`𝐆_c^<=𝐆_c^r𝚺_0^<𝐆_c^a`$ with the lesser self-energy matrix

$$𝚺_0^<=i\mathrm{\Gamma }\left(\begin{array}{cc}f_{}(ϵ)\mathrm{cos}^2\varphi +f_{}(ϵ)\mathrm{sin}^2\varphi \hfill & \mathrm{sin}\varphi \mathrm{cos}\varphi [f_{}(ϵ)f_{}(ϵ)]\hfill \\ \mathrm{sin}\varphi \mathrm{cos}\varphi [f_{}(ϵ)f_{}(ϵ)]\hfill & f_{}(ϵ)\mathrm{sin}^2\varphi +f_{}(ϵ)\mathrm{cos}^2\varphi \hfill \end{array}\right).$$
(7)
Substituting the expressions for $`𝐆_c^<`$ into (5), we obtain the expression for the tunneling spin current

$$J_{}=\frac{e\mathrm{\Gamma }}{\mathrm{}}\frac{dϵ}{2\pi }\{\mathrm{𝐮𝐆}_0^r(ϵ)\overline{𝚺}_0^<𝐆_0^a(ϵ)𝐮^T\}_{}\{f_{}(ϵ)f_{}(ϵ)\},$$
(8)
where

$$\overline{𝚺}_0^<=\mathrm{\Gamma }\left(\begin{array}{cc}\mathrm{sin}^2\varphi \hfill & \mathrm{sin}\varphi \mathrm{cos}\varphi \hfill \\ \mathrm{sin}\varphi \mathrm{cos}\varphi \hfill & \mathrm{cos}^2\varphi \hfill \end{array}\right).$$
(9)
Since the electrons are flowing from spin-down channel to spin-up channel, we define the total spin current as $`J_s=J_{}J_{}=2J_{}`$. Although the coherent spin superposition is strongest at two-level resonance $`\omega =\omega _R`$, the spin current is also influenced by the spin levels $`\stackrel{~}{\epsilon }_{d\sigma }`$relative to $`\mu _\sigma `$. Figure 2(a) shows the non-interacting spin current versus rotating frequency at different temperatures. The undressed spin levels are set to lie deeply below $`E_F`$. One can see that at low temperature (solid line), the spin current peak is not at, but far from $`\omega =\omega _R`$. This is because at $`\omega =\omega _R`$ the dressed spin-up level $`\stackrel{~}{\epsilon }_d`$ is still lower than $`\mu _{}`$ and thus the electron in the dot cannot tunnel out via small thermal excitation. Due to the state exchange in QD, as shown in the inset of Fig. 2(a), when $`\omega `$ increases further, crossing the avoided crossing between the two spin levels, the upper eigenstate is dominated by spin component $`|`$ while the lower state by $`|`$. Thus when the upper level is higher than $`\mu _{}`$, spin tunneling processes occurs, and a peak develops in the current spectrum. When the temperature is increased to $`k_BT|\mu _{}\stackrel{~}{\epsilon }_d|`$, the spin-current peak shifts to $`\omega =\omega _R`$ with a larger peak amplitude. Further increasing temperature such that $`k_BT>\omega `$ will smear out the mismatch between $`\mu _{}`$ and $`\mu _{}`$, and thus the spin current amplitude begins to decrease, as shown in Fig. 2(a). The dot levels can be conveniently controlled by the gate voltage $`V_g`$. When $`V_g`$ is modulated such that $`\stackrel{~}{\epsilon }_{d\sigma }`$ are at the interval between $`\mu _{}`$ and $`\mu _{}`$, then the spin current peak will occur exactly at $`\omega =\omega _R`$ \[see Fig. 2(b)\].

The model can also illustrate the spin-current generation in the Coulomb blockade regime when Coulomb interaction $`U`$ is treated in the mean-field formalism. As an example, we show in Fig. 3 the spin current versus gate voltage for several values of $`\omega `$. Two Coulomb peaks with interval $`U`$ can be resolved, and as in Fig. 2, the current amplitude increases when the magnetic field frequency is tuned towards to two-level resonance $`\omega =\omega _R`$.

The electronic correlation effect.—In the low temperature, strong coupling, and large charging energy regime, correlation effects enter and are expected to significantly influence the spin transport. To illustrate correlation effects, we again use the equation of motion method to solve the retarded Green function in (2), which generates higher-order Green functions that have to be truncated to close the equationMeir . After a straightforward calculation, in the infinite-$`U`$ limit we obtain

$$𝐆^r(ϵ)=[ϵ𝐈\widehat{\stackrel{~}{\epsilon }}_d𝚺_0^r(ϵ)𝚺_1^r(ϵ)]^1(𝐈𝐧_c),$$
(10)
where $`(\widehat{\stackrel{~}{\epsilon }}_d)_{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}\stackrel{~}{\epsilon }_{d\sigma }`$ and the interaction self-energies are given by $`[𝚺_1^r(ϵ)]_{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}_k[\frac{\mathrm{sin}^2\varphi |V_k|^2f(ϵ_{k\sigma })}{ϵϵ_{k\sigma }\epsilon _{c\sigma }+\epsilon _{c\overline{\sigma }}}+\frac{\mathrm{cos}^2\varphi |V_k|^2f(ϵ_{k\overline{\sigma }})}{ϵϵ_{k\overline{\sigma }}\epsilon _{c\sigma }+\epsilon _{c\overline{\sigma }}}]`$. $`(𝐧_c)_{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}c_{\overline{\sigma }}^+c_{\overline{\sigma }}`$ needs to be calculated self-consistently via the relation $`c_\sigma ^+c_\sigma =\mathrm{Im}\frac{dϵ}{2\pi }𝐆_{\sigma \sigma }^<(ϵ)`$. $`𝐆^<(ϵ)`$ is difficult to obtain since the lesser self-energy cannot be given exactly in the interacting regime. However, the situation simplifies in the steady-state transportSun2 , in which the occupation on the dot is time-independent, $`i\frac{d}{dt}c_\sigma ^+(t)c_\sigma (t)=0`$. Writing down the Heisenberg equations for $`c_\sigma ^+(t)`$ and $`c_\sigma (t)`$, and integrating out the lead electron operators, we find the self-consistent integral equations for the retarded Green functions as follows:

$$c_{}^+c_{}=\mathrm{Im}\frac{dϵ}{2\pi }𝐆_{}^<(ϵ)=\frac{dϵ}{2\pi }\mathrm{Im}(𝐆_{}^a𝐆_{}^r)[\mathrm{cos}^2\varphi f_{}+\mathrm{sin}^2\varphi f_{}],$$
(11a)
$$c_{}^+c_{}=\mathrm{Im}\frac{dϵ}{2\pi }𝐆_{}^<(ϵ)=\frac{dϵ}{2\pi }\mathrm{Im}(𝐆_{}^a𝐆_{}^r)[\mathrm{sin}^2\varphi f_{}(ϵ)+\mathrm{cos}^2\varphi f_{}(ϵ)],$$
(11b)
which now, together with (10), close the equations of motion for $`𝐆^r`$. Similarly, the quantities $`c_\sigma ^+c_{\overline{\sigma }}=i\frac{dϵ}{2\pi }G_{\sigma \overline{\sigma }}^<(ϵ)`$ can also be expressed by the retarded Green functions. Thus we can calculate the spin current in (5) directly after $`𝐆^r(ϵ)`$ is numerically obtained.

The magnetic-field-induced excitation properties are shown in Fig. 4 by plotting local spin-resolved spectral densities $`\rho _\sigma (ϵ)=[\mathrm{𝐮𝐆}_c^r(ϵ)𝐮^\text{T}]_{\sigma \sigma }`$ for different values of $`\omega `$. In the absence of a rotating magnetic field ($`\omega ,B_1=0`$), as observed from Fig. 4(a), the spectral density for each spin component is characterized by a broad single-particle peak around $`\epsilon _{d\sigma }`$, and a sharp Kondo peak at Zeeman energies $`ϵ=\omega _R`$ for up spin and $`ϵ=\omega _R`$ for down spin. The spectral weight of $`\rho _{}`$ is enhanced by the magnetic field, while $`\rho _{}`$ is greatly suppressed. Therefore, a net spin moment develops in the QD. Physically, the Kondo peak at $`ϵ=\omega _R`$ is due to the cotunneling processes (see inset) in which one spin-down electron in the lead at the Fermi level tunnels into the dot and occupies energy level $`\stackrel{~}{\epsilon }_{}`$, followed by another electron on $`\stackrel{~}{\epsilon }_{}`$ tunneling out to the spin-up channel with the energy $`\omega _R`$ below from Fermi level $`E_F`$. Such cotunneling processes transfer electrons from spin-up channel to spin-down channel of the lead. Thus a spin current may be induced in an efficient way by these Kondo-type cotunneling processes. Similarly, the Kondo peak at $`ϵ=\omega _R`$ in $`\rho _{}(ϵ)`$ is due to a superposition of many cotunneling processes in which one spin down electron with the energy $`\omega _R`$ above the Fermi level tunnels into the dot and simultaneously another spin-up electron in the dot tunnels out to the Fermi level of the lead. Both kinds of cotunneling processes contribute to the spin-current generation. When the rotating magnetic field is switched on, the spin Kondo resonances and spectral structures begin to change. At $`\omega <\omega _R`$, the broad single-spin and sharp Kondo peaks shift toward higher energies for up spin, while they shift to lower energies for down spin \[Fig. 4(b)\]. At two-level resonance $`\omega =\omega _R`$, as shown in Fig. 4(c): (i) each spin Kondo singularity splits into two prominent peaks around $`\mu _{}`$ and $`\mu _{}`$; (ii) $`\rho _{}(ϵ)`$ and $`\rho _{}(ϵ)`$ completely overlap with identical spectral weight. This overlap is caused by the fact that at $`\omega =\omega _R`$, each spin eigenstate in the dot is a strong superposition of $`|`$ and $`|`$ states with equal probability amplitude. When $`\omega `$ increases further, the two-level superposition is suppressed, and the state exchange occurs. As a response, the spectral densities of two spins also exchange their structures; as can be seen from Fig. 4(d), $`\rho _{}(ϵ)`$ is enhanced while $`\rho _{}(ϵ)`$ is greatly suppressed.

To see the influence of the many-body correlation effect on the spin transport, we show in Figure 5 the spin current as a function of rotating frequency. Compared to the non-interacting case, it can be seen that the spin current is greatly enhanced by the cotunneling processes in the strong coupling, low temperature regime. A detailed study of cotunneling effects on the spin current is in progress.

In summary, we have analyzed the spin current properties of an interacting QD exposed to a rotating magnetic field. The spin-flip process in the QD and the effective difference of two spin chemical potentials suffices the generation of spin current with no charge current. The interplay of the many-body cotunneling process and coherent two-level resonance in the QD produces a remarkable enhancement of spin current. Since a controllable generation of spin current provides an efficient source for spin injection, spin detection, and virous variety of spin-based devices, we expect the present results may have practical applications in the field of spintronics.

This work was supported by grants CNSF 90103024, MOST-G001CB3095, and U.S.-DOE-DE-FG03-98ER45687.

Figure captions

Fig. 1. Model setup for spin transport through a QD. The degenerate spin levels in the dot are split and coupled by the external magnetic fields $`B_0`$ and $`B_1`$, respectively. Due to the imbalance between spin-up and spin-down chemical potentials (see text), the electron tunnels from the spin-down channel, through the QD, to the spin-up channel.

Fig. 2. (a) Non-interacting spin current as a function of rotating frequency at different temperatures with the undressed spin levels set $`\epsilon _d=2.5`$, $`\epsilon _d=1.5`$; (b) Contour plot of spin current as a function of gate voltage and rotating frequency. The spin levels in the dot are modulated by the gate voltage as $`\epsilon _{d\sigma }(V_g)=\epsilon _{d\sigma }+eV_g`$ with $`\epsilon _d=1`$, $`\epsilon _d=2`$. Other parameters are $`\mathrm{\Gamma }=0.1`$ and $`g\mu _BB_1=0.2`$.

Fig. 3. Spin current versus gate voltage for different rotating frequencies with $`U=2`$. Other parameters are the same as in Fig. 2(b).

Fig. 4. The spin up (solid line) and spin down (dashed line) spectral densities at different rotating frequency. Parameters are $`\epsilon _d=2.5\mathrm{\Gamma }`$, $`\epsilon _d=1.5\mathrm{\Gamma }`$, $`T=0.001\mathrm{\Gamma }`$, and $`g\mu _BB_1=0.2\mathrm{\Gamma }`$.

Fig. 5. Spin current as a function of rotating frequency in the Kondo-type cotunneling regime (solid line). The non-interacting spin current is also shown (dotted line) for comparison. Other parameters are the same as in Fig. 4.