text:
# Nonlinear Hydromagnetic Wave Support of a Stratified Molecular Cloud

## 1. Introduction

Interstellar molecular clouds, the sites of current star formation in our Galaxy, have long been known to yield supersonic line widths of molecular spectral lines (e.g., see Zuckerman & Palmer 1974). Objects classified as molecular clouds span a large range of mean radii ($`R1\mathrm{pc}100\mathrm{pc}`$), masses ($`M10^2M_{}10^6M_{}`$) and mean number density ($`n10^1\mathrm{cm}^310^3\mathrm{cm}^3`$). In fact, these quantities are correlated with the one-dimensional velocity dispersion $`\sigma `$ through the well known line-width-size-density relations (e.g., Solomon et al. 1987):

$`\sigma `$ $`=`$ $`0.72(R/\mathrm{pc})^{0.5}\mathrm{km}\mathrm{s}^1,`$ (1)
$`n`$ $`=`$ $`2.3\times \mathrm{\hspace{0.17em}10}^3(R/\mathrm{pc})^1\mathrm{cm}^3.`$ (2)
Thus, the velocity dispersion is typically supersonic since the sound speed $`c_\mathrm{s}`$ is only $`0.2\mathrm{km}\mathrm{s}^1`$ for the typical molecular cloud temperature $`T10\mathrm{K}`$ (e.g., Goldsmith & Langer 1978).

The largest clouds, of mass $`M10^4M_{}`$, often classified as Giant Molecular Clouds (GMC’s) are in fact complexes of smaller clouds, since the volume averaged density may be lower than the excitation density of CO (e.g., Blitz & Williams 1999), and also lower than allowed from thermal stability arguments (Falgarone & Puget 1986). Hence, the basic building blocks of interstellar molecular clouds, which contain most of mass of molecular material, are the dark (or dwarf) molecular clouds, which have $`R1\mathrm{pc}10\mathrm{pc}`$, $`M10^2M_{}10^4M_{}`$ and $`n10^2\mathrm{cm}^310^3\mathrm{cm}^3`$. These clouds have velocity dispersion $`\sigma 12\mathrm{km}\mathrm{s}^1`$, and represent the class of objects that we are interested in modeling in this study. We also note that observed smaller scale ($`R0.1`$ pc) dense cores are a separate object class which are embedded within dark clouds and collectively contain only a small fraction of the total cloud mass.

While even dark clouds have masses that significantly exceed the thermal Jeans (1928) mass

$`M_\mathrm{J}`$ $`=`$ $`\left({\displaystyle \frac{\pi kT}{mG}}\right)^{3/2}{\displaystyle \frac{1}{\rho ^{1/2}}}`$ (3)
$`=`$ $`17\left({\displaystyle \frac{T}{10\mathrm{K}}}\right)^{3/2}\left({\displaystyle \frac{n}{10^3\mathrm{cm}^3}}\right)^{1/2}M_{}`$
(where we have used $`\rho =mn`$, and $`m=2.33m_\mathrm{H}`$, in which $`m_\mathrm{H}`$ is the mass of a hydrogen atom), the line-width-size-density relations do imply that molecular clouds are individually in an approximate virial balance between turbulent and gravitational energies. In this paper, we equate the presence of nonthermal line widths with the presence of a random superposition of nonlinear (presumably hydromagnetic) waves, and refer loosely to the latter as “turbulence”. Within each cloud, the turbulence is expected to collectively exert a force (e.g., Chandrasekhar 1951) which resists the inward pull of gravity.

The origin and persistence of turbulent motions in molecular clouds remain an active area of investigation. It was long understood (see Mestel 1965; Goldreich & Kwan 1974) that supersonic hydrodynamic motions would decay rapidly through shocks, thereby creating a mystery of why the turbulence was commonly observed throughout the lifetime of molecular clouds. For GMC’s, the lifetime is estimated to be a few $`\times \mathrm{\hspace{0.17em}10}^7`$ yr (e.g., Blitz & Williams 1999), which is at least a few times longer than the crossing time $`t_\mathrm{c}2R/\sigma `$ of the complexes. The smaller dark molecular clouds may have even longer lifetimes $`10^8`$ yr (Shu, Adams. & Lizano 1987). Since slow and fast mode MHD waves are also compressive and can be highly dissipative, it was suggested by Arons & Max (1975) that the transverse Alfvén mode might be a long-lived component, thereby preventing a rapid overall collapse of the clouds and explaining the observed persistence of the turbulence over their inferred lifetime. Mouschovias (1975) made the related suggestion that the long-lived component may be due to standing waves, i.e., normal mode oscillations, of magnetized clouds. Such global magneto-gravitational motions cannot be studied with a simple plane wave analysis in an infinite uniform medium.

The magnetohydrodynamic (MHD) wave picture has been strengthened by the detection of large-scale magnetic fields within molecular clouds, through maps of polarized absorption and emission (e.g., Vrba, Strom, & Strom 1976; Goodman et al. 1990; Schleuning 1998; Matthews & Wilson 2002), and Zeeman measurements of the line-of-sight magnetic field strengths (Crutcher 1999 and references therein). The latter imply that the magnetic energy (like the turbulent energy) is comparable to the gravitational potential energy; equivalently, the mass-to-flux ratio is close a critical value

$$\left(\frac{M}{\mathrm{\Phi }}\right)_{\mathrm{crit}}=c_1G^{1/2}.$$
(4)
In the above equation, the constant $`c_1`$ has been calculated to be in the range 0.13-0.17 based on detailed two-dimensional equilibrium states (Mouschovias & Spitzer 1976; Tomisaka, Ikeuchi, & Nakamura 1988) and is equal to $`1/(2\pi )`$ in the case of a flattened one-dimensional layer (Nakano & Nakamura 1978). Collectively, the magnetic field measurements allow the possibility that clouds are supported lateral to the large-scale field $`B`$ by its associated Lorentz force, while these magnetic field lines act as a carrier of Alfvénic disturbances with $`\delta BB`$, explaining the observed spectral line widths and preventing the clouds from assuming a very flattened configuration,

The dynamical effect of propagating MHD waves using analytic or semi-analytic means has been studied by several authors. Dewar (1970) developed a formalism for calculating the effect of small amplitude hydromagnetic waves on a slowly varying background medium, using the WKB approach. Assuming ideal MHD and no dissipation of the waves, this leads to a steady-state relation between wave pressure $`P_\mathrm{w}`$ and gas density $`\rho `$ of the form $`P_\mathrm{w}\rho ^{1/2}`$ (McKee & Zweibel 1995). Simplified calculations of the effect of small amplitude MHD waves on a molecular cloud have been performed by Fatuzzo & Adams (1993) and Martin, Heyvaerts, & Priest (1997). Both satisfy the above scaling of $`P_\mathrm{w}`$ in the ideal MHD limit. In particular, the WKB model of Martin et al. (1997) yields a steady-state density structure of an infinite one-dimensional cloud supported by short-wavelength Alfvén waves in the ideal MHD limit, and also when accounting for damping of linear waves by ion-neutral friction. In addition to ion-neutral friction, which damps even linear waves (Kulsrud & Pearce 1969), there are several nonlinear effects which will work to enhance dissipation. The second order effect of a gradient in the magnetic pressure $`\delta B^2/8\pi `$ will in general lead to steepening of the waves followed by dissipation (Cohen & Kulsrud 1974). Zweibel & Josafatsson (1983) state the form of this dissipation rate, which can dominate the process of ion-neutral friction for nonlinear and/or long wavelength modes. There are also other known nonlinear avenues for wave dissipation, such as the conversion of a parallel-propagating Alfvén wave into an acoustic wave and another Alfvén wave traveling in the opposite direction (Sagdeev & Galeev 1969).

Recently, several studies have resulted in a solution of the full set of nonlinear equations of ideal MHD using finite difference approximations. Gammie & Ostriker (1996) performed a one-dimensional numerical simulation of MHD turbulence in a periodic domain. This has been followed up by several multi-dimensional simulations, also in a periodic domain (e.g., Stone, Ostriker, & Gammie 1998; Ostriker, Gammie & Stone 1999; Mac Low et al. 1998; Mac Low 1999; Padoan & Nordlund 1999; Ostriker, Stone, & Gammie 2001). These models impose prescribed velocity fluctuations, at the initial time and sometimes also throughout the computed time, and investigate the dissipation rate of turbulence, as well as various properties of turbulent fluctuations. One of the important results of these papers is that the decay time of the MHD turbulence is comparable to the crossing time over the driving scale of the turbulence. The nonlinear coupling of the Alfvén to the fast and slow mode MHD waves is considered to be a significant source of the dissipation which reduces the lifetime of MHD turbulence relative to the ideal of pure Alfvénic turbulence. We note that periodic models represent only a local region of a much larger molecular cloud, and maintain a fixed mean density. They cannot study global effects associated with the density stratification in a cloud or the existence of a cloud boundary.

In this paper, we perform a different type of numerical simulation of MHD turbulence. We concentrate on one self-gravitating cloud and study the effect of the turbulence on the mechanical structure of the cloud. This corresponds to an extension of the model studied by Martin et al. (1997) into a fully nonlinear counterpart. As an initial condition, we use a hydrostatic equilibrium between thermal pressure and self-gravity in a cloud that is bounded by an external high temperature medium. We input turbulent energy into the system continuously and see how the mechanical equilibrium changes. Similar numerical simulations have been performed to study the propagation of Alfvén waves in the solar chromosphere and corona (e.g., Hollweg, Jackson & Galloway 1982; Mariska & Hollweg 1985; Hollweg 1992; Kudoh & Shibata 1999; Saito, Kudoh, & Shibata 2001). We extend this class of model to a self-gravitating cloud, and study the relation between the strength of the turbulence and various global properties of a molecular cloud. This is the first of a series of papers on global models of MHD wave support in molecular clouds. In this paper, we develop high-resolution one-dimensional models under the assumption of ideal MHD.

The paper is organized as follows. The numerical model we used for the simulation is summarized in § 2. The results of the simulation are described in § 3. We add some discussion of the results in § 4 and summarize the paper in § 5.

## 2. The Numerical Model

### 2.1. Schematic View of our Model

Figure 1 shows a schematic picture of our model. We consider a molecular cloud that is threaded by a large-scale magnetic field, and concentrate on a local region of the molecular cloud enclosed by the rectangle in the figure. We assume a driving force near the midplane of the cloud and follow the dynamical evolution of the vertical structure of the cloud.

Fig. 1.  Schematic picture of our model. A molecular cloud that is threaded by a large-scale magnetic field is considered. A local region of the molecular cloud enclosed by the rectangle in the figure can represent the region we are modeling in this paper. We input a driving force near the midplane of the cloud and follow the dynamical evolution of the vertical structure of the cloud.

### 2.2. Assumptions

For simplicity, we assume ideal MHD, 1.5-dimensions, and isothermality for each Lagrangian fluid element. The 1.5-dimensional approximation means that physical quantities depend on only one coordinate, but we evolve nonzero components of vectors in one additional direction. Isothermality for each Lagrangian fluid element means that the temperature does not change in time for each fluid element as it moves through Eulerian space, which is different from the assumption of a uniform time-independent temperature throughout the region. These assumptions appear in § 2.3 more concretely.

### 2.3. Basic Equations

We use local Cartesian coordinates $`(x,y,z)`$ on the molecular cloud, where we set $`z`$ to be the direction of the large-scale magnetic field. According to the symmetry of the 1.5-dimensional approximation, we set

$$\frac{}{x}=\frac{}{y}=0.$$
(5)
The above symmetry and the divergence-free condition on the magnetic field imply

$$B_z=\mathrm{constant},$$
(6)
where $`B_z`$ is the $`z`$-component of the magnetic field that threads the molecular cloud. Moreover, from the assumption of linear polarization of the waves, we can set

$$v_x=B_x=0$$
(7)
without loss of generality, where $`v_x`$ and $`B_x`$ are the $`x`$-components of the velocity and magnetic field, respectively.

Therefore, the basic MHD equations we use in this paper are as follows: mass conservation,

$$\frac{\rho }{t}+v_z\frac{\rho }{z}=\rho \frac{v_z}{z};$$
(8)
the $`z`$-component of the momentum equation,

$$\frac{v_z}{t}+v_z\frac{v_z}{z}=\frac{1}{\rho }\frac{P}{z}\frac{1}{4\pi \rho }B_y\frac{B_y}{z}+g_z;$$
(9)
the $`y`$-component of the momentum equation,

$$\frac{v_y}{t}+v_z\frac{v_y}{z}=\frac{1}{4\pi \rho }B_z\frac{B_y}{z};$$
(10)
the equation of energy,

$$\frac{T}{t}+v_z\frac{T}{z}=(\gamma 1)T\frac{v_z}{z};$$
(11)
the $`y`$-component of the induction equation,

$$\frac{B_y}{t}=\frac{}{z}(v_zB_y+v_yB_z);$$
(12)
the equation of state,

$$P=\rho \frac{kT}{m}=c_s^2\rho ,$$
(13)
where

$$c_s\sqrt{\frac{kT}{m}}$$
(14)
is the isothermal sound speed; and the Poisson equation,

$$\frac{g_z}{z}=4\pi G\rho .$$
(15)
In these equations, $`t`$ is the time, $`G`$ is the gravitational constant, $`k`$ is Boltzmann’s constant, $`m`$ is the mean molecular mass, $`\gamma `$ is the specific heat ratio, $`\rho `$ is the density, $`P`$ is the pressure, $`T`$ is the temperature, $`g_z`$ is the $`z`$-component of the gravitational field, $`v_z`$ is the $`z`$-component of velocity, $`v_y`$ is the $`y`$-component of velocity, and $`B_y`$ is the $`y`$-component of the magnetic field.

In equation (11), we assume $`\gamma =1`$, so that the energy equation becomes

$$\frac{T}{t}+v_z\frac{T}{z}=0,$$
(16)
which quantifies the assumption of isothermality for each Lagrangian fluid element.

### 2.4. Initial Conditions

As an initial condition, we assume hydrostatic equilibrium of a self-gravitating one-dimensional cloud. The hydrostatic equilibrium is calculated from the equations

$$\frac{1}{\rho }\frac{dP}{dz}=g_z,$$
(17)
$$\frac{dg_z}{dz}=4\pi G\rho ,$$
(18)
and

$$P=\rho \frac{kT}{m},$$
(19)
subject to the boundary conditions

$$g_z(z=0)=0,$$
(20)
$$\rho (z=0)=\rho _0,$$
(21)
and

$$P(z=0)=\rho _0\frac{kT_0}{m},$$
(22)
where $`\rho _0`$ and $`T_0`$ are the initial density and initial temperature at $`z=0`$, respectively.

In order to solve the above equations, we need to assume an initial temperature distribution. If the temperature is uniform throughout the region, we have the following analytic solution $`\rho _\mathrm{S}`$ found by Spitzer (1942);

$$\rho _\mathrm{S}(z)=\rho _0\mathrm{sech}^2(z/H_0),$$
(23)
where

$$H_0=\frac{c_{s0}}{\sqrt{2\pi G\rho _0}}$$
(24)
is the scale height, and

$$c_{s0}=\sqrt{\frac{kT_0}{m}}.$$
(25)
However, an isothermal molecular cloud is usually surrounded by warm or hot material, such as neutral hydrogen or ionized gas. Therefore, we assume the initial temperature distribution to be

$$T(z)=T_0+\frac{1}{2}(T_cT_0)\left[1+\mathrm{tanh}\left(\frac{|z|z_c}{z_d}\right)\right],$$
(26)
where we take $`T_c=100T_0`$, $`z_c=3H_0`$, and $`z_d=0.2H_0`$ throughout the paper. This distribution shows that the temperature is uniform and equal to $`T_0`$ in the region of $`0z<z_c=3H_0`$ and smoothly increases to another uniform value $`T_c=100T_0`$ at $`zz_c=3H_0`$. By using this temperature distribution, we can solve the ordinary differential equations (17)-(19) numerically. The numerical solution of these equations shows that the initial distribution of density is almost the same as Spitzer’s solution in the region of $`0z<z_c`$ (see Fig. 2).

We also assume the following initial conditions;

$$v_z(z)=v_y(z)=0,$$
(27)
$$B_y(z)=0,$$
(28)
$$B_z(z)=B_0,$$
(29)
where $`B_0`$ is a constant. According to equation (6), $`B_z`$ is spatially uniform and independent of time throughout the calculations.

### 2.5. Driving Force

We introduce a perturbation into the initially hydrostatic cloud by adding a driving force, $`F(z,t)`$, into the $`y`$-component of the momentum equation (10) as follows:

$$\rho (\frac{v_y}{t}+v_z\frac{v_y}{z})=\frac{1}{4\pi }B_z\frac{B_y}{z}+F(z,t),$$
(30)
where

$$F(z,t)=\{\begin{array}{cc}\rho a_d(\frac{t}{10t_0})\mathrm{sin}(2\pi \nu _0t)\mathrm{exp}[(\frac{z}{z_a})^2]\hfill & (t<10t_0)\hfill \\ \rho a_d\mathrm{sin}(2\pi \nu _0t)\mathrm{exp}[(\frac{z}{z_a})^2]\hfill & (10t_0t40t_0)\hfill \\ 0\hfill & (t>40t_0),\hfill \end{array}$$
(31)
and

$$t_0=H_0/c_{s0}.$$
(32)
Furthermore, $`a_d`$ is the amplitude of the induced acceleration, $`\nu _0`$ is the frequency of the driving force, and $`z_a`$ represents the region in which we input the driving force. The equations show that we input the sinusoidal driving force near the midplane of the cloud, and increase the maximum driving force linearly with time until $`t=10t_0`$, and maintain it to be constant during $`10t_0t40t_0`$. After $`t=40t_0`$, we terminate the driving force.

### 2.6. Boundary Conditions

We used a mirror symmetric boundary condition at $`z=0`$ and a free boundary at $`z=z_{\mathrm{out}}`$, the outer boundary of the calculation. In order to remove the reflection of waves at the outer boundary, we set $`z_{\mathrm{out}}`$ to be a large value, i.e.,

$$z_{\mathrm{out}}2.6\times 10^4H_0,$$
(33)
and use nonuniform grid spacing for large $`z`$ (see § 2.8).

In addition to this, we set the gravitational field to be zero at large distances by introducing an artificial acceleration ($`g_{}`$) into equations (9) and (17), i.e.,

$$g_{}(z)=\frac{g_z(z)}{2}\left[1+\mathrm{tanh}\left(\frac{zz_g}{5H_0}\right)\right],$$
(34)
where we take $`z_g=100H_0`$ in this paper. This force is zero until $`zz_g=100H_0`$, but it becomes essentially equal to $`g_z`$ with negative sign for $`z>z_g=100H_0`$ and compensates the gravity there.

The above boundary conditions show that we have effectively two outer boundaries. The first one is $`z_{\mathrm{out}}`$ and the second one is $`z_g`$, which is the boundary for the gravitational field. In order to remove the reflection of waves at the outer boundary, it is useful to have a large $`z_{\mathrm{out}}`$. However, a very large $`z_{\mathrm{out}}`$ introduces numerical problems because the density decreases exponentially at large $`z`$, according to the stratification of a self-gravitating cloud. A very low density leads to a very large Alfvén speed, making simulations inefficient by forcing very small time steps for an explicit calculation. By setting the net gravitational field to be zero at large distances, we can avoid an extremely low density for large $`z`$ and therefore take a large $`z_{\mathrm{out}}`$ in order to remove the reflection of Alfvén waves at the outer boundary. We also note that for any slab of finite extent, as opposed to the infinite slab implied by our one-dimensional model, the net gravitational field should indeed decrease at large distances, though not in the specific manner prescribed here.

Due to the added artificial acceleration, our effective numerical boundary is located at $`zz_g=100H_0`$, and we cannot trust the results beyond this region. However, as we show later, all dynamical events we are interested in occur within $`z<50H_0`$, where most of the mass and energy are concentrated. Hence, the effect of the artificial acceleration is negligible for the main results in this simulation.

### 2.7. Numerical Parameters

A natural set of fundamental units for this problem are $`c_{s0}`$, $`H_0`$, and $`\rho _0`$. These yield a time unit $`t_0=H_0/c_{s0}`$. The initial magnetic field strength introduces one dimensionless free parameter, i.e.

$$\beta _0\frac{8\pi P_0}{B_0^2}=\frac{8\pi \rho _0c_{s0}^2}{B_0^2},$$
(35)
which is the initial ratio of gas to magnetic pressure at $`z=0`$.

In this cloud, $`\beta _0`$ is related to the mass-to-flux ratio. For Spitzer’s self-gravitating cloud, the mass-to-flux ratio normalized to the critical value is

$$\mu _\mathrm{S}2\pi G^{1/2}\frac{\mathrm{\Sigma }_\mathrm{S}}{B_0},$$
(36)
where

$$\mathrm{\Sigma }_\mathrm{S}=_{\mathrm{}}^{\mathrm{}}\rho _S𝑑z=2\rho _0H_0$$
(37)
is the column density of the Spitzer’s self-gravitating cloud. Therefore,

$$\beta _0=\mu _\mathrm{S}^2.$$
(38)
The column density of the cloud we used in this paper is almost equal to that of Spitzer’s cloud. If we define the column density of the cloud, $`\mathrm{\Sigma }`$, as the integral of density within $`z_c<z<z_c`$, then

$$\mathrm{\Sigma }=_{z_c}^{z_c}\rho (t=0)𝑑z0.988\mathrm{\Sigma }_\mathrm{S}.$$
(39)
This means that we can use the value of $`\mu _\mathrm{S}`$ an excellent approximation to the dimensionless mass-to-flux ratio of the model cloud.

Observations show that the dimensionless mass-to-flux ratios of molecular clouds are close to unity over a range of scales (Crutcher 1999; Shu et al. 1999). Hence, we take $`\beta _0=1`$ in the models presented in this paper.

The driving force introduces three more important free parameters: $`\stackrel{~}{a}_d=a_d(H_0/c_{s0}^2)`$, the dimensionless amplitude of the acceleration due to driving, $`\stackrel{~}{\nu }_0=\nu _0t_0`$, the dimensionless frequency of driving, and $`\stackrel{~}{z}_a=z_a/H_0`$, the dimensionless scale of the driving region. For simplicity, we take $`\stackrel{~}{\nu }_0=1`$ and $`\stackrel{~}{z}_a=0.1`$ throughout this paper, and adjust the strength of driving by varying $`\stackrel{~}{a}_d`$ from $`10`$ to $`50`$.

Dimensional values of all quantities can be found through a choice of $`T_0`$ and $`\rho _0`$, along with the values of the dimensionless free parameters. For example, if $`T_0=10\mathrm{K}`$ and $`n_0=\rho _0/m=10^4\mathrm{cm}^3`$, then $`c_{s0}=0.2\mathrm{km}\mathrm{s}^1`$, $`H_0=0.05\mathrm{pc}`$, $`N_S=\mathrm{\Sigma }_S/m=3\times 10^{21}\mathrm{cm}^2`$, $`t_0=2.5\times 10^5`$ yr, and $`B_0=20\mu \mathrm{G}`$ if $`\beta _0=1`$.

### 2.8. Numerical Technique

In order to solve the equations numerically, we use the CIP method (Yabe & Aoki 1991) for equations (8)-(10) and (16), and the MOCCT method (Stone & Norman 1992) for equation (12). The combination of the CIP and MOCCT methods is summarized in Kudoh, Matsumoto & Shibata (1999). The CIP method is a useful method to solve an advection equation such as equation (16) accurately, and is also applicable to advection terms of equations (8)-(10).

Due to the mirror-symmetric boundary condition at $`z=0`$, the Poisson equation (15) can be simply integrated from the midplane of the cloud;

$$g_z(z)=4\pi G_0^z\rho (z)𝑑z.$$
(40)
We solve this equation by numerical integration.

In this simulation, we actually use variants of the original CIP method. We use a conservative-CIP method for the mass conservation equation (8), which was recently developed by Xiao et al. (2002). This scheme assures exact conservation of mass and less numerical oscillations. The original CIP method does not assure the exact conservation of mass, so that a systematic deviation during each time step can cause a problem for long time integrations; this is especially a problem if the Poisson equation, which utilizes the mass distribution, is being solved simultaneously. Hence, the conservative-CIP method significantly improves the accuracy of our solution. We also use the monotonic-CIP method (Xiao, Yabe & Ito 1996) for the advection terms of equations (9), (10) and (16), and use the CCUP method (Yabe & Wang 1991) for the calculation of gas pressure, in order to get more numerically stable results. The recent developments of CIP related schemes are summarized in Yabe, Xiao & Utsumi (2001).

We used a uniform grid size, $`\mathrm{\Delta }z_i=0.02H_0`$, from $`z=0`$ to $`z=80H_0`$, where $`\mathrm{\Delta }z_i`$ is the grid size of the $`i`$-th grid point. For $`80H_0<z<120H_0`$, we gradually increase the grid size as $`\mathrm{\Delta }z_{i+1}=Min(1.05\mathrm{\Delta }z_i,0.2H_0)`$, where $`Min`$ is a function that picks the minimum from the two values; hence, the maximum grid size in this region is $`0.2H_0`$. For $`120H_0<z<z_{\mathrm{out}}`$, where the gravitational field is compensated by the artificial force, we further increase the grid size according to $`\mathrm{\Delta }z_{i+1}=1.05\mathrm{\Delta }z_i`$. In this paper, we used a total of 4310 grid points, most of which are concentrated in the region of the uniform grid, i.e., $`z80H_0`$.

## 3. Results

### 3.1. Typical Result

In this subsection, we show the results of $`\stackrel{~}{a}_d=30`$ as the typical case.

#### 3.1.1 Density Structure in the Cloud

Figure 2 shows the density and temperature as a function of $`z`$. The dashed lines show the initial distributions and the solid lines show the distribution at $`t=30t_0`$.

Fig. 2.  Density (upper panel) and temperature (lower panel) as a function of $`z`$. The dashed lines show the initial distributions and the solid lines show the distribution at $`t=30t_0`$. The dotted line in the density plot shows the distribution of density in Spitzer’s infinite self-gravitating equilibrium with a uniform temperature $`T_0`$.

The dotted line in the density shows the distribution of Spitzer’s self-gravitating cloud with uniform temperature of $`T_0`$. The initial density deviates from Spitzer’s solution around $`z=z_c=3H_0`$ where the initial temperature increases up to $`100T_0`$. The density decreases rapidly around $`z=z_c=3H_0`$ because of the pressure balance between low temperature cloud and high temperature external medium. For $`z>z_c=3H_0`$, the scale height of the density is very large and the density decreases gradually.

The snapshot of density at $`t=30t_0`$ shows that the density has a complicated structure including many shock fronts. This is caused by the driving force in equation (30). The driving force generates nonlinear Alfvén waves in the cloud which produce a magnetic pressure gradient. The magnetic pressure gradient and thermal pressure gradient usually push the cloud upward, but the self-gravity of the cloud always pulls it down. These up and down motions create the complicated structure in the cloud. On the other hand, the temperature shows a smooth structure. This is due to isothermality for each Lagrangian fluid element. Only the position of the temperature transition region changes in time.

Figure 3 shows the time evolution of the density. The density plots at various times are stacked with time increasing upward in uniform increments of $`0.2t_0`$. Because the driving force increases linearly with time up to $`t=10t_0`$, the density changes gradually at first. After $`t=10t_0`$, the density structure shows many shock waves propagating in the cloud, and significant upward and downward motions of the outer portion of the cloud, including the transition region. After terminating the driving force at $`t=40t_0`$, the shock waves are dissipated in the cloud and the transition region falls to around the initial position, although it is still oscillating.

Fig. 3.  Time evolution of the density. The density plots at various times are stacked with time increasing upward in uniform increments of $`0.2t_0`$.

Figure 4a shows the column density as a function of the time averaged position of several Lagrangian fluid elements, which are equally spaced at time $`t=0`$ with spacing $`\mathrm{\Delta }z=0.1H_0`$ starting at $`z=0.01H_0`$. The time average is calculated between $`t=10t_0`$ and $`t=40t_0`$, while the driving force is input with constant amplitude. The dashed line shows the initial distribution. The Lagrangian fluid elements have constant enclosed column density as a function of time. By using this property, we evaluate the location of Lagrangian fluid elements from the surface density. The difference between the initial distribution and that of the time average shows that the cloud is lifted up. Figure 4b shows the time average of the density for each element as a function of the time averaged position of each element. These time averages were also calculated between $`t=10t_0`$ and $`t=40t_0`$. In contrast to the snapshot in Figure 2, the time averaged density structure shows a smooth distribution. The dashed line shows the initial density distribution of the cloud. The scale height of the time averaged distribution is about 3 times larger than the initial value.

Fig. 4.  (a) Column density as a function of the time averaged position of each Lagrangian fluid element. The dashed line shows the initial distribution. The time average is calculated between $`t=10t_0`$ and $`t=40t_0`$. The Lagrangian fluid elements have constant enclosed column density as a function of time. By using this property, we evaluate the location of Lagrangian fluid elements from the surface density. Each Lagrangian fluid element is equally spaced at time $`t=0`$ with spacing $`\mathrm{\Delta }z=0.1H_0`$ starting at $`z=0.01H_0`$. (b) Time average of the density for each fluid element as a function of the time averaged position of the element. The dashed line shows the initial distribution.

#### 3.1.2 Velocities in the Cloud

Figure 5a shows the time evolution of $`v_y`$ at $`z=0`$. Initially, it is oscillating about zero sinusoidally with the frequency of the driving force. However, as time goes on, the oscillation shows sharp structures, and it becomes nonsymmetric around the mean. The sharp structures are caused by nonlinear effects in the cloud such as shock waves. Moreover, the mean value shows a deviation from zero as time goes on. This implies that the net $`y`$-momentum of the cloud is not zero and the cloud has a mean motion in the $`y`$-direction as well as an oscillatory motion.

Fig. 5.  Transverse velocities versus time. (a) The $`y`$-component of the velocity at $`z=0`$ as a function of time. (b) The mean $`y`$-component of velocity of the cloud as a function of time. (c) The time average of the mean $`y`$-component of velocity for each cycle of the sinusoidal period of the driving force as a function of time. (d) The $`y`$-component of the velocity at $`z=0`$ minus the time average of the mean $`y`$-component of velocity, as a function of time.

This mean motion is ultimately caused by the driving force. Prior to $`t=10t_0`$, the driving force is not symmetric in time, because it is proportional to $`t\mathrm{sin}(2\pi \nu _0t)`$. However, even if the driving force is symmetric in time, the net momentum in the cloud never remains exactly zero if we use a driving force like the one in equation (31) which is confined to a region near the midplane. The net momentum of the cloud is generated by the nonlinear and nonsymmetric propagation of this disturbance into the stratified cloud, the nonsymmetric restoring forces in the system, and nonsymmetric reflection at the cloud boundary. As disturbances propagate, the restoring force due to the magnetic field is not exactly symmetric in space and time. Also, some of the wave momentum is reflected at the boundary between the cold cloud and hot external medium, and part of it escapes from the cloud. Hence, a net $`y`$-component of momentum can and does appear within the cloud, creating a mean motion. We note that if we were modeling a two-dimensional cloud, i.e., the system had finite extent in the $`y`$-direction, we could reduce the net drift with a $`y`$-distribution of driving force such that the net driving force is zero in the $`y`$-direction. Also, if the system was closed in the $`z`$-direction, as in a periodic boundary system, it is possible to add momentum in a regulated way so that the net transverse momentum remains exactly zero (e.g., Gammie & Ostriker 1996). However, the combined effects of our 1.5-dimensional approximation, cloud stratification, and open cloud boundary make it difficult to regulate the net drift of the cloud in the $`y`$-direction.

Therefore, for our analysis, we divide $`v_y`$ into two parts. The first is the mean velocity which shows the mean motion of the entire cloud. The second is the oscillating component of the velocity. We calculate the mean $`y`$-component of velocity of the cloud as

$$v_\mathrm{m}(t)=\frac{_0^{z_f(t)}\rho v_y𝑑z}{_0^{z_f(t)}\rho 𝑑z}.$$
(41)
In this equation, $`z_f(t)`$ is the full mass position of the cloud, which is defined by

$$_0^{z_f(t)}\rho 𝑑z=0.998\frac{\mathrm{\Sigma }_\mathrm{S}}{2},$$
(42)
and $`z_f(t)`$ corresponds to the position of the Lagrangian fluid element which is initially located at $`z_c`$, the initial position of the transition region of the temperature. The time evolution of $`v_m`$ is shown in Figure 5b.

Figure 5b shows that the mean velocity is still oscillating while the driving force is input ($`t=040t_0`$). In order to remove the oscillation, we take time average of $`v_m`$ for each cycle of the sinusoidal period ($`t_0`$) of the driving force and define $`v_m`$ in the following manner. For example, $`v_m`$ between $`nt_0`$ and $`(n+1)t_0`$ is calculated as

$$v_m(t=nt_0(n+1)t_0)=\frac{1}{t_0}_{nt_0}^{(n+1)t_0}v_{\mathrm{mean}}(t^{})𝑑t^{},$$
(43)
where $`n`$ is an integer. This calculation is done from $`n=0`$ to $`n=59`$. Figure 5c shows $`v_m`$ as a function of time. According to the definition, $`v_m`$ has the same value between $`nt_0`$ and $`(n+1)t_0`$.

Finally, we define the oscillating component of the $`y`$-velocity as

$$v_y^{}=v_yv_m.$$
(44)
Figure 5d shows the time evolution of $`v_y^{}`$. In contrast to Figure 5a, it is oscillating around $`v_y^{}=0`$.

Figure 6 shows the time average of various velocities for each Lagrangian fluid element plotted versus the time average of the position of each element. Open circles show the time averaged root mean square of $`v_z`$, filled circles show that of $`v_y^{}`$, and the dashed line shows the time average of the sound speed. The solid line shows the root mean square of the time average of the effective one-dimensional velocity dispersion

$$\sigma =\sqrt{\frac{1}{2}[v_z^2+(v_y^{})^2]+c_s^2}$$
(45)
for each fluid element. The time average is taken between $`t=10t_0`$ and $`t=40t_0`$. This figure shows that the $`y`$-component of the velocity is the dominant component in the cloud. The averaged $`y`$-component of the velocity increases as a function of $`z_t`$, except for near the midplane. This means that the largest velocity dispersion occurs in the low density region. This is a similar tendency to that of linear Alfvén waves. If we assume the WKB approximation and no wave dissipation, the energy flux of the waves is constant, i.e.,

$$\rho (v_y^{})^2V_A=\mathrm{constant},$$
(46)
where $`V_A`$ is the Alfvén velocity of the background magnetic field, i.e.,

$$V_A=\frac{B_z}{(4\pi \rho )^{1/2}}.$$
(47)
This leads to

$$v_y^{}\rho ^{1/4}.$$
(48)
Although this relation is not exactly applicable to our nonlinear result, our result shows a similar tendency. We discuss this further in the next section.

Fig. 6.  Time averaged root mean square of various velocities for each Lagrangian fluid element versus the time average of the position of each element. Open circles show the time averaged root mean square of $`v_z`$, filled circles show that of $`v_y^{}`$, and the dashed line shows the time average of the sound speed. The solid line shows the root mean square of the time average of the effective one-dimensional velocity dispersion.

#### 3.1.3 Pressures in the Cloud

Figure 7a shows the time averages of thermal pressure, magnetic pressure, and dynamic pressures for each Lagrangian fluid element as a function of the time average of the position of each element. The thick dashed line shows the thermal pressure, the thick solid line shows the magnetic pressure of the $`y`$-component of the magnetic field, i.e., $`\frac{B_y^2}{8\pi }`$, the dash-dotted line shows the dynamic pressure for the $`y`$-component of velocity, i.e., $`\frac{1}{2}\rho (v_y^{})^2`$, and the dotted line shows the dynamic pressure of the $`z`$-component of velocity, i.e., $`\frac{1}{2}\rho v_z^2`$. The thin straight line shows the magnetic pressure of background magnetic field, $`\frac{B_0^2}{8\pi }`$, and the thin dashed line shows the initial thermal pressure.

Fig. 7.  (a) Time averages of thermal pressure, magnetic pressure, and dynamic pressures for each Lagrangian fluid element as a function of the time averaged position of each element. The thick dashed line shows the thermal pressure, the thick solid line shows the magnetic pressure of the $`y`$-component of the magnetic field, the dash-dotted line shows the dynamic pressure of the $`y`$-component of velocity, and the dotted line shows the dynamic pressure of the $`z`$-component of velocity. The thin solid line shows the magnetic pressure of background magnetic field, and the thin dashed line shows the initial thermal pressure. (b) The dynamic pressure for the $`y`$-component of velocity (dash-dotted line) and the magnetic pressure of the $`y`$-component of the magnetic field (solid line) as a function of density. The dotted line shows the scaling expected from the WKB theory.

In contrast to the WKB theory, the dynamic pressure of the $`z`$-component of the velocity is nonzero in our results, although it is smaller than the dynamic pressure of the $`y`$-component of the velocity. According to the small-amplitude theory, the dynamic pressure of the $`y`$-component of the velocity and the magnetic pressure obey equipartition, i.e.,

$$\frac{1}{2}\rho (v_y^{})^2=\frac{B_y^2}{8\pi }.$$
(49)
Figure 7a shows that such an equipartition is almost satisfied between $`z=1H_0`$ and $`z=5H_0`$, although there are deviations near the midplane. This deviation comes from the effect of a nonzero driving force at $`z=0`$ and the symmetrical boundary condition $`B_y(z=0)=0`$. For $`z>5H_0`$, a region which contains a small fraction of the total cloud mass (see Figure 4), the two energies have distinct spatial profiles. Near the interface between the cold and hot material (the cloud boundary), the magnetic pressure decreases more rapidly than expected from the WKB theory. We believe this is due primarily to a standing wave that is set up in the outer cloud, with the cloud boundary acting as a node for $`B_y`$ and an antinode for $`v_y^{}`$. Although many different wave modes are generated by the turbulence in the cloud, only those which satisfy the boundary condition for a standing wave will interfere constructively upon reflection. A transverse standing wave can be set up even though the boundary itself is moving, since the Alfvén speed is much greater in the outer cloud than the vertical speed of the boundary. The influence of the boundary reaches well inside the cloud, since the effective wavelength of Alfvén waves (scaling as $`\rho ^{1/2}`$ for waves of fixed frequency) is also quite large in the low density region. For example, at $`z=5H_0`$, where $`\rho 0.1\rho _0`$, Alfvén waves of the input frequency $`\nu _0=c_{s0}/H_0`$ have wavelength $`\lambda 4.5H_0`$; at $`z=10H_0`$, it increases to $`\lambda 8H_0`$. An important result is that even though the waves are quasilinear in the outer cloud ($`B_yB_0`$), the transverse velocity amplitude is much larger than expected from equipartition arguments.

In Figure 7b, we show the dynamic pressure of the $`y`$-component of velocity and the magnetic pressure of the $`y`$-component of the magnetic field as a function of density in order to compare with the WKB model predictions. We find that the energy in transverse velocities, $`\frac{1}{2}\rho (v_y^{})^2`$, scales approximately as $`\rho ^{1/2}`$ in the outer cloud, the same as the WKB theory, although there is a noticeable upward turn at $`z9H_0`$ due to the standing wave effect. However, the wave magnetic energy $`B_y^2/(8\pi )`$, which is directly responsible for vertical support, decreases more rapidly than in the WKB model in the outer cloud, scaling approximately as $`\rho `$. We note that strict adherence to the WKB predictions is not expected in our model due to nonlinearity, wave dissipation, and the low frequency $`\nu _0=c_{s0}/H`$ of the input turbulence.

#### 3.1.4 Energies in the Cloud

Figure 8 shows the time evolution of various energies in the cloud. These are calculated as follows: kinetic energy of the $`z`$-component of velocity,

$$E_{kz}(t)=_0^{z_f(t)}\frac{1}{2}\rho v_z^2𝑑z;$$
(50)
kinetic energy of the $`y`$-component of velocity,

$$E_{ky}(t)=_0^{z_f(t)}\frac{1}{2}\rho v_y^2𝑑zE_{km},$$
(51)
where $`E_{km}`$ is the kinetic energy of the mean motion of the cloud, i.e.,

$$E_{km}(t)=\frac{1}{2}v_m^2_0^{z_f(t)}\rho 𝑑z;$$
(52)
magnetic energy of the $`y`$-component of magnetic field,

$$E_m(t)=_0^{z_f(t)}\frac{B_y^2}{8\pi }𝑑z;$$
(53)
sum of the above terms,

$$E_T(t)=E_{kz}(t)+E_{ky}(t)+E_m(t).$$
(54)
In the above equations, the integration was done from 0 to $`z_f(t)`$ because we are interested in energies of the cold material. Strictly speaking, the total energy in each case is twice the value we calculate due to the mirror symmetric boundary condition at $`z=0`$. Equation (51) shows that the mean kinetic energy of the cloud is subtracted from the kinetic energy of $`y`$-component of velocity, so that $`E_{ky}`$ means the kinetic energy of the oscillating velocity.

Fig. 8.  Time evolution of various energies in the cloud. The thick solid line shows $`E_{ky}`$, the dotted line shows $`E_{kz}`$, the thin solid line shows $`E_m`$, and dash-dotted line shows $`E_T`$. The values of each energy are smoothed out over every one cycle of the driving force in order to remove periodic oscillations originating from the driving force.

In Figure 8, the thick solid line shows $`E_{ky}`$, the dotted line shows $`E_{kz}`$, the thin solid line shows $`E_m`$, and dashed-dotted line shows $`E_T`$. The values of each energy are smoothed out over every one cycle of the driving force in order to remove periodic oscillations originating from the driving force. Among the energies, $`E_{ky}`$ and $`E_m`$ are comparable to each other, but $`E_{kz}`$ is significantly smaller than the others. The sum of energies $`E_T`$ is almost constant in a logarithmic scale while the input energy is constant ($`t=10t_040t_0`$), though it has some fluctuations. After the driving force is terminated at $`t=40t_0`$, the energies decrease almost exponentially. The decreasing time is $`8.5t_0`$.

### 3.2. Parameter Dependence on the Strength of the Driving Force

#### 3.2.1 Density, Velocity and Pressure

In this paper, we study the effect of changing the strength of the driving force, by changing $`a_d`$ in equation (31).

Figure 9 shows the time evolution of densities for $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. The $`y`$-velocities at $`z=0`$, $`v_y^{}`$, are also shown in the figure. This figure shows that a stronger driving force causes a larger turbulent velocity, which results in a more dynamic evolution of the molecular cloud, including stronger shock waves, and larger excursions of the cloud boundary. However, after terminating the driving force at $`t=40t_0`$, the shock waves dissipate and the clouds shrink in both cases.

Fig. 9.  Evolution of clouds with $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. (a) Time evolution of densities for $`\stackrel{~}{a}_d=20`$. (b) The oscillating part of $`y`$-velocity at $`z=0`$, $`v_y^{}`$, as a function of time for $`\stackrel{~}{a}_d=20`$. (c) Time evolution of densities for $`\stackrel{~}{a}_d=40`$. (d) The oscillating part of $`y`$-velocity at $`z=0`$, $`v_y^{}`$, as a function of time for $`\stackrel{~}{a}_d=40`$. The density plots at various times are stacked with time increasing upward in uniform increments of $`0.2t_0`$.

Figure 10 shows the time average of the density and velocities for each Lagrangian fluid element for $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. The time average is taken between $`t=10t_0`$ and $`t=40t_0`$. This figure also shows that a stronger driving force causes the cloud to move further outward and a larger velocity dispersion within the cloud.

Fig. 10.  Density and velocity dispersions for the cases $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. (a) Time average of the density for each element as a function of the time averaged position of each element for $`\stackrel{~}{a}_d=20`$. The time average is calculated between $`t=10t_0`$ and $`t=40t_0`$. Each Lagrangian fluid element is equally spaced at time $`t=0`$ with spacing $`\mathrm{\Delta }z=0.1H_0`$ starting at $`z=0.01H_0`$. (b) The time averaged root mean square of various velocities for each Lagrangian fluid element plotted versus the time average of the position of each element for $`\stackrel{~}{a}_d=20`$. Open circles show the time averaged root mean square of $`v_z`$, filled circles show that of $`v_y^{}`$, and the dashed line shows the time average of the sound speed $`c_s`$. The solid line shows the root mean square of the time average of the effective one-dimensional velocity dispersion. (c) The time average of the density for each element as a function of the time averaged position of each element for $`\stackrel{~}{a}_d=40`$. (d) The time averaged root mean square of various velocities for each Lagrangian fluid element plotted versus the time average of the position of each element for $`\stackrel{~}{a}_d=40`$.

Figure 11 shows the time average of pressures for $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. When the driving force is strong, the magnetic pressure and dynamic pressure of the $`z`$-component of velocity become significantly larger than the thermal pressure.

Fig. 11.  Time averages of the thermal pressure, magnetic pressure, and dynamic pressures for each Lagrangian fluid element as a function of the time averaged position of each element. The upper panel is for $`\stackrel{~}{a}_d=20`$ and the lower panel is for $`\stackrel{~}{a}_d=40`$. The thick dashed line shows the thermal pressure, the thick solid line shows the magnetic pressure of the $`y`$-component of the magnetic field, the dash-dotted line shows the dynamic pressure of the $`y`$-component of velocity, and the dotted line shows the dynamic pressure of the $`z`$-component of velocity. In both panels, the thin solid line shows the magnetic pressure of background magnetic field, and the thin dashed line shows the initial thermal pressure.

#### 3.2.2 Energy

Figure 12 shows the time evolution of $`E_T`$ for both $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. These values are also smoothed out over every one cycle of the driving force to remove periodic oscillations. In the case of $`\stackrel{~}{a}_d=20`$, the energy is almost constant in a logarithmic scale while the driving amplitude is constant ($`t=10t_040t_0`$). However, in the case of $`\stackrel{~}{a}_d=40`$, the energy is still gradually increasing until $`t=30t_0`$, but afterwards becomes almost constant until $`t=40t_0`$. After terminating the driving force at $`t=40t_0`$, both energies decrease almost exponentially. The energy decreasing times, $`t_d`$, for each parameter, which are estimated by fitting an exponential function, are listed in Table 1. In this study, we terminate the driving force at $`t=40t_0`$ in every case for simplicity. However, we found that the energy decreasing time can vary somewhat depending on when we terminate the driving force. Accounting for this as well as the fitting error of the exponential function, we conclude that the energy decreasing times have a range of variation about $`\pm 2t_0`$ of the values listed in Table 1.

Fig. 12.  Time evolution of $`E_T`$ for both $`\stackrel{~}{a}_d=20`$ and $`\stackrel{~}{a}_d=40`$. The values of each energy are smoothed out over every one cycle of the driving force in order to remove periodic oscillations originating from the driving force.

#### 3.2.3 Correlations between Velocities and Sizes

Here we investigate the correlation between the velocity dispersion and the height of the cloud. Figure 13a shows the time averaged velocity dispersions $`\sigma ^2_t^{1/2}`$ of different Lagrangian fluid elements for different $`\stackrel{~}{a}_d`$, as a function of $`z_t`$. The open circles correspond to Lagrangian fluid elements whose initial positions are located at $`z=2.51`$, which is close to the edge of the cold cloud. The filled circles corresponds to Lagrangian fluid elements whose initial positions are located at $`z=0.61`$, which is approximately the half-mass position of the cold cloud. Each circle corresponds to a different value of $`\stackrel{~}{a}_d`$. These values are summarized in Table 1. The dotted line shows

$$\sigma ^2_t^{1/2}z_t^{0.5}.$$
(55)
This figure shows that the velocity dispersions have a good correlation with the heights of the molecular clouds.

Fig. 13.  Global properties of an ensemble of clouds with different driving strengths $`\stackrel{~}{a}_d`$. (a) Time averaged velocity dispersions of different Lagrangian fluid elements for different $`\stackrel{~}{a}_d`$, as a function of time averaged positions. The open circles correspond to Lagrangian fluid elements whose initial positions are located at $`z=2.51`$, which is close to the edge of the cold cloud. The filled circles correspond to Lagrangian fluid elements whose initial positions are located at $`z=0.61`$, which is approximately the half-mass position of the cold cloud. The dotted line shows $`\sigma ^2_t^{1/2}z_t^{0.5}`$. Each circle can be associated with a particular model in our study by comparison with the numbers in Table 1. (b) Time averaged velocity dispersions as a function of the mean Alfvén velocity of the cloud. The dotted line shows $`\sigma ^2_t^{1/2}\overline{V}_A`$.

Figure 13b shows the correlation between the velocity dispersion and mean Alfvén velocity of the cloud $`\overline{V}_A`$ at the mean position $`z_t`$ of the Lagrangian elements. The dotted line shows

$$\sigma ^2_t^{1/2}\overline{V}_A.$$
(56)
This figure shows that the velocity dispersions have a good correlation with the mean Alfvén velocity, defined by

$$\overline{V}_A\frac{B_0}{\sqrt{4\pi \overline{\rho }}},$$
(57)
where

$$\overline{\rho }=\frac{\mathrm{\Sigma }}{2z_t}$$
(58)
is the mean density and $`\mathrm{\Sigma }`$ is the column density for each Lagrangian element. We discuss the meaning of these correlations in the next section.

## 4. Discussion

### 4.1. Velocity Dispersion and Equilibrium

In Figure 13a, we found that the velocity dispersion obeys the linewidth-size relation

$$\sigma ^2_t^{1/2}z_t^{0.5}.$$
(59)
This relation is satisfied (at both the full-mass and half-mass Lagrangian positions) for an ensemble of the clouds with different strengths of the driving force. It is consistent with observational results of molecular clouds (Larson 1981; Myers 1983; Solomon et al. 1987). Within any individual cloud, $`\sigma ^2_t^{1/2}`$ also increases toward the cloud boundary.

A similar but not identical relation, relating the velocity dispersion at $`z=0`$ to the size of the cloud, is expected in virtually any one-dimensional model, regardless of the form of spatial variation of the pressure within the cloud. This is because an integral over the vertical force-balance equation reveals that the total pressure at $`z=0`$, $`P_{\mathrm{tot},0}`$ must equal the weight of the accumulated gas above, $`\pi G\mathrm{\Sigma }^2/2`$, assuming that the surface pressure is negligible. If $`P_{\mathrm{tot},0}=\rho _0\sigma _{\mathrm{eff},0}^2`$, where $`\sigma _{\mathrm{eff},0}`$ is the effective velocity dispersion at $`z=0`$, and $`\mathrm{\Sigma }=2\rho _0Z`$, in which $`Z`$ is the typical size scale of the cloud, then one naturally obtains

$$\sigma _{\mathrm{eff},0}Z^{1/2}$$
(60)
for an ensemble of clouds of fixed total column density $`\mathrm{\Sigma }`$ but varying $`\sigma _{\mathrm{eff},0}`$. This relation applies to the Spitzer equilibrium state, in which $`\sigma _{\mathrm{eff},0}=c_s`$ and $`Z=H_0`$. It also applies to the equilibrium state calculated by Martin et al. (1997), in which $`\sigma _{\mathrm{eff},0}`$ is the effective velocity dispersion associated with Alfvén waves at $`z=0`$. While a similar relation will inevitably apply to our nonlinear model as well, we note that the relation (59) is more general, in that it relates the velocity dispersion at the position (not $`z=0`$) of a Lagrangian mass element to the position of that element. This more general correlation is a less predictable property of our time-averaged equilibrium state. It is also closer to what is often measured, since optical depth effects may mean that a measured velocity dispersion samples the largest scales of an observed cloud rather than the center.

Equation (59) is strongly related to another relation

$$\sigma ^2_t^{1/2}\overline{V}_A$$
(61)
which is shown in Figure 13b, once more for both the half-mass and full-mass positions for our ensemble of clouds. This relation is also consistent with observational results of molecular clouds (Crutcher 1999; Basu 2000). For example, Figure 1b of Basu (2000) shows an excellent correlation between the line-of-sight component of the large-scale magnetic field $`B_{\mathrm{los}}`$ and $`\sigma \rho ^{1/2}`$ for observed clouds of widely varying values of $`\rho `$, $`B_{\mathrm{los}}`$, and $`\sigma `$, essentially the same correlation as equation (61) if $`B_{\mathrm{los}}B`$ on average.

Here we would like to point out that the relation (61) does not necessarily imply that turbulent motions are due to Alfvén waves, although that is primarily the case in our simulation. The relation is far more general in that it is obtained when the velocity dispersion is caused by any mechanism which results in the clouds being in virial balance between gravitational and turbulent energies, as well as having a mass-to-flux ratio that is approximately uniform from cloud to cloud (Mouschovias 1987; Shu et al. 1987). Rather than imagining a time-independent mean Alfvén speed in a cloud (strictly true in an incompressible medium or a periodic box simulation with a fixed mean density), and turbulent motions becoming comparable to that speed, our simulation reveals that it is the mean Alfvén speed of any cloud which readjusts to a new value as $`\overline{\rho }`$ drops to accommodate a particular level of turbulent driving. In this view, the turbulent dispersion $`\sigma `$ is the more fundamental quantity, dependent on the particular source of driving, and $`\overline{V}_A`$ is a quantity which readjusts to become comparable to $`\sigma `$.

Finally, we note that the relation (61) depends on the definition of the mean Alfvén velocity. If we use the local Alfvén velocity of each particle, we could not get a clear correlation like equation (61). The mean Alfvén velocity is often calculated observationally using the mean density, as in equations (57) and (58), because it is difficult to directly measure a local Alfvén velocity. The local Alfvén velocity in a gravitationally stratified cloud often takes on values much different than the overall mean quantity $`\overline{V}_A`$ used in equation (61). For example, the time averaged local Alfvén velocity near the edge of the cloud is about 2.6 times greater than $`\overline{V}_A`$ for the case of $`\stackrel{~}{a}_d=30`$.

### 4.2. Relation to Linear Model and Chandrasekhar-Fermi Formula

In our simulation, the velocity component parallel to the background magnetic field ($`v_z`$) is generated by the nonlinear effect of the waves. However, the time-averaged magnitude of $`v_z`$ is significantly less than the time-averaged magnitude of $`v_y^{}`$ in all our models (see Figures 6, 7, 10, and 11). Although the waves are nonlinear, the coupling between the Alfvén and slow mode MHD waves is not so strong as to destroy an approximate equipartition between $`B_y^2/(8\pi )`$ and $`\frac{1}{2}\rho (v_y^{})^2`$ throughout much of the cloud, which is the expected result for purely transverse linear Alfvén waves. However, this equipartition does break down in the outer part of the cloud, where standing wave motions are established in which $`B_y`$ has a node at the cloud boundary, and $`v_y^{}`$ has an antinode.

The breakdown of the $`v_y^{}`$ v.s. $`B_y`$ relation in outer cloud means that the original Chandrasekhar-Fermi formula breaks down there as well. If we assume that the dispersion of polarization angle of the magnetic field ($`\delta \theta `$) is related to the velocity dispersion as

$$\delta \theta =\frac{|B_y|}{B_0}=\alpha \frac{|v_y^{}|}{V_A},$$
(62)
we can estimate the strength of the magnetic field of the cloud by using the observational values of polarization angle, density, and velocity dispersion. This yields

$$B_0=\alpha (4\pi \rho )^{1/2}|v_y^{}|(\delta \theta )^1,$$
(63)
where $`\alpha `$ is a nondimensional factor and equals 1 for linear Alfvén waves. The relation (63) with $`\alpha =1`$ was proposed by Chandrasekhar & Fermi (1953) in order to estimate the strength of the background magnetic field in the interstellar medium. The generalized form with $`\alpha 1`$ can be fit to our simulation if we define $`\alpha `$ in our simulation as the square root of the ratio of time averaged magnetic pressure to dynamic pressure for each Lagrangian fluid element, i.e.,

$$\alpha =\left(\frac{B_y^2/(8\pi )_t}{\rho (v_y^{})^2/2_t}\right)^{1/2}.$$
(64)
The resulting distribution of $`\alpha `$ as a function of $`z_t`$ for the case of $`\stackrel{~}{a}_d=30`$ is shown in Figure 14. The time averaged density, which was shown in Figure 4, is also shown as a dashed line to clarify the edge of the cloud. It shows that $`\alpha `$ is close to unity inside the cloud in the region where most of the mass is enclosed (see Figure 4), but it decreases near the surface of the cold cloud, where the standing wave effect becomes important; the minimum value is about $`\alpha =0.23`$. Outside of the cloud ($`z_t>12H_0`$ in Figure 14), $`\alpha =1`$ since the waves are linear there due to the low density and high ambient Alfvén speed. Therefore, while our nonlinear model supports the use of the Chandrasekhar-Fermi formula throughout most of a stratified cloud, we caution against its use with $`\alpha =1`$ near the surface of a cloud.

Fig. 14.  Parameter $`\alpha `$ which appears in the Chandrasekhar-Fermi formula (eq. ), as a function of the time averaged position of fluid elements within the cloud for the case of $`\stackrel{~}{a}_d=30`$. The time averaged density, which was shown in Figure 4, is also shown as a dashed line, mainly to clarify the position of the edge of the cloud.

### 4.3. Global Oscillations

Figure 15 shows the time evolution of the position of a Lagrangian fluid element whose initial position is $`z=2.51`$, for the standard model with driving strength $`\stackrel{~}{a}_d=30`$. This corresponds to the motion of the outer edge of the cloud. The motions resemble longitudinal normal mode oscillations, with the excursions of the outer cloud very similar to free-fall. The dotted line in the figure shows the trajectory for free-fall motion of the fluid element for several different time intervals. The peak of the trajectories are chosen to coincide with a peak of the oscillation in the computed model. The paths are parabolic due to the constant gravity acting on a comoving mass shell in the one-dimensional approximation. Similar to the transverse standing wave that is set up in the outer cloud, the longitudinal motions in this region also resemble a standing wave pattern with an antinode of $`v_z`$ at the cloud boundary. The outermost part of the cloud suffers the greatest displacements, as in the case of a pulsating star. Once the internal driving is discontinued ($`t>40t_0`$), the outer surface also moves inward in nearly a free-fall manner.

Fig. 15.  Time evolution of the position of a Lagrangian fluid element whose initial position is located at $`z=2.51`$ for the case of $`\stackrel{~}{a}_d=30`$. The dotted line shows the trajectories of free-fall motion for the element.

Figure 3 and Figure 15 also show that some residual oscillations remain on the largest scale even after the damping of most of the internal turbulent energy. This effect is reminiscent of the recent observation of Lada et al. (2003) that implied a global oscillation of Barnard 68, a cloud with only subsonic internal turbulence. We hypothesize that Barnard 68 may have dissipated most of the internal turbulence left over from its formation, but that only a large-scale oscillation remains, as in the late stage of our simulation after turbulent driving is discontinued. However, we caution that ours is a one-dimensional model, and globally coherent standing-wave motions need to be demonstrated in a multi-dimensional model.

### 4.4. Dissipation of Energy

The dissipation of energy that we input into nonlinear transverse Alfvén modes is caused though either shock waves or grid-scale dissipation, which is similar to the results of Stone et al. (1998). In addition to the dissipation, a part of the energy escapes from the cloud in our simulation. We measured the Poynting flux, which is the dominant energy flux, at the full-mass position and the 20-percent mass position, which is just outside the region we input the driving force. The time integral of the Poynting flux at the full-mass position is about 30 percent of that at the 20-percent mass position.

In our model, the majority of dissipation occurs due to the transfer of energy to small scales via nonlinear steepening and/or a turbulent cascade, followed by damping on the grid scale though numerical resistivity. This is because the kinetic energy of longitudinal motions, which are the cause of shock waves, is smaller than that of transverse motions. The primary dissipation of Alfvénic motions then, is not due to coupling to slow modes, although this is present at a significant level. This is similar to the finding of the one-dimensional numerical model of Gammie & Ostriker (1996), although the same authors claim that in a multi-dimensional simulation, compressive effects become a major source of dissipation (Stone et al. 1998); however, see Cho & Lazarian (2003), who find that the primary dissipation mechanism in a three-dimensional simulation is not the coupling of Alfvén modes to slow and fast modes.

So if the primary source of dissipation in this simulation is grid-scale dissipation after nonlinear steepening of waves (and perhaps some cascade to smaller wavelengths), then is the dissipation rate a physical effect or an artifact of the simulations? To answer this, we have increased the resolution of our numerical simulations until the dissipation rate is nearly independent of grid size for the case of the typical parameter. Hence, the numerical diffusion is no longer resolution-determined at our current resolution of 50 points per initial scale height in the region where most of the mass is located. If we use a smaller grid size, the final scale of the dissipation is smaller, but the rate of dissipation stays nearly the same. In reality, there must be a physical effect such as ion-neutral friction which would damp the waves at small scales.

The dissipation time of our results are a few crossing times of the time-averaged scales of the clouds (see Table 1). They are a bit longer than estimated from periodic box simulations. We think that this comes from the generation of longer-wavelength modes as the waves travel to low density regions near the cloud’s surface. However, it is already known that one-dimensional simulations have lower dissipation rates than two or three-dimensional ones (Ostriker et al. 2001). Therefore, higher dimensional global simulations will give the final answer in the future.

### 4.5. Future Work

This is the first in a series of papers. In the next paper, we will conduct a complete survey of the effect of important parameters such as $`\beta _0`$ and $`\nu _0`$. Additionally, the study of random, rather than sinusoidal, disturbances will be considered. Moreover, we will study the case of circularly polarized Alfvén waves by including an $`x`$-component of motions in the simulation. A circularly polarized wave is possibly less dissipative than a linearly polarized wave, since in the theoretical limit of an infinite wave train, a flux of circularly polarized waves has no associated magnetic pressure gradient and thereby induces no compression of gas. Although conversion of Alfvén modes to slow modes is already not a dominant source of dissipation in our model, it will be useful to see if the dissipation rate is even less than measured here when circularly polarized waves propagate in a finite-sized cloud.

We also expect to include ion-neutral friction in a future simulation. This effect is expected to be important in molecular clouds, especially in this model as a damping mechanism of Alfvén waves (Kulsrud & Pearce 1969; Zweibel & Josafatsson 1983). It would also change the density structure of a cloud lifted up by wave pressure, as demonstrated in the case of linear Alfvén waves by Martin et al. (1997).

Finally, we believe that a multi-dimensional turbulent simulation in a gravitationally stratified medium is very necessary for the future, especially to obtain more definitive dissipation rates in a stratified cloud and to study the feedback to turbulence due to gravitational contraction and differential rotation.

## 5. Summary

We have performed a numerical simulation of nonlinear MHD waves in a stratified molecular cloud. Our main results are as follows:

1. Due to the effective pressure of MHD turbulence, our one-dimensional cloud is lifted upward and establishes a steady-state characterized by oscillations about a time-averaged equilibrium state. The outer, low density parts of the cloud undergo the largest oscillations, which have the character of free-fall motions. After turbulent driving is discontinued, the cloud falls back toward the initial equilibrium, but some large-scale oscillations of the cloud are the longest-lived modes.
2. The nonlinear effect of the wave propagation results in a significant conversion of energy from Alfvén waves to compressive slow (acoustic) modes and consequent shock formation. However, the energy in transverse modes always remains significantly greater than that in longitudinal modes. Further dissipation of transverse modes takes place through the generation of smaller scale structure through nonlinear steepening or a cascade, culminating in resistive dissipation on the grid-scale.
3. Within each cloud, the magnetic pressure associated with wave motion, $`B_y^2/(8\pi )_t`$, is in approximate equipartition with the kinetic energy in transverse motions $`\rho (v_y^{})^2/2_t`$ in the region containing most of the mass. However, in the low density outer regions, the wave amplitudes have the character of standing waves, such that $`B_y`$ has a node and $`v_y^{}`$ has an antinode. All of this means that the Chandrasekhar-Fermi model applied to a stratified cloud will have a multiplicative constant $`\alpha `$ which is $`1`$ with weak spatial dependence through much of the cloud, but which drops to much lower values near the edge of the cloud.
4. After turbulent driving is discontinued, the dissipation time of the cloud turbulence is $`10t_0`$ or several crossing times of the time-averaged equilibrium state during turbulent driving. These times are up to a few times longer than those found in multi-dimensional periodic simulations. This would be due partly to fewer dissipation avenues in a one-dimensional approximation, but also partly to the generation of long-wavelength modes due to the cloud stratification.
5. For an ensemble of clouds with different levels of internal driving, we find the relation $`\sigma Z^{0.5}`$, where $`\sigma `$ is the time averaged one-dimensional velocity dispersion measured at a size scale $`Z`$ of the time averaged equilibrium state of the cloud.
6. For the same ensemble of clouds, the above relation also implies that $`\sigma \overline{V}_A`$, where $`\overline{V}_A`$ is the mean Alfvén velocity determined from the mean density of the cloud. This relation is in agreement with observed properties of magnetized clouds and cloud cores.

TK acknowledges support due to a Fellowship from the Canadian Galactic Plane Survey, which is supported by a Natural Sciences and Engineering Research Council of Canada (NSERC) Collaborative Research Opportunities Grant. TK also benefited from a Fellowship from SHARCNET, a high-performance computing project funded by ORDCF and CFI/OIT. SB was supported by an individual research grant from NSERC. Numerical computations were carried out mainly on the VPP5000 at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan.


format: text

text:
# Excited heavy baryon masses to order Λ_{𝑄𝐶𝐷}/𝑚_𝑄 from QCD sum rules

## I Introduction

In the past decade continuous progress has been made in the investigation of excited heavy baryons. The lowest lying orbitally excited charmed states $`\mathrm{\Lambda }_c(2593)`$ and $`\mathrm{\Lambda }_c(2625)`$ have been observed by several collaborationslc-exp , the excited states of $`\mathrm{\Xi }_c`$ have also been reported recentlyxc-exp , and more data are expected in the near future. From the theoretical prospect of view those data need to be studied comprehensively. Furthermore, with the collection of more experimental data for exited heavy baryon states it is useful to make some theoretical predictions on their spectroscopies.

Heavy baryons containing a single heavy quark can be described exactly by heavy quark effective theory (HQET) HQET ; review ; KKP  in the heavy quark limit. This fact should be contributed to the spin-flavor symmetry of system comprised of infinitely heavy quarks. HQET has been applied successfully to learn about the properties of heavy mesons and baryons, including the spectroscopy and weak decays. The mass formula for a spin symmetry doublet of heavy baryons up to order $`1/m_Q`$ corrections can be written as

$$M=m_Q+\overline{\mathrm{\Lambda }}\frac{1}{2m_Q}(\lambda _1+d_M\lambda _2),$$
(1)
where parameter $`\overline{\mathrm{\Lambda }}`$ is the effective mass of the light degrees of freedom in the $`m_Q\mathrm{}`$ limit, $`\lambda _1`$ and $`\lambda _2`$ are related to the heavy quark kinetic energy and the chromomagnetic energy of HQET, respectively

$`\lambda _1`$ $`=`$ $`B(v)\overline{h}_v(iD^{})^2h_vB(v),`$
$`d_M\lambda _2`$ $`=`$ $`B(v)\overline{h}_v\sigma _{\mu \nu }{\displaystyle \frac{g_s}{2}}G^{\mu \nu }h_vB(v).`$ (2)
The constant $`d_M`$ characterizes the spin-orbit interaction of the heavy quark and the gluon field, it is zero for singlet and $`1`$, $`\frac{1}{2}`$ for spin 1/2, spin 3/2 doublet, respectively. Thus the splitting of the spin $`1/2`$ and $`3/2`$ doublets is

$$M_{B_Q^{}}^2M_{B_Q}^2=\frac{3}{2}\lambda _2,$$
(3)
where $`B_Q`$, $`B_Q^{}`$ denote spin 1/2 and 3/2 doublet, respectively.

The heavy baryon mass parameters $`\lambda _1`$ and $`\lambda _2`$ play a most significant role in our understanding of the spectroscopy FaNe  and inclusive decay rates Bigi . They must be estimated in some nonperturbative approaches due to the asymptotic freedom property of QCD. A viable approach is the QCD sum rules svzsum  formulated in the framework of HQET hqetsum . This method allows us to relate hadronic observable to QCD parameters via the operator product expansion (OPE) of the correlator. In the case of ground state heavy baryon, predictions on the mass spectroscopy have been made to leading and next-to-leading order in $`\alpha _s`$ Grozin ; Groote  and to order $`1/m_Q`$ Dai ; wang  using QCD sum rule method. As to the excited baryon mass spectroscopy, only results to leading order in $`1/m_Q`$ expansion have been obtained from QCD sum rule CSH ; zhu . In wang  we have calculated the baryonic parameters $`\lambda _1`$ and $`\lambda _2`$ for the ground state $`\mathrm{\Lambda }_Q`$ and $`\mathrm{\Sigma }_Q`$ baryons using QCD sum rules in the HQET. Employing the baryonic currents from CSH  we now derive these parameters for excited $`\mathrm{\Lambda }`$\- and $`\mathrm{\Sigma }`$-type baryons following the same procedure.

The remainder of this paper is organized as follows. In Sec. II.1 we introduce the interpolating currents for excited state heavy baryons and briefly present the two-point sum rules. The direct Laplace sum rules analysis for the matrix elements is presented in Sec. II.2. The Sec. III is devoted to numerical results and our conclusions. Some comments are also available in Sec. III.

## II derivation of the sum rules

### II.1 Heavy baryonic currents and two-point sum rules

In this work we adopt those currents constructed from Bethe-Salpeter equation in Ref.CSH  as

$`j_{\mathrm{\Sigma }_{Qk1}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau C\gamma _5D_t^{\mu _1}q_2^b)\mathrm{\Gamma }^{}h_v^c,`$ (4a)
$`j_{\mathrm{\Lambda }_{Qk0}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau C\gamma _t^\mu D_{t\mu }q_2^b)\mathrm{\Gamma }^{}h_v^c,`$ (4b)
$`j_{\mathrm{\Lambda }_{Qk1}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau Cϵ_{\mu _1\nu \sigma \rho }\gamma _t^\nu v^\rho D_t^\sigma q_2^b)\mathrm{\Gamma }^{}h_v^c,`$ (4c)
$`j_{\mathrm{\Lambda }_{QK1}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau C\gamma _5q_2^b)\mathrm{\Gamma }^{}D_t^{\mu _1}h_v^c,`$ (4d)
$`j_{\mathrm{\Sigma }_{QK0}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau C\gamma _t^\mu q_2^b)\mathrm{\Gamma }^{}D_{t\mu }h_v^c,`$ (4e)
$`j_{\mathrm{\Sigma }_{QK1}}`$ $`=`$ $`ϵ_{abc}(q_1^{Ta}\tau Cϵ_{\mu _1\nu \sigma \rho }\gamma _t^\nu v^\rho q_2^b)\mathrm{\Gamma }^{}D_t^\sigma h_v^c,`$ (4f)
in which $`C`$ is the charge conjugation matrix, $`\tau `$ is the flavor matrix which is antisymmetric for $`\mathrm{\Lambda }`$-type baryon and symmetric for $`\mathrm{\Sigma }`$-type baryon, $`\mathrm{\Gamma }^{}`$s are some gamma matrices, and a, b, c denote the color indices. $`\mathrm{\Gamma }^{}`$ can be chosen co-variantly as

$$\mathrm{\Gamma }^{}=\gamma _{t\mu _1}\gamma _5$$
(5)
for $`\mathrm{\Sigma }_{Qk1},\mathrm{\Lambda }_{Qk1},\mathrm{\Lambda }_{QK1},\mathrm{\Sigma }_{QK1}`$ doublets’ spin 1/2 baryon, and

$$\mathrm{\Gamma }^{}=\mathrm{\Gamma }_{\mu _1\rho _1}=\frac{1}{3}(g_{t\mu _1\rho _1}+\gamma _{t\mu _1}\gamma _{t\rho _1}),$$
(6)
for $`\mathrm{\Sigma }_{Qk1},\mathrm{\Lambda }_{Qk1},\mathrm{\Lambda }_{QK1},\mathrm{\Sigma }_{QK1}`$ doublets’ spin 3/2 baryon, in which $`g_{t\mu _1\rho _1}`$ and $`\gamma _{t\mu _1}`$ are perpendicular to heavy quark velocity $`v`$, defined as $`g_{t\mu \nu }=g_{\mu \nu }v_\mu v_\nu ,\gamma _{t\mu }=\gamma _\mu \text{/}vv_\mu `$. For the singlets $`\mathrm{\Lambda }_{Qk0}`$ and $`\mathrm{\Sigma }_{QK0}`$ the $`\mathrm{\Gamma }^{}`$ is simply unit matrix $`I`$. Notations used here to describe excited state heavy baryons are the same as those used in Ref. KKP ; CSH : k denotes $`l_k=1`$ and $`l_K=0`$ whereas K denotes $`l_k=0`$ and $`l_K=1`$, in which $`l_k`$ is the orbital angular momentum describes the relative motion of the two light quarks and the orbital angular momentum $`l_K`$ describes orbital motion of the center of mass of the two light quarks relative to the heavy quark. Q denotes heavy quark and the number in subscript is the total angular momentum of the light diquark system. As that for ground state baryons, the flavor configuration of $`\mathrm{\Lambda }`$ type baryon is antisymmetric and $`\mathrm{\Sigma }`$ type is symmetric. Depending on the number of derivatives or the form of $`\mathrm{\Gamma }^{}`$, interpolating currents can have forms different from those listed above and may also be used in applicationsCSH ; zhu .

In the following analysis we would use those currents to interpolate excited heavy baryon states. At the leading order of the $`1/m_Q`$ expansion they do not mix with each other even with the same quantum number, but to the next-to-leading order the mixing of interpolating currents will appear. In our subsequent calculations we would only use those currents and did not consider the effect resulting from the mixing of interpolating currents for references Grozin ; mixing1 ; mixing2  have shown that the stability criterion for QCD sum rule applications excludes the existence of interpolating currents mixing or though there does exist the mixing the numerical result will not change drastically compared with the case without mixing.

The baryonic coupling constant in HQET are also needed in our calculation, they are defined in form as follows

$$0jB(v)=Fu,$$
(7)
where $`B(v)`$ denotes excited baryon state and $`u`$ can be the ordinary spinor u or the Rarita-Schwinger spinor $`u_\alpha `$ in the HQET corresponding to spin 1/2 or spin 3/2 doublet, respectively. Irrespective of an irrelevant constant factor in the leading order, the coupling constants for spin 1/2 and spin 3/2 doublet are equivalent for their identical spin-parity of the light degrees of freedom.

In order to determine the effective mass of the excited baryons, we analyze the two-point correlator defined as

$$id^4xe^{ikx}0T\{j(x)\overline{j}(0)\}0=\frac{1+\text{/}v}{2}Tr[\tau \tau ^+]\mathrm{\Pi }(\omega ),$$
(8)
where $`k`$ is the residual momentum and $`\omega =vk`$. For large negative value of $`\omega `$ $`\mathrm{\Pi }(\omega )`$ can be expressed in terms of perturbative and nonperturbative contributions. The nonperturbative effects can be accounted for by including quark and gluon condensates ordered by increasing dimension, which are the series of power corrections in the ”small” $`1/\omega `$ variable. The Borel transformation in the variable $`\omega `$ can help to improve the convergence of these nonperturbative series.

With those interpolating currents listed in Eq. (4) it is straightforward to obtain the two-point sum rules:

$`F_{\mathrm{\Lambda }_{Qk0}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Lambda }_{Qk0}}/T}`$ $`=`$ $`{\displaystyle \frac{18T^8}{\pi ^4}}\delta _7(\omega _c/T)+{\displaystyle \frac{T^4}{4\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{16}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (9a)
$`F_{\mathrm{\Sigma }_{Qk1}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Sigma }_{Qk1}}/T}`$ $`=`$ $`{\displaystyle \frac{90T^8}{\pi ^4}}\delta _7(\omega _c/T){\displaystyle \frac{3T^4}{2^5\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{16}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (9b)
$`F_{\mathrm{\Lambda }_{Qk1}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Lambda }_{Qk1}}/T}`$ $`=`$ $`{\displaystyle \frac{216T^8}{\pi ^4}}\delta _7(\omega _c/T){\displaystyle \frac{T^4}{4\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{8}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (9c)
$`F_{\mathrm{\Lambda }_{QK1}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Lambda }_{QK1}}/T}`$ $`=`$ $`{\displaystyle \frac{216T^8}{\pi ^4}}\delta _7(\omega _c/T){\displaystyle \frac{T^4}{8\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{4}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (9d)
$`F_{\mathrm{\Sigma }_{QK1}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Sigma }_{QK1}}/T}`$ $`=`$ $`{\displaystyle \frac{288T^8}{\pi ^4}}\delta _7(\omega _c/T){\displaystyle \frac{3T^4}{2\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{2}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (9e)
$`F_{\mathrm{\Sigma }_{QK0}}^2e^{\overline{\mathrm{\Lambda }}_{\mathrm{\Sigma }_{QK0}}/T}`$ $`=`$ $`{\displaystyle \frac{504T^8}{\pi ^4}}\delta _7(\omega _c/T){\displaystyle \frac{T^4}{8\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2}{4}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}}.`$ (9f)
In calculations we adopted the gaussian ansatz for the nonlocal quark condensate to get the dimension 6 condensate contribution. Dimension $`D>6`$ condensates are not included. The functions $`\delta _n(\omega _c/T)`$ arise from the continuum subtraction and are defined in wang .

### II.2 Sum rules for $`\lambda _1`$ and $`\lambda _2`$

For the evaluation of the matrix elements $`\lambda _1`$ and $`\lambda _2`$ we consider the three-point correlators as follows

$`i^2{\displaystyle d^4xd^4ye^{ikxik^{}y}0T\{j(x)\overline{h}_v(iD^{})^2h_v(0)\overline{j}(y)\}0}`$ $`=`$ $`{\displaystyle \frac{1+\text{/}v}{2}}Tr[\tau \tau ^+]T_K(\omega ,\omega ^{}),`$
$`i^2{\displaystyle d^4xd^4ye^{ikxik^{}y}0T\{j(x)\overline{h}_v\sigma _{\mu \nu }\frac{g_s}{2}G^{\mu \nu }h_v(0)\overline{j}(y)\}0}`$ $`=`$ $`d_M{\displaystyle \frac{1+\text{/}v}{2}}Tr[\tau \tau ^+]T_S(\omega ,\omega ^{}),`$ (10)
where the coefficients $`T_K(\omega ,\omega ^{})`$ and $`T_S(\omega ,\omega ^{})`$ are analytic functions in the “off-shell energies” $`\omega =vk`$ and $`\omega ^{}=vk^{}`$ with discontinuities for positive values of these variables. Saturating the three-point functions with complete set of baryon states, one can isolate the part of interest, the contribution of the lowest-lying baryon states associated with the heavy-light currents, as one having poles in both the variables $`\omega `$ and $`\omega ^{}`$ at the value $`\omega =\omega ^{}=\overline{\mathrm{\Lambda }}`$.

Confining us to take into account these leading contributions of perturbation and the operators with dimension $`D6`$ in OPE, the relevant diagrams in our theoretical calculations for the kinetic energy are shown in Fig. 1. The relevant diagrams for the chromo-magnetic interaction do not differ from those for the ground state baryons wang , so we do not show them here. Using dispersion relations $`T_K(\omega ,\omega ^{})`$ and $`T_S(\omega ,\omega ^{})`$ can be casted into the form of integrals of the double spectral densities. Following Refs. IWfun ; DHHL ; Neubert92D , introduce new variables $`\omega _+=\frac{1}{2}(\omega +\omega ^{})`$ and $`\omega _{}=\omega \omega ^{}`$, perform the integral over $`\omega _{}`$, assume quark–hadron duality in $`\omega _+`$, and employ Borel transformation $`B_\tau ^\omega `$, $`B_\tau ^{}^\omega ^{}`$ to suppress the continuum contributions and subtractions, we then obtained the desired sum rules. Considered the symmetry of the correlator it is natural to set the parameters $`\tau `$, $`\tau ^{}`$ to be the same and equal to $`2T`$, where $`T`$ is the Borel parameter of the two-point functions. We ended up with the set of sum rules

$`\lambda _2^{\mathrm{\Sigma }_{QK1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^{11}\mathrm{\hspace{0.33em}3}}{\pi ^4}}{\displaystyle \frac{\alpha _s}{\pi }}T^{10}\delta _9(\omega _c/T){\displaystyle \frac{2T^6}{\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{16m_0^2T^2\alpha _s}{3\pi }}\overline{q}q^2e^{\frac{m_0^2}{16T^2}},`$ (11a)
$`\lambda _2^{\mathrm{\Lambda }_{QK1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^{10}\mathrm{\hspace{0.33em}3}}{\pi ^4}}{\displaystyle \frac{\alpha _s}{\pi }}T^{10}\delta _9(\omega _c/T){\displaystyle \frac{3T^6}{\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{4m_0^2T^2\alpha _s}{3\pi }}\overline{q}q^2e^{\frac{m_0^2}{16T^2}},`$ (11b)
$`\lambda _2^{\mathrm{\Sigma }_{Qk1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^8\mathrm{\hspace{0.33em}3\hspace{0.33em}5}}{\pi ^4}}{\displaystyle \frac{\alpha _s}{\pi }}T^{10}\delta _9(\omega _c/T){\displaystyle \frac{T^6}{2\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{m_0^2T^2\alpha _s}{3\pi }}\overline{q}q^2e^{\frac{m_0^2}{16T^2}},`$ (11c)
$`\lambda _2^{\mathrm{\Lambda }_{Qk1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^9\mathrm{\hspace{0.33em}3\hspace{0.33em}5}}{\pi ^4}}{\displaystyle \frac{\alpha _s}{\pi }}T^{10}\delta _9(\omega _c/T){\displaystyle \frac{3T^6}{\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{4m_0^2T^2\alpha _s}{3\pi }}\overline{q}q^2e^{\frac{m_0^2}{16T^2}},`$ (11d)
$`\lambda _1^{\mathrm{\Lambda }_{QK1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^6\mathrm{\hspace{0.33em}3}^3\mathrm{\hspace{0.33em}5}T^{10}}{\pi ^4}}\delta _9(\omega _c/T)+{\displaystyle \frac{19T^6}{2\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{5m_0^4}{16}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (11e)
$`\lambda _1^{\mathrm{\Sigma }_{QK0}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^6\mathrm{\hspace{0.33em}3}^4\mathrm{\hspace{0.33em}5}T^{10}}{\pi ^4}}\delta _9(\omega _c/T)+{\displaystyle \frac{113T^6}{4\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{5m_0^4}{16}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (11f)
$`\lambda _1^{\mathrm{\Sigma }_{QK1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^8\mathrm{\hspace{0.33em}3}^2\mathrm{\hspace{0.33em}5}T^{10}}{\pi ^4}}\delta _9(\omega _c/T){\displaystyle \frac{3T^6}{\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{5m_0^4}{8}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (11g)
$`\lambda _1^{\mathrm{\Sigma }_{Qk1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^4\mathrm{\hspace{0.33em}3}^2\mathrm{\hspace{0.33em}17}T^{10}}{\pi ^4}}\delta _9(\omega _c/T)+{\displaystyle \frac{21T^6}{8\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{3m_0^4}{64}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (11h)
$`\lambda _1^{\mathrm{\Lambda }_{Qk0}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^4\mathrm{\hspace{0.33em}3}^2\mathrm{\hspace{0.33em}7}T^{10}}{\pi ^4}}\delta _9(\omega _c/T){\displaystyle \frac{17T^6}{8\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{3m_0^4}{64}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}},`$ (11i)
$`\lambda _1^{\mathrm{\Lambda }_{Qk1}}F^2e^{\overline{\mathrm{\Lambda }}/T}`$ $`=`$ $`{\displaystyle \frac{2^6\mathrm{\hspace{0.33em}3}^2\mathrm{\hspace{0.33em}11}T^{10}}{\pi ^4}}\delta _9(\omega _c/T){\displaystyle \frac{5T^6}{2\pi ^2}}{\displaystyle \frac{\alpha _s}{\pi }}G^2+{\displaystyle \frac{3m_0^4}{32}}\overline{q}q^2e^{\frac{m_0^2}{8T^2}}.`$ (11j)
The unitary normalization of flavor matrix $`Tr[\tau \tau ^+]=1`$ has been applied to get those sum rules, as what has been done for the two-point sum rules.

## III numerical results and conclusions

In the following analysis, the standard value for the condensates are adopted svzsum . From two-point sum rules the effective mass can be obtained via a derivative of the Borel parameter as $`\overline{\mathrm{\Lambda }}=T^2d\mathrm{ln}E/dT`$, where $`E`$ denotes the right hand side of the obtained two-point sum rule. Then comply to the standard procedure of sum rule analysis, we changed the continuum threshold $`\omega _c`$ and Borel parameter $`T`$ to find the optimal stability window and the numerical value of the effective mass will be determined within this window. For those sum rules obtained above, we found the typical value for the continuum threshold is $`\omega _c1.6\text{GeV}`$ and typical interval for the Borel parameter is $`\mathrm{\Delta }T0.4\text{GeV}`$, which is a narrower one than the window for the ground state baryon. The only exception of this assertion is $`\mathrm{\Lambda }_{Qk0}`$, for which we found a much lower continuum threshold $`\omega _c1.2\text{GeV}`$ and a narrower interval $`\mathrm{\Delta }T0.3\text{GeV}`$. Also it is worth noting that for $`\mathrm{\Lambda }_{Qk0}`$, the main contribution does not come from the perturbative part. The condensate ones play an important role in the determination of stability window. Indeed, if one keeps only gluon condensate contribution, there will be no stability window at all, just like the case in Ref.CSH . But if assuming the perturbative dominance and omitting condensate contributions, we only ended up with a better stability and the numerical value is almost exactly the same. For the other sum rules for the effective mass the case is different. The dominant contribution to those sum rules comes from the perturbative part, and the dimension 4 and dimension 6 operators in the OPE only amount to 20% and 10% within the stability windows, respectively. The numerical results are presented in Fig. 3. For the aim of clarity we give our numerical average for the effective masses in Table 1.

In order to get the numerical results for two $`1/m_Q`$ order parameters, we divide our three-point sum rules by two-point functions to obtain $`\lambda _1`$ and $`\lambda _2`$ as functions of the continuum threshold $`\omega _c`$ and the Borel parameter $`T`$. This procedure can eliminate the systematic uncertainties and cancel the parameter $`\overline{\mathrm{\Lambda }}`$. From the experiences of QCD sum rule applications in the field of heavy quark physics it is well known that three-point sum rule receives heavier contamination from the continuum modes than two-point one, and the stability is not as good as that for two-point sum ruleBaBr ; Grozin ; wang ; mixing2 . Our results for $`\lambda _1`$ sum rules of those excited heavy baryons are shown in Fig. 3 and in Fig. 4 results for $`\lambda _2`$ sum rules are presented. The typical value of continuum threshold is $`\omega _c1.8\text{GeV}`$, $`\omega _c2.4\text{GeV}`$ for $`\lambda _1`$ and $`\lambda _2`$ sum rules, respectively. And the typical interval of the stability window is also $`\mathrm{\Delta }T0.3\text{GeV}`$ for both. The $`\lambda _1`$ sum rule for the $`\mathrm{\Sigma }_{QK0}`$ baryon is an exception of above statement, for which the typical value of continuum threshold is $`\omega _c1.4\text{GeV}`$ and the interval is $`\mathrm{\Delta }T0.2\text{GeV}`$. As for the convergence of the OPE within the stability windows for those sum rules we would like to state some facts. For the $`\lambda _1`$ sum rules, except the case of $`\mathrm{\Lambda }_{Qk1}`$, in which the dimension 4 and 6 operators give rise to contributions almost equal to that of the perturbative part, all other sum rules behave well, in which the contributions of the dimension 4 and 6 operators amount to 50% and 20% of that of the perturbative part. On the other hand, for the $`\lambda _2`$ sum rules, the dimension 4 operator still contributes to almost 50% compared to that of the perturbative part, but the dimension 6 operator gives only a negligible contribution, typically less than 10%. The numerical results are listed in Table 1, too.

With those values and heavy quark masses given in wang , $`m_c=1.41\text{GeV}`$ and $`m_b=4.77\text{GeV}`$, masses of excited heavy baryons to order $`1/m_Q`$ can be obtained immediately. We give those masses in Table 2. Our results are comparable to the prediction on the excited heavy baryon masses obtained by using quark potential model mass .

The splitting between spin 1/2 and spin 3/2 doublet can be obtained by multiplying $`\lambda _2`$ by a factor $`3/2`$, cf. Eq. (3). It is $`0.13\pm 0.02\text{GeV}^2`$, $`0.16\pm 0.05\text{GeV}^2`$, $`0.20\pm 0.02\text{GeV}^2`$ and $`0.19\pm 0.04\text{GeV}^2`$ for $`\mathrm{\Lambda }_{QK1}`$, $`\mathrm{\Lambda }_{Qk1}`$, $`\mathrm{\Sigma }_{QK1}`$ and $`\mathrm{\Sigma }_{Qk1}`$ doublets, respectively. The approximately equal value for $`\mathrm{\Lambda }_{QK1}`$ and $`\mathrm{\Lambda }_{Qk1}`$, $`\mathrm{\Sigma }_{QK1}`$ and $`\mathrm{\Sigma }_{Qk1}`$ may be interpreted as the signal which implies that current mixing effect cannot be large. If taking the middle value as theoretical predictions for the physical state, then the splitting for excited baryon states are

$`\mathrm{\Lambda }_{Q1}^{\mathrm{\hspace{0.17em}2}}\mathrm{\Lambda }_{Q1}^2`$ $`=`$ $`0.15\pm 0.03\text{GeV}^2,`$
$`\mathrm{\Sigma }_{Q1}^{\mathrm{\hspace{0.17em}2}}\mathrm{\Sigma }_{Q1}^2`$ $`=`$ $`0.20\pm 0.03\text{GeV}^2.`$ (12)
Based on the current experimental datadata , the splitting of excited $`\mathrm{\Lambda }_c`$ doublet is $`0.17\text{GeV}^2`$, which is in agreement with our theoretical result. When it is scaled up to the bottom quark mass scale there will be a factor $`0.8`$ approximately due to the renormalization group improvement.

To conclude, we have calculated the $`1/m_Q`$ order corrections to the excited heavy baryon masses from QCD sum rules within the framework of the HQET. From thus obtained spectrum for the c quark case, we found that $`\mathrm{\Lambda }_{Qk0}`$ and $`\mathrm{\Sigma }_{Qk1}`$ baryons lie $`600\text{MeV}`$ above the ground state baryon $`\mathrm{\Lambda }_c`$, while $`\mathrm{\Lambda }_{Qk1}`$, $`\mathrm{\Lambda }_{QK1}`$, $`\mathrm{\Sigma }_{QK1}`$ lie $`800\text{MeV}`$ above $`\mathrm{\Lambda }_c`$ and for $`\mathrm{\Sigma }_{QK1}`$ baryon this value is $`900\text{MeV}`$, typically with an error $`300\text{MeV}`$. When it comes to the b quark case, the result is that $`\mathrm{\Lambda }_{Qk0}`$ and $`\mathrm{\Sigma }_{Qk1}`$ baryons lie $`300\text{MeV}`$ above the ground state baryon $`\mathrm{\Lambda }_b`$, while $`\mathrm{\Lambda }_{Qk1}`$, $`\mathrm{\Lambda }_{QK1}`$, $`\mathrm{\Sigma }_{QK1}`$ and $`\mathrm{\Sigma }_{QK1}`$ lie $`550\text{MeV}`$ above $`\mathrm{\Lambda }_b`$, for which the typical error is $`200\text{MeV}`$. For the c quark case, $`1/m_Q`$ order corrections are $`400\text{MeV}`$, which is not a small one due to the large value of the kinetic energy. When it comes to the b quark case those corrections will be suppressed by the still larger b quark mass. Our theoretical predictions for the doublet splitting are in agreement with the current experimental data.

###### Acknowledgements.

We would like to thank Y.B. Dai for helpful discussions. This work was supported in part by the National Natural Science Foundation of China under Contract No. 10275091.

Figure Captions

Fig. 1. Non-vanishing diagrams for the kinetic energy. The kinetic energy operator is denoted by a white square, the interpolating baryon currents by black circles. Heavy-quark propagators are drawn as double lines.

Fig. 2. Sum rules of the effective mass $`\overline{\mathrm{\Lambda }}`$ for: (a) $`\mathrm{\Lambda }_{Qk0}`$, (b) $`\mathrm{\Lambda }_{Qk1}`$, (c) $`\mathrm{\Sigma }_{Qk1}`$, (d) $`\mathrm{\Lambda }_{QK1}`$, (e) $`\mathrm{\Sigma }_{QK1}`$ and (f) $`\mathrm{\Sigma }_{QK0}`$ baryons. The different choices of the continuum threshold $`\omega _c`$ corresponding to different curves are designated in individual figures respectively. Curves are plotted against the Borel parameter $`T`$.

Fig. 3. Sum rules of the kinetic energy for: (a) $`\mathrm{\Lambda }_{Qk0}`$, (b) $`\mathrm{\Lambda }_{Qk1}`$, (c) $`\mathrm{\Sigma }_{Qk1}`$, (d) $`\mathrm{\Lambda }_{QK1}`$, (e) $`\mathrm{\Sigma }_{QK1}`$ and (f) $`\mathrm{\Sigma }_{QK0}`$ baryons. Others are the same as those in Fig. 2.

Fig. 4. Sum rules of the chromo-magnetic interaction for: (a) $`\mathrm{\Lambda }_{Qk1}`$, (b) $`\mathrm{\Sigma }_{Qk1}`$, (c) $`\mathrm{\Lambda }_{QK1}`$ and (d) $`\mathrm{\Sigma }_{QK1}`$ baryons. Others are the same as those in Fig. 2.