text:
# Time-Dependent, Multifluid, Magnetohydrodynamic Shock Waves with Grain Dynamics I. Formulation and Numerical Tests

## 1 Introduction

Shock waves in weakly-ionized plasmas usually exhibit a multifluid structure, in which the charged and neutral particles behave as separate, interpenetrating fluids. If the shock speed is less than the speed at which the charged fluid communicates compressive disturbances, the charged particles will be accelerated ahead of the neutrals in a “magnetic precursor,” and the neutrals accelerated, compressed and heated downstream by collisions with streaming charged particles (Mullan 1971; Draine 1980, henceforth D80). The flow is continuous (a C- or C shock) if the length scale for acceleration of the neutrals is large compared to the scale for cooling; otherwise the flow undergoes a discontinuous transition (a J shock) at a viscous subshock (D80; Chernoff 1987; Roberge & Draine 1990). The dynamics, chemistry, and emission from multifluid shock waves can differ profoundly from those of their single-fluid counterparts. Consequently, they have been the subject of numerous investigations (see Draine & McKee 1993 for the most recent comprehensive review).

The energy and momentum of the “charged fluid” are usually dominated by the magnetic field, the charged particles acting mainly to provide energy and momentum transfer between the field and the neutral fluid. At low densities the coupling is dominated by ion-neutral scattering and shock waves can be modeled as a three-fluid system: the neutral gas, the ion/magnetic field fluid, and the fluid of electrons (which move with the ions but need not have the same temperature; see D80). However when the density of the preshock medium exceeds $`n10^5`$cm-3,111We use $`n`$ to denote the density of H nuclei in all forms. the energy and momentum transfer is dominated by collisions between neutral particles and charged dust grains (Elmegreen 1979; D80; Nakano & Umebayashi 1980; Ciolek & Mouschovias 1993). If the grains are well coupled to the magnetic field, they can be incorporated into the 3-fluid picture as an additional species of ion. This is a reasonable approximation for PAHs and very small grains with radii $`a_{\mathrm{sg}}10^7`$ cm (henceforth referred to collectively as “small grains”) at densities up to $`10^9`$ (Neufeld & Hollenbach 1994). If the small grains are present in dense gas with a fractional abundance not much less than $`10^7`$ (the abundance inferred for the diffuse ISM, see Li & Draine 2001 and references therein), they will dominate the field-neutral coupling and a 3-fluid model of the dynamics is appropriate (e.g., Kaufman & Neufeld 1996a,b).

However if small grains are depleted in dense gas, as seems likely, then scattering by large grains with $`a_\mathrm{g}0.1`$ µm will dominate the field-neutral coupling. Because the large grains are only partially tied to the magnetic field for $`n10^5`$cm-3, they must be treated as an additional fluid (Chernoff 1985; Havnes, Hartquist & Pilipp 1987; Ciolek & Mouschovias 1993). The problem is further complicated by the fact that large grains have different charge states, a continuous spectrum of sizes, and masses many orders of magnitude larger than those of the ions. Draine (D80) developed an approximate description of the grain dynamics by treating the grains as test particles and neglecting their inertia; in this approximation each grain drifts through the neutrals with a speed that can be determined by balancing the gas drag and electromagnetic forces. Chernoff (1985) treated the grains as a separate fluid, included their inertia, and demonstrated that the grains can undergo a J shock. Pilipp, Hartquist & Havnes (1990) and Pilipp & Hartquist (1994) carried out numerical simulations of steady 4-fluid shocks, emphasizing the importance of grain charging processes and the dynamical effects of the grain charge and current. Wardle (1998, see also Wardle & Ng 1999) has formulated equations of motion valid for steady multifluid shocks with an arbitrary number of grain fluids, in a form that is valid for very large densities, and has discussed optimal procedures for representing the spectrum of large grain sizes by a single, effective size.

This is the first in a series of papers on the effects of dust on the formation, propagation, and structure of nonlinear MHD waves and shocks in dense molecular clouds, bipolar outflows, and protoplanetary disks. Here we describe a dynamical model that incorporates and extends previous work in several ways: (i) We allow for time-dependent flow. Prior work on time-dependent, multifluid shocks (Tóth 1994; Smith & Mac Low 1997; Mac Low & Smith 1997; Stone 1997; Chieze, Pineau des Forets & Flower 1998) has been based on the 3-fluid model. We are interested in time-dependent solutions because they permit us to study the formation of magnetic precursors, because of their possible relevance to bipolar outflows (Flower & Pineau des Forets 1999), and because time-dependent simulations are a robust way to find steady solutions (e.g., steady J shocks) that are otherwise hard to compute. (ii) We include the inertia of large grains, which has been neglected in all prior studies except Chernoff’s (1985). Neglecting the inertia is a good approximation when the gas drag time for large grains is small compared to all other dynamical time scales. While this is usually appropriate for the dynamical time of the neutrals, it is often a poor approximation for the charged fluid (see §2.3.3). A realistic treatment of grain dynamics is necessary for studies on the formation of magnetic precursors: Draine & McKee (1993) pointed out that magnetic precursors could not form (and C shocks therefore could not exist) for shocks faster than $`20\text{km s}\text{-1}`$ if large grains are able to “load” the magnetic field lines efficiently. The code we describe was designed in part to address this issue. (iii) We include small grains as a separate fluid. This permits us to model grain charging realistically and to examine various scenarios for small-grain depletion in dense gas. (iv) We include charge fluctuations by treating the grains in different charge states (3 each for small- and large grains) as separate fluids. Charge fluctuations exchange mass, momentum, and energy between the charged and neutral grain fluids, with implications for the loading of magnetic field lines. Nevertheless, the scope of our calculations is limited in some important respects: We do not attempt to describe the continuous size distribution of the large grains, which are assumed to have uniform radii. We have also restricted our initial investigations to one-dimensional geometry so that, for example, time-dependent simulations of Wardle instabilities with dust are beyond our present capabilities.

The plan of this paper is as follows. In § 2 we give the equations of motion and discuss the conditions where various approximations are likely to be valid. In § 3 we describe our numerical methods and summarize a suite of benchmark calculations that establish the accuracy of our code. Our results are summarized in § 4.

## 2 Formulation

### 2.1 Modeling Assumptions

We consider plane-parallel disturbances propagating along the $`+z`$ direction of a Cartesian frame, $`(x,y,z)`$, defined so that the undisturbed plasma at $`z=+\mathrm{}`$ is at rest. Fluid variables are assumed to depend only on $`z`$ and $`t`$. The plasma is threaded by a magnetic field of the form $`𝑩=B_x(z,t)\widehat{𝒙}`$. That is we restrict the discussion to perpendicular shocks; oblique shocks will be described elsewhere. Fluid motions in the $`x`$-direction (i.e., along the field) are uncoupled from those in the $`y`$ and $`z`$ directions and are ignored.

We model the plasma as a system of 9 interacting fluids. Fluid variables are labeled by a subscript f, where $`\mathrm{f}=\mathrm{n}`$, i, e, sg-, sg0, sg+, g-, g0, and g+ denote the neutral fluid, ions, electrons, small grains with charge $`e`$, $`0`$, or $`+e`$, and large grains with charge $`e`$, $`0`$, or $`+e`$, respectively. The particles of fluid f have mean mass $`m_\mathrm{f}`$, number density $`n_\mathrm{f}`$, mass density $`\rho _\mathrm{f}`$, flow velocity $`𝒗`$f, pressure $`P_\mathrm{f}`$, and internal energy density $`E_\mathrm{f}`$, where $`E_\mathrm{f}`$ includes the contributions of thermal motions, rotational and vibrational excitation, dissociation, etc.

The neutral and ion fluids are composed of $`N_\mathrm{n}`$ and $`N_\mathrm{i}`$ distinct chemical species, respectively, so that

$$n_\mathrm{n}=\underset{\alpha =1}{\overset{N_\mathrm{n}}{}}n_{n\alpha }$$
(1)
and similarly for the ions. Realistic simulations require chemical networks elaborate enough to predict the abundances of important coolants (primarily H2O, CO, and H2 in dense plasmas; see Neufeld & Kaufman 1993; Neufeld, Lepp & Melnick 1995) and to model changes in the fractional ionization (which can profoundly affect the hydrodynamics; see Flower, Pineau des Forêts & Hartquist 1985). In subsequent papers we include chemical networks tailored to the particular applications of interest to us. Here we are interested mainly in establishing the accuracy of our numerical code. The benchmark calculations we present were obtained using a simple analytic expression for the radiative cooling rate (see §3); the chemistry of coolants is therefore not required. To model the ionization we adopt a drastically simplified “chemistry” with just 4 neutral species and 2 ions. The neutral fluid is composed of molecular hydrogen with $`n_{\mathrm{H}_2}=0.5n`$, atomic helium with $`n_{\mathrm{He}}=0.1n`$, plus a generic heavy atom, a, and molecule, m. The ion fluid is composed of the atomic ion a+ and molecular ion m+. Since the dominant ions in the dense astrophysical plasmas relevant to star formation and protostellar disks are HCO+, Na+, and Mg+, we take the ion mass to be $`m_{\mathrm{a}^+}=m_{\mathrm{m}^+}=m_\mathrm{i}=25`$ amu. The chemical abundances are determined by a network of 20 gas-phase and grain surface reactions (Appendix A). The fractional ionization predicted by our simple network is in good agreement with more elaborate models for $`n10^5`$cm-3 (Fig. 1).

The large and small grains are assumed to have uniform radii $`a_\mathrm{g}`$ and $`a_{\mathrm{sg}}`$, where $`a_\mathrm{g}10^5`$ cm and $`a_{\mathrm{sg}}10^7`$ cm are adjustable parameters. Modifying our code to include a range of grain sizes would be straightforward but is beyond our current computational resources. The total abundances of large and small grains, summed over all charge states, are assumed constant.

### 2.2 Governing Equations

#### 2.2.1 Neutral Particles

The mass, momentum, and energy equations for the neutral fluid are

$$\frac{\rho _\mathrm{n}}{t}+\frac{}{z}\left(\rho _\mathrm{n}v_{\mathrm{n}z}\right)=S_\mathrm{n},$$
(2)
$$\frac{}{t}\left(\rho _\mathrm{n}v_{\mathrm{n}z}\right)+\frac{}{z}\left(\rho _\mathrm{n}v_\mathrm{n}^2+P_\mathrm{n}+Q_{\mathrm{visc}}\right)=F_{\mathrm{n}z},$$
(3)
$$\frac{}{t}\left(\rho _\mathrm{n}v_{\mathrm{n}y}\right)+\frac{}{z}\left(\rho _\mathrm{n}v_{\mathrm{n}y}v_{\mathrm{n}z}\right)=F_{\mathrm{n}y},$$
(4)
and

$$\frac{}{t}\left(\frac{1}{2}\rho _\mathrm{n}v_\mathrm{n}^2+E_\mathrm{n}\right)+\frac{}{z}\left[\left(\frac{1}{2}\rho _\mathrm{n}v_\mathrm{n}^2+E_\mathrm{n}+P_\mathrm{n}+Q_{\mathrm{visc}}\right)v_{\mathrm{n}z}\right]=𝑭\text{n}\mathbf{}𝒗\text{n}+G_\mathrm{n}\mathrm{\Lambda }_\mathrm{n}$$
(5)
respectively, where $`v_\mathrm{n}^2=v_{\mathrm{n}z}^2+v_{\mathrm{n}y}^2`$; the “source terms” $`S_\mathrm{n}`$, $`𝑭`$n and $`G_\mathrm{n}`$ are the net rates per volume at which mass, momentum and thermal energy are added to the neutrals by interactions with the other fluids, and $`\mathrm{\Lambda }_\mathrm{n}`$ is the rate of radiative cooling. The benchmark tests presented in §3 were carried out for an idealized fluid with $`\mathrm{\Lambda }_\mathrm{n}`$ prescribed analytically and $`E_\mathrm{n}=P_\mathrm{n}/\left(\gamma 1\right)`$, with constant $`\gamma `$. Future papers will incorporate a realistic calculation of radiative cooling by CO, H2O, and H2 (Neufeld & Kaufman 1993; Neufeld, Lepp & Melnick 1995) and the internal energy of H2 accounting for its time-dependent excitation and H2 dissociation. Our model includes artificial viscosity to smear out jump fronts in the neutral flow. We use the von Neumann-Richtmyer prescription for the pseudo-viscous stress,

$$Q_{\mathrm{visc}}=\{\begin{array}{cc}q\left(\mathrm{\Delta }z\right)^2\rho _\mathrm{n}\left|v_n/z\right|^2& \mathrm{if}v_n/z<0,\hfill \\ & \\ 0& \mathrm{otherwise}\hfill \end{array}$$
(6)
(see Richtmyer & Morton 1967), where $`\mathrm{\Delta }z`$ is the grid spacing and $`q`$ is an adjustable dimensionless parameter. If present, a jump front will be smeared over $`q`$ grid cells.

#### 2.2.2 Electrodynamics

We assume that the magnetic field is frozen into the fluid of ions and electrons. (The validity of this approximation is discussed in §2.3.) In the frozen field approximation, the electrons and ions move with a common velocity,

$$𝒗\text{e}=𝒗\text{i},$$
(7)
the electric field satisfies

$$𝑬\frac{𝒗\text{i}}{c}\mathbf{\times }𝑩,$$
(8)
and the Maxwell equations reduce to the induction equation,

$$\frac{B_x}{t}+\frac{}{z}\left(B_xv_{\mathrm{i}z}\right)=0,$$
(9)
plus the requirement of macroscopic charge neutrality,

$$n_\mathrm{e}=n_\mathrm{i}+n_{\mathrm{g}^+}+n_{\mathrm{sg}^+}n_\mathrm{g}^{}n_{\mathrm{sg}^{}}.$$
(10)
#### 2.2.3 Electrons

We do not solve mass or momentum equations for the electrons because $`n_\mathrm{e}`$ and $`𝒗`$e are functions of the other fluid variables in the frozen-field approximation (cf. eqs.  and ). The electron energy equation is omitted in favor of Chernoff’s (1987) approximation,

$$G_\mathrm{e}=\mathrm{\Lambda }_\mathrm{e},$$
(11)
which says that the electron fluid rapidly reaches thermal balance due to its small heat capacity. The shock structure depends on the electron temperature via the rate coefficients for electron impact excitation of molecules (which contribute to the radiative cooling) and the rate coefficients for recombination and grain charging (which influence the ionization; see Appendix A). The solutions presented in §3 are insensitive to $`T_\mathrm{e}`$ because the cooling is prescribed analytically and the flow time is much less than the timescale for changes in the ionization. We therefore adopt the (unrealistic) approximation $`T_\mathrm{e}=0.15T_\mathrm{i}`$ for expediency; this approximation is based on the results of Draine, Roberge, & Dalgarno (1983, hereafter DRD83). Future papers on realistic simulations will incorporate a detailed description of electron heating and cooling including the atomic and molecular processes described in DRD83.

#### 2.2.4 Ions

We solve the mass equation for the ions,

$$\frac{\rho _\mathrm{i}}{t}+\frac{}{z}\left(\rho _\mathrm{i}v_{\mathrm{i}z}\right)=S_\mathrm{n},$$
(12)
but not the momentum or energy equations. We calculate the momentum of the ion fluid by (i) neglecting the inertia of the ions and electrons and (ii) requiring the total electric current to be consistent with the magnetic field gradient. We assume that the ion and electron velocities are determined solely by the balance between electromagnetic forces and drag due to elastic collisions with the neutrals. Then requirement (i) implies

$$0=en_\mathrm{e}\left(𝑬+\frac{𝒗\text{e}}{c}\mathbf{\times }𝑩\right)+\frac{\rho _\mathrm{e}}{\tau _{\mathrm{e},\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗\text{e}\right)$$
(13)
and

$$0=en_\mathrm{i}\left(𝑬+\frac{𝒗\text{i}}{c}\mathbf{\times }𝑩\right)+\frac{\rho _\mathrm{i}}{\tau _{\mathrm{i},\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗\text{i}\right),$$
(14)
where $`\tau _{\mathrm{i},\mathrm{j}}^\mathrm{d}`$ is the drag time for species i due to elastic scattering with species j (a function of the fluid variables, see DRD83).222Notice that equations (13) and (14) would yield the nonsensical solution $`𝒗\text{e}=𝒗\text{i}=𝒗\text{n}`$ if one naively inserted expression (8) for the electric field. The meaning of this result is that the drag force must vanish in the frozen-field approximation, where gas drag is taken to be negligible compared with electromagnetic forces. Adding expressions (13) and (14) gives

$$e\left(n_\mathrm{i}n_\mathrm{e}\right)𝑬+\frac{e}{c}\left(n_\mathrm{i}𝒗\text{i}n_\mathrm{e}𝒗\text{e}\right)\mathbf{\times }𝑩+\frac{\rho _\mathrm{i}}{\tau _{\mathrm{i},\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗\text{i}\right)0,$$
(15)
where we have neglected the drag force on the electron fluid compared to the drag force on the ions.333The ratio of the drag forces is of order $`(m_\mathrm{e}/m_\mathrm{n})\sigma _{\mathrm{en}}\overline{v}_\mathrm{e}/\sigma v_{\mathrm{in}}10^3`$, where $`\sigma _{\mathrm{en}}`$ is the geometrical cross section, $`\overline{v}_\mathrm{e}`$ is the electron thermal velocity, and $`\sigma v_{\mathrm{in}}`$ is the ion-neutral momentum transfer rate coefficient (DRD83). The quantities $`n_\mathrm{i}n_\mathrm{e}`$ and $`n_\mathrm{i}𝒗\text{i}n_\mathrm{e}𝒗\text{e}`$ can be written in terms of the densities and velocities of the other charged particles using the neutrality condition, (10), and the definition of the total current,

$$𝒋=\left(n_\mathrm{i}𝒗\text{i}n_\mathrm{e}𝒗\text{e}+n_{\mathrm{g}^+}𝒗_{\mathrm{g}^+}n_\mathrm{g}^{}𝒗_\mathrm{g}^{}+n_{\mathrm{sg}^+}𝒗_{\mathrm{sg}^+}n_{\mathrm{sg}^{}}𝒗_{\mathrm{sg}^{}}\right)e.$$
(16)
Making these substitutions in eq. (15) we find, after a little algebra, that

$`{\displaystyle \frac{\rho _\mathrm{i}}{\tau _{\mathrm{i},\mathrm{n}}^\mathrm{d}}}\left(𝒗\text{n}𝒗\text{i}\right)`$ $`=`$ $`{\displaystyle \frac{\left(\mathbf{}\mathbf{\times }𝑩\right)\mathbf{\times }𝑩}{4\pi }}{\displaystyle \frac{en_{\mathrm{g}^+}}{c}}\left[\left(𝒗_{\mathrm{g}^+}𝒗\text{i}\right)\mathbf{\times }𝑩\right]{\displaystyle \frac{en_{\mathrm{sg}^+}}{c}}\left[\left(𝒗_{\mathrm{sg}^+}𝒗\text{i}\right)\mathbf{\times }𝑩\right]`$ (17)
$`+{\displaystyle \frac{en_\mathrm{g}^{}}{c}}\left[\left(𝒗_{\mathrm{sg}^{}}𝒗\text{i}\right)\mathbf{\times }𝑩\right]+{\displaystyle \frac{en_{\mathrm{sg}^{}}}{c}}\left[\left(𝒗_{\mathrm{sg}^{}}𝒗\text{i}\right)\mathbf{\times }𝑩\right],`$
where we used equation (8) to express the electric field in terms of $`𝒗`$i and Ampère’s Law to write the current in terms of $`\mathbf{}`$$`\mathbf{\times }`$$`𝑩`$. Equation (17) says that the electromagnetic and drag forces on the combined ion+electron fluid are in balance, as they must be when ion and electron inertia are neglected. Equation (17) is solved to find $`𝒗\text{i}(z,t)`$ in terms of the other fluid variables.

The ion energy equation is omitted in favor of the approximation

$$k_\mathrm{B}T_\mathrm{i}k_\mathrm{B}T_\mathrm{n}+\frac{1}{3}m_\mathrm{n}\left|𝒗\text{i}𝒗\text{n}\right|^2$$
(18)
(Chernoff 1987), which says that the ion fluid is in thermal balance with heating and cooling due to ion-neutral scattering.

#### 2.2.5 Small Grains

Our code solves the equations for mass conservation,

$$\frac{\rho _{\mathrm{sg}^{}}}{t}+\frac{}{z}\left(\rho _{\mathrm{sg}^{}}v_{\mathrm{sg}^{}z}\right)=S_{\mathrm{sg}^{}},$$
(19)
$$\frac{\rho _{\mathrm{sg}^0}}{t}+\frac{}{z}\left(\rho _{\mathrm{sg}^0}v_{\mathrm{sg}^0z}\right)=S_{\mathrm{sg}^0},$$
(20)
and

$$\frac{\rho _{\mathrm{sg}^+}}{t}+\frac{}{z}\left(\rho _{\mathrm{sg}^+}v_{\mathrm{sg}^+z}\right)=S_{\mathrm{sg}^+}.$$
(21)
The source terms on the r.h.s. of equations (19)-(21) describe the transfer of grains from one charge state to another due to the capture of ions and electrons, and are computed by integrating the appropriate cross sections (Draine & Sutin 1987) over the distribution of relative velocities for charged particles and drifting grains (see Appendix A). We neglect the inertia of the small grains and compute their velocities algebraically by requiring the total force on each small-grain fluid to vanish. Thus we require

$$0=\rho _{\mathrm{sg}^{}}\left(𝒗_{\mathrm{sg}^{}}𝒗\text{i}\right)\mathbf{\times }𝛀_{\mathrm{sg}^{}}+\frac{\rho _{\mathrm{sg}^{}}}{\tau _{\mathrm{sg}^{},\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗_{\mathrm{sg}^{}}\right)\frac{\rho _{\mathrm{sg}^{}}}{\tau _{\mathrm{sg}^{},\mathrm{i}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^{}}+\frac{\rho _{\mathrm{sg}^0}}{\tau _{\mathrm{sg}^0,\mathrm{e}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^0},$$
(22)
$$0=\frac{\rho _{\mathrm{sg}^0}}{\tau _{\mathrm{sg}^0,\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗_{\mathrm{sg}^0}\right)\frac{\rho _{\mathrm{sg}^0}}{\tau _{\mathrm{sg}^0,\mathrm{i}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^0}\frac{\rho _{\mathrm{sg}^0}}{\tau _{\mathrm{sg}^0,\mathrm{e}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^0}+\frac{\rho _{\mathrm{sg}^{}}}{\tau _{\mathrm{sg}^{},\mathrm{i}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^{}}+\frac{\rho _{\mathrm{sg}^+}}{\tau _{\mathrm{sg}^+,\mathrm{e}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^+},$$
(23)
and

$$0=\rho _{\mathrm{sg}^+}\left(𝒗_{\mathrm{sg}^+}𝒗\text{i}\right)\mathbf{\times }𝛀_{\mathrm{sg}^+}+\frac{\rho _{\mathrm{sg}^+}}{\tau _{\mathrm{sg}^+,\mathrm{n}}^\mathrm{d}}\left(𝒗\text{n}𝒗_{\mathrm{sg}^+}\right)+\frac{\rho _{\mathrm{sg}^0}}{\tau _{\mathrm{sg}^0,\mathrm{i}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^0}\frac{\rho _{\mathrm{sg}^+}}{\tau _{\mathrm{sg}^+,\mathrm{e}}^{\mathrm{cf}}}𝒗_{\mathrm{sg}^+},$$
(24)
where $`𝛀`$f ($`=Q𝑩/m_\mathrm{f}c`$) is the cyclotron frequency for particles of type f with charge $`Q`$ and $`\tau _{\mathrm{i},\mathrm{j}}^{\mathrm{cf}}`$ is the mean time for a small grain to change its affiliation from fluid i to fluid j by capturing an ion or electron. As noted in previous studies (e.g., Pilipp et al. 1987; Ciolek & Mouschovias 1993), the capture of charged particles by the grains transfers mass and momentum between the different grain subfluids. Our code includes these effects via the charge fluctuation terms in the momentum equations and the source terms in the mass equations. The velocity dispersion of the small grains does not affect other properties of the flow; consequently, we do not solve the small-grain energy equation.

#### 2.2.6 Large Grains

We neglect the velocity dispersion of the large grains and omit the energy equation. We include the large-grain inertia and solve the equations for mass conservation,

$$\frac{\rho _\mathrm{g}^{}}{t}+\frac{}{z}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^{}z}\right)=S_\mathrm{g}^{},$$
(25)
$$\frac{\rho _{\mathrm{g}^0}}{t}+\frac{}{z}\left(\rho _{\mathrm{g}^0}v_{\mathrm{g}^0z}\right)=S_{\mathrm{g}^0},$$
(26)
$$\frac{\rho _{\mathrm{g}^+}}{t}+\frac{}{z}\left(\rho _{\mathrm{g}^+}v_{\mathrm{g}^+z}\right)=S_{\mathrm{g}^+},$$
(27)
and momentum conservation,

$`{\displaystyle \frac{}{t}}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^{}y}\right)+{\displaystyle \frac{}{z}}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^{}y}v_{\mathrm{g}^{}z}\right)`$ $`=`$ $`\rho _\mathrm{g}^{}\left(v_{\mathrm{g}^{}z}v_{\mathrm{i}z}\right)\mathrm{\Omega }_\mathrm{g}+{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}y}v_{\mathrm{g}^{}y}\right)`$ (28)
$`{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^{}y}+{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^0y},`$
$`{\displaystyle \frac{}{t}}\left(\rho _{\mathrm{g}^0}v_{\mathrm{g}^0y}\right)+{\displaystyle \frac{}{z}}\left(\rho _{\mathrm{g}^0}v_{\mathrm{g}^0y}v_{\mathrm{g}^0z}\right)`$ $`=`$ $`{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}y}v_{\mathrm{g}^0y}\right)+{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^{}y}{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^0y}`$ (29)
$`{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^0y}+{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^+y},`$
$`{\displaystyle \frac{}{t}}\left(\rho _{\mathrm{g}^+}v_{\mathrm{g}^+y}\right)+{\displaystyle \frac{}{z}}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^+y}v_{\mathrm{g}^+z}\right)`$ $`=`$ $`\rho _{\mathrm{g}^+}\left(v_{\mathrm{g}^+z}v_{\mathrm{i}z}\right)\mathrm{\Omega }_\mathrm{g}+{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}y}v_{\mathrm{g}^+y}\right)`$ (30)
$`+{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^0y}{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^+y},`$
$`{\displaystyle \frac{}{t}}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^{}z}\right)+{\displaystyle \frac{}{z}}\left(\rho _\mathrm{g}^{}v_{\mathrm{g}^{}z}^2\right)`$ $`=`$ $`\rho _\mathrm{g}^{}\left(v_{\mathrm{g}^{}y}v_{\mathrm{i}y}\right)\mathrm{\Omega }_\mathrm{g}+{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}z}v_{\mathrm{g}^{}z}\right)`$ (31)
$`{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^{}z}+{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^0z},`$
$`{\displaystyle \frac{}{t}}\left(\rho _{\mathrm{g}^0}v_{\mathrm{g}^0z}\right)+{\displaystyle \frac{}{z}}\left(\rho _{\mathrm{g}^0}v_{\mathrm{g}^0z}^2\right)`$ $`=`$ $`{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}z}v_{\mathrm{g}^0z}\right)+{\displaystyle \frac{\rho _\mathrm{g}^{}}{\tau _{\mathrm{g}^{},\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^{}z}{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^0z}`$ (32)
$`{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^0z}+{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^+z},`$
and

$`{\displaystyle \frac{}{t}}\left(\rho _{\mathrm{g}^+}v_{\mathrm{g}^+z}\right)+{\displaystyle \frac{}{z}}\left(\rho _{\mathrm{g}^+}v_{\mathrm{g}^+z}^2\right)`$ $`=`$ $`\rho _{\mathrm{g}^+}\left(v_{\mathrm{g}^+y}v_{\mathrm{i}y}\right)\mathrm{\Omega }_\mathrm{g}+{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{n}}^\mathrm{d}}}\left(v_{\mathrm{n}z}v_{\mathrm{g}^+z}\right)`$ (33)
$`+{\displaystyle \frac{\rho _{\mathrm{g}^0}}{\tau _{\mathrm{g}^0,\mathrm{i}}^{\mathrm{cf}}}}v_{\mathrm{g}^0z}{\displaystyle \frac{\rho _{\mathrm{g}^+}}{\tau _{\mathrm{g}^+,\mathrm{e}}^{\mathrm{cf}}}}v_{\mathrm{g}^+z},`$
where $`\mathrm{\Omega }_\mathrm{g}\left|𝛀_\mathrm{g}^{}\right|=\left|𝛀_{\mathrm{g}^+}\right|`$.

#### 2.2.7 Chemistry

The ion and neutral mass density depend on the chemistry, which is calculated in parallel with the hydrodynamics. Our code solves a set of chemical equations for the neutrals,

$$\frac{n_{\mathrm{n}\alpha }}{t}+\frac{}{z}\left(n_{\mathrm{n}\alpha }v_{\mathrm{n}z}\right)=𝒞_{\mathrm{n}\alpha }𝒟_{\mathrm{n}\alpha },$$
(34)
and ions,

$$\frac{n_{\mathrm{i}\alpha }}{t}+\frac{}{z}\left(n_{\mathrm{i}\alpha }v_{\mathrm{i}z}\right)=𝒞_{\mathrm{i}\alpha }𝒟_{\mathrm{i}\alpha },$$
(35)
where $`𝒞_{\mathrm{f}\alpha }`$ and $`𝒟_{\mathrm{f}\alpha }`$ are the net rates per volume at which species $`\alpha `$ in fluid f is created and destroyed, respectively.

### 2.3 Validity and Relevance

#### 2.3.1 Frozen Field Approximation

The frozen field approximation is valid when $`\mathrm{\Gamma }_\mathrm{i}1`$, where

$$\mathrm{\Gamma }_\mathrm{i}=\frac{eB}{m_\mathrm{i}c}36\frac{B_3}{n_8}$$
(36)
is the ion Hall parameter, $`B=B_3`$ mG is the magnetic field, and $`n10^8n_8`$cm-3. The line-of-sight magnetic field components inferred from Zeeman splitting in thermally-excited transitions (summarized by Crutcher 1999) imply that $`\mathrm{\Gamma }_\mathrm{i}1`$ for $`n10^7`$cm-3 but observations of denser regions are rare. Extrapolating the thermal Zeeman data to higher densities using the observed correlation $`Bn^\kappa `$ (with $`\kappa 0.5`$, see Crutcher 1999; Basu 2000) would imply that the frozen-field approximation is valid for $`n10^9`$cm-3. The extrapolated field values are consistent with the fields estimated in H2O masers with $`n10^9`$cm-3 but the field strengths inferred from maser observations are controversial (Sarma, Troland & Romney 2001 and references therein).

#### 2.3.2 Inertia of Ions and Small Grains

Neglecting the inertia of ions and small grains is an excellent approximation for steady shocks. Let $`L_{\mathrm{drag}}(\mathrm{sg})`$ be the length scale for gas drag to decelerate the small grains, $`L_{\mathrm{neut}}`$ the scale for compression of the neutrals, and $`L_{\mathrm{field}}`$ the scale for compression of the magnetic field. Figure 2 compares these and other length scales for a shock with speed $`v_\mathrm{s}=10\text{km s}\text{-1}`$ and various preshock densities. The dynamical scales were estimated using dimensional analysis of the equations of motion for a steady shock (DRD83; see eq. –):

$$L_{\mathrm{neut}}=\frac{\rho _{\mathrm{n0}}v_\mathrm{s}^2}{F_\mathrm{n}},$$
(37)
and

$$L_{\mathrm{field}}=\frac{B_0^2}{8\pi F_\mathrm{i}},$$
(38)
where the subscript “0” denotes preshock quantities. We took the preshock magnetic field to be

$$B_0=b\left(n_0/\text{cm}\text{-3}\right)^{1/2}\mu \mathrm{G}$$
(39)
with $`b=0.3`$. We included the drag due to electron-neutral and ion-neutral scattering (DRD83, eq. -) as well as the gas drag on small and large grains (D80, eq. ). We computed $`L_{\mathrm{field}}`$ at a typical point in the magnetic precursor where the neutral particles are unperturbed and the magnetic field has been compressed by a factor of 2 (see D80, eq. ). For $`L_{\mathrm{neut}}`$ we took a representative downstream point, with unperturbed neutrals and the magnetic field compressed by a factor of $`M_{\mathrm{nA}}`$, the neutral Alfvén Mach number (see Kaufman & Neufeld 1996a and references therein). Preshock chemical and grain abundances were taken from Figure 1 and the drift velocities of large and small grains were computed by neglecting their inertia for expediency. We set

$$L_{\mathrm{drag}}(\mathrm{sg})=v_{\mathrm{sg}^{}z}\tau _{\mathrm{sg}^{},\mathrm{n}}^\mathrm{d},$$
(40)
with the grain velocity evaluated in a frame where the shock is steady and the drag time computed from the velocity-dependent drag force on a drifting sphere (Draine & Salpeter 1979). ¿From Figure 2 and analogous calculations for other shock speeds we conclude that the inertia of small grains (and hence the ions) can be neglected in steady, nondissociative shocks with $`v_\mathrm{s}<30\text{km s}\text{-1}`$.

In time-dependent shocks, the inertia of light particles is a subtler issue. The formation of a magnetic precursor depends on the speed at which compressive disturbances are communicated to the charged fluid, and this “signal speed” is modified when the inertia of light particles is neglected. As a simple illustration, consider the propagation of small-amplitude waves in a plasma of ions, electrons, and neutrals. In Figure 3 we plot the phase velocity of waves propagating perpendicular to $`𝑩`$, as predicted by the exact dispersion relation for a cold plasma (solid curves) and the dispersion relation with ion inertia neglected (open circles). The exact dispersion relation predicts the well-known result (Kulsrud & Pearce 1969) that compressive waves propagate in the ion-electron fluid for wavelengths

$$\lambda <\lambda _1=4\pi v_{\mathrm{ims}}\tau _{\mathrm{i},\mathrm{n}}^\mathrm{d}$$
(41)
and in the bulk plasma for

$$\lambda >\lambda _2=\lambda _1\sqrt{\rho _\mathrm{n}/4\rho _\mathrm{i}},$$
(42)
where $`v_{\mathrm{ims}}=B/\sqrt{4\pi \rho _\mathrm{i}}`$ is the ion magnetosound speed for a cold plasma. When the ion inertia is neglected, the short-wavelength mode disappears and the fastest “signal speed” for propagating waves drops from $`v_{\mathrm{ims}}`$ to the bulk magnetosound speed, $`v_{\mathrm{ms}}v_{\mathrm{ims}}\sqrt{\rho _\mathrm{i}/\rho _\mathrm{n}}`$. This seems to suggest that the “inertialess equations” would fail to admit C type solutions for $`v_\mathrm{s}>v_{\mathrm{ms}}`$. However Smith and Mac Low (1997) have shown, by direct numerical integration of the inertialess equations, that a J shock with $`v_\mathrm{s}>v_{\mathrm{ms}}`$ can evolve into a C shock. We describe a similar result in §3. The development of a magnetic precursor in both examples is initiated by rapid ion motions with $`v_\mathrm{i}v_{\mathrm{ims}}`$. Evidently the formation of a C shock does not require propagating waves with phase velocities larger than the shock speed. We speculate that the ions are accelerated by evanescent waves in the “forbidden zone” $`\lambda _1<\lambda <\lambda _2`$, which have long damping times (Fig. 4). If correct, our conjecture may have interesting implications for maximum speed of a C shock. This and related issues will be discussed elsewhere (Roberge & Ciolek, in preparation).

#### 2.3.3 Relevance of Large-Grain Inertia and Charge Fluctuations

The drag length for large grains \[$`L_{\mathrm{drag}}`$(g)\] is shown in Figure 2. For $`n10^6`$cm-3, the large grains are partially coupled to $`𝑩`$ and may affect the formation of C shocks by “loading” the field lines (Draine & McKee 1993; see §1). Over this density range the drag length is comparable to $`L_{\mathrm{field}}`$ (Fig. 2b), so including the large-grain inertia is necessary. Length scales for the capture of charged particles by small (sg0+e, sg-+i) and large (g0+e, g-+i) grains are also plotted in Figure 2. The charge fluctuation times were computed using abundances from Figure 1 and rate coefficients from Appendix A. Because the shock structure is sensitive to changes in the ionization (Flower, Pineau des Forets & Hartquist 1985; Pilipp et al. 1990), it is necessary to include the charge fluctuations of small grains for $`n10^9`$cm-3 (Fig. 2a). Figure 2b shows that large grains undergo many charge fluctuations as the magnetic field is compressed, partially coupling the neutral grains to the field lines. Consequently the charge fluctuations of large grains are also important.

## 3 Solutions

### 3.1 Initial and Boundary Conditions

We model two types of problem with our code. In the first type (“free propagation problems”), the initial conditions are some disturbance, linear or nonlinear, specified on a finite computational domain, $`0<z<Z`$. For free propagation problems we impose “non-reflecting boundary conditions” (discussed below) at $`z=0`$ and $`z=Z`$. These boundary conditions allow the disturbance to propagate off the domain as if the boundaries were at infinity. In the other type of problem (“inflow problems”), the initial conditions are a uniform, static “undisturbed state” for $`0<z<Z`$ and a steady flow is incident at the boundary $`z=0`$. For inflow problems we specify the flow variables at $`z=0`$ and impose nonreflecting boundary conditions at $`z=Z`$.

For inflow problems, the undisturbed state is specified by choosing the density $`n_0`$ and magnetic field strength, $`B_0`$. In the exemplary calculations presented in this paper, the temperature of the undisturbed gas is not of interest. We assign the initial neutral, ion, and electron temperatures to reasonable but somewhat arbitrary values instead of, say, requiring the fluids to be in thermal balance. Once the gas density and temperature are specified, the initial densities $`n_{\mathrm{f}\alpha 0}`$ of the various chemical species are determined by solving the chemical rate equations for a static system (see §2.2.7). To determine the abundances uniquely we must specify the fractional abundance of heavy nuclei ($`x_\mathrm{a}`$), and the initial abundances of large ($`x_{\mathrm{g},0}`$), and small ($`x_{\mathrm{sg},0}`$) grains. We typically set $`x_\mathrm{a}10^6`$. The initial grain abundances are specified in practice by choosing the grain radii, $`a_\mathrm{g}`$ and $`a_{\mathrm{sg}}`$, the densities $`\rho _\mathrm{g}`$ and $`\rho _{\mathrm{sg}}`$ of the grain materials, and the initial dust-to-gas mass ratios, $`\chi _{\mathrm{g},0}`$ and $`\chi _{\mathrm{sg},0}`$. For diffuse regions of molecular clouds or protostellar cores, we typically set $`a_\mathrm{g}=0.1`$ µm $`\rho _\mathrm{g}=3.3`$g cm-3, and $`\chi _{\mathrm{g},0}=10^2`$, yielding $`x_\mathrm{g}=2\times 10^{12}`$. For the small grains we usually take $`a_{\mathrm{sg}}=4`$ Å, $`\rho _{\mathrm{sg}}=2.2`$g cm-3, and $`\chi _{\mathrm{sg},0}=5\times 10^5`$, yielding an abundance, $`x_{\mathrm{sg},0}=2\times 10^7`$, consistent with the constraints on small grains in diffuse gas (Puget & Leger 1989; Weingartner & Draine 2001; Li & Draine 2001). For models of shocks in dense gas, we adjust the grain radii and dust-to-gas ratios to account approximately for the possible depletion of small grains, coagulation of large grains, (e.g., Ossenkopf 1993), and depletion of large grains due to ambipolar diffusion during core formation (Ciolek & Mouschovias 1994, 1996).

The actual problem we solve is a dimensionless one. The governing equations (§ 2) are put into nondimensional form by using as a unit of distance some specified length $`L`$, with $`Z=L`$, $``$ a constant. Velocities are normalized to $`C_{\mathrm{n0}}`$, the isothermal sound speed in the undisturbed neutral gas. (For free propagation problems, “undisturbed gas” means the gas outside the initial disturbance.) The unit of time is therefore $`L/C_{\mathrm{n0}}`$. Temperatures and densities are normalized to the temperature, $`T_{\mathrm{n0}}`$, and mass density, $`\rho _{\mathrm{n0}}`$, of the undisturbed neutral gas, and the unit of magnetic field strength is taken to be $`B_0`$. A consequence of our selection of units is that the dimensionless Alfvén speed of the undisturbed gas, $`v_{\mathrm{A},\mathrm{n0}}/C_{\mathrm{n0}}`$ \[$`=B_0/(4\pi \rho _{\mathrm{n0}})^{1/2}C_{\mathrm{n0}}`$\], appears as a characteristic parameter in the problem.

The use of non-reflecting boundary conditions (also referred to as “transmissive” or “radiation” boundary conditions) allows the flow to pass unimpeded through the computational boundaries, and prevents the reflection of spurious waves from the boundaries back into the computational domain. A comprehensive review of non-reflecting boundary conditions in numerical hydrodynamics is presented by Givoli (1991). We use the following formulation, due to Thompson (1987) and Vanajakshi, Thompson, & Black (1989). Examination of the governing equations (§2) reveals that they are all conservative transport equations of the form

$$\frac{𝓤}{t}+\frac{𝓕}{z}+𝓒=0,$$
(43)
where

$$𝓤=\left[\begin{array}{cc}\rho _\mathrm{n}& \\ \rho _\mathrm{n}v_{\mathrm{n}y}& \\ \rho _\mathrm{n}v_{\mathrm{n}z}& \\ \frac{1}{2}\rho _\mathrm{n}v_\mathrm{n}^2+E_\mathrm{n}& \\ \mathrm{}& \end{array}\right]$$
(44)
is the (column) vector of conserved variables, $`𝓕`$ is the corresponding vector flux of conserved quantities, and $`𝓒`$ is the vector of source terms. The system (43) is transformed into the quasilinear form

$$\frac{𝓤^{\mathbf{}}}{t}+𝗔\mathbf{}\frac{𝓤^{\mathbf{}}}{z}+𝓒^{\mathbf{}}=0,$$
(45)
where

$$𝓤^{\mathbf{}}=\left[\begin{array}{c}\rho _\mathrm{n}\\ v_{\mathrm{n}y}\\ v_{\mathrm{n}z}\\ E_\mathrm{n}/\rho _\mathrm{n}\\ \mathrm{}\end{array}\right],$$
(46)
is the vector of primitive variables and the Jacobian matrix, $`𝗔(𝓤^{\mathbf{}})`$, is readily determined. The eigenvalues $`\left\{\xi _k\right\}`$ of $`𝗔`$ are the characteristic speeds of the flow. Non-reflecting boundary conditions are implemented by imposing the set of equations

$$𝒍_k\mathbf{}\frac{𝓤^{\mathbf{}}}{t}+_k+𝒍_k\mathbf{}𝓒^{\mathbf{}}=0,$$
(47)
at the appropriate boundary, where $`𝒍_k`$ is the left (row) eigenvector associated with the eigenvalue $`\xi _k`$ and

$$_k=\{\begin{array}{c}\xi _k𝒍_k\mathbf{}𝓤^{\mathbf{}}/z\text{for outgoing waves,}\\ \\ 0\text{for incoming waves.}\end{array}$$
(48)
As demonstrated by Thompson (1987) and Vanajakshi et al. (1989), the conditions (48) for outgoing waves ($`\xi _k<0`$ at $`z=0`$, $`\xi _k>0`$ at $`z=Z`$) and incoming waves ($`\xi _k>0`$ at $`z=0`$, $`\xi _k<0`$ at $`z=Z`$) eliminate spurious reflections and produce accurate numerical solutions at the artificial boundaries.

### 3.2 Numerical Method

The dimensionless governing equations with nonreflective and/or inflow boundary conditions comprise a set of coupled, nonlinear, partial differential equations (PDEs). We solve the PDEs numerically using a technique known as the “method of lines” (e.g., Liskovets 1965; Schiesser 1991). The computational domain is divided into $`N`$ cells of uniform width $`\mathrm{\Delta }z=Z/N`$ with cell boundaries at $`\{z_j=j\mathrm{\Delta }z,j=0,1,\mathrm{},N\}`$. Mean values of the variables are assigned to the center of each cell, i.e., to the points

$$\overline{z}_j\frac{(z_j+z_{j1})}{2}=z_{j1}+\frac{1}{2}\mathrm{\Delta }z=(j\frac{1}{2})\mathrm{\Delta }z,$$
(49)
for $`j1`$. The mean value of variable $`f`$ in cell $`j`$ will be denoted by $`\overline{f}_jf(\overline{z}_j)`$. In the method of lines, spatial derivatives are approximated using finite differences, thereby converting the PDEs into a system of ordinary differential equations (ODEs) of the form

$$\frac{d𝒚}{dt}=𝒇[𝒚(t),t],$$
(50)
where $`d𝒚/dt`$, $`𝒚`$, and $`𝒇`$ are vectors of length $`WN`$, and $`W`$ is the number of variables defined at each cell. The ODEs are solved using standard stiff ODE solving subroutines readily available in the public domain.

For cells inside the computational domain, spatial derivatives are approximated with the fourth-order differencing scheme

$$\left[\frac{f}{z}\right]_j\frac{(1/12)(\overline{f}_{j+2}\overline{f}_{j2})+(2/3)(\overline{f}_{j+1}\overline{f}_{j1})}{\mathrm{\Delta }z}.$$
(51)
For cells adjacent to the boundaries, we use the three-point scheme

$`\left[{\displaystyle \frac{f}{z}}\right]_1`$ $``$ $`{\displaystyle \frac{(1/3)f_2+f_1(4/3)f_{b0}}{\mathrm{\Delta }z}},`$ (52)
$`\left[{\displaystyle \frac{f}{z}}\right]_N`$ $``$ $`{\displaystyle \frac{(4/3)f_{bN}f_N(1/3)f_{N1}}{\mathrm{\Delta }z}},`$ (53)
where $`f_{b0}`$ and $`f_{bN}`$ refer to the values of variables at the boundaries $`z=0`$ and $`z_N=Z`$, respectively. Thus, derivatives at the boundaries depend only on variables inside the computational domain. They are accurate to $`𝒪(\mathrm{\Delta }z)^2`$, while the interior derivatives (eq. ) are accurate to $`𝒪(\mathrm{\Delta }z^4)`$.

Advection equations, of the form $`f/t=(fv)/z`$, are discretized using the monotonic, upwind algorithm of van Leer (1979). Following this algorithm, advection terms are approximated as

$$\left[\frac{(fv)}{z}\right]_j\frac{(f__Rv__Rf__Lv__L)}{\mathrm{\Delta }z},$$
(54)
where the flux $`f__R`$$`v__R`$ at the right cell face is approximated by setting $`v__R=(\overline{v}_j+\overline{v}_{j+1})/2`$, and

$$f__R=\{\begin{array}{c}\overline{f}_j+\frac{1}{2}(\overline{\mathrm{\Delta }f})_j\mathrm{if}v__R0,\\ \\ \overline{f}_{j+1}\frac{1}{2}(\overline{\mathrm{\Delta }f})_{j+1}\mathrm{if}v__R<0,\end{array}$$
(55)
and the flux $`f__L`$$`v__L`$ at the left cell face is approximated by setting $`v__L=(\overline{v}_{j1}+\overline{v}_j)/2`$ and

$$f__L=\{\begin{array}{c}\overline{f}_{j1}+\frac{1}{2}(\overline{\mathrm{\Delta }f})_{j1}\mathrm{if}v__L0,\\ \\ \overline{f}_j\frac{1}{2}(\overline{\mathrm{\Delta }f})_j\mathrm{if}v__L<0.\end{array}$$
(56)
Here we have introduced the “monotonic derivative,”

$$(\overline{\mathrm{\Delta }f})_j=\{\begin{array}{cc}\frac{(\overline{f}_{j+1}\overline{f}_j)(\overline{f}_j\overline{f}_{j1})}{\overline{f}_{j+1}\overline{f}_{j1}}\mathrm{if}(\overline{f}_{j+1}\overline{f}_j)(\overline{f}_j\overline{f}_{j1})>0,& \\ & \\ 0\mathrm{otherwise}.& \end{array}$$
(57)
For the boundary cell adjacent to $`z=0`$ we set $`v__L=v_{b0}`$. If $`v__L0`$, $`f_L=f_{b0}`$, otherwise $`f_L`$ is given by equation (56). Similarly, for the cell adjoining $`z=Z`$, $`v_R=v_{bN}`$. If $`v_R0`$, equation (55) is used for $`f_L`$; if $`v_R<0`$, $`f_R=f_{bN}`$ instead.

### 3.3 Benchmark Tests

To establish the accuracy of our numerical code, in this section we present the results of various test runs and compare them with known solutions. In all of the tests presented here, we set $`E_\mathrm{n}=P_\mathrm{n}/(\gamma 1)`$, where $`\gamma `$ is the adiabatic index of the neutral gas, assumed constant. In models that include cooling, we use the analytical cooling function $`\mathrm{\Lambda }_\mathrm{n}`$ defined by Roberge & Draine (1990, eq. \[3.18\]).

#### 3.3.1 Shock Tube

A standard test for numerical codes is the well-known shock tube problem (e.g., Sod 1978). In the absence of magnetic fields, the solution for this type of Riemann problem can be found exactly using analytic methods (Toro 1997 and references therein). The gas is adiabatic in the shock tube test, with $`\gamma =5/3`$. The magnetic field strength is made vanishingly small; in this limit, the ions, electrons, and grains move with the neutrals (due to collisional forces) and do not appreciably affect the flow, because the density of charged particles is much smaller than the neutral gas density.

The shock tube test falls into the category of “free propagation” problems with nonreflecting boundary conditions at both edges of the spatial domain (see §3.1). The neutral gas is initially at rest. The initial density and temperature profiles are constant on both halves of the domain with a discontinuity at $`z=Z/2`$. In our test we chose $`\rho _\mathrm{n}(z)=12\rho _{\mathrm{n0}}`$ and $`T_\mathrm{n}(z)=1.875T_{\mathrm{n0}}`$ for $`z<Z/2`$ and $`\rho _\mathrm{n}(z)=\rho _{\mathrm{n0}}`$ and $`T_\mathrm{n}(z)=T_{\mathrm{n0}}`$ for $`z>Z/2`$, where $`\rho _{\mathrm{n0}}`$ and $`T_{\mathrm{n0}}`$ are scaling values that merely define the dimensionless units (see §3.1).

For this test we turned off the source terms in the continuity equations for all fluid species. Because the gas flow is taken to be adiabatic, the cooling function for the neutral gas is set to zero.

The results of the shock tube test for a dimensionless time $`t_1>0`$ are displayed in Figures $`4a`$-$`4d`$ (circles). This particular model used 400 mesh points; for visual clarity in each plot we display only every eighth data point. Figure $`4a`$ shows $`\rho _\mathrm{n}(z,t_1)/\rho _{\mathrm{n0}}`$ as a function of $`z/Z`$. The ratio $`v_{\mathrm{n}z}(z,t_1)/C_{\mathrm{n0}}`$ is exhibited in Figure $`4b`$. The pressure $`P_\mathrm{n}(z,t_1)`$ relative to $`P_{\mathrm{n},0}`$, and $`T_\mathrm{n}(z,t_1)/T_{\mathrm{n},0}`$ are shown in Figures $`4c`$ and $`4d`$, respectively. For comparison, also displayed in these figures are the exact Riemann solutions (solid curves). Examination reveals that our multifluid code effectively reproduces the full structure of the shock tube problem, including the shock wave, rarefaction wave, and contact discontinuity which are present at this time. The r.m.s. relative errors of the numerically-computed variables are typically $`10^2`$ — except at the oscillations near the base of the rarefaction wave and at the contact discontinuity.

Figures $`5a`$-$`d`$ present numerical results for the same model run (circles) at a later time $`t_2`$, which are again compared to the exact solution (solid curves). By this time, both the forward-moving shock wave and contact discontinuity, and also the backward propagating rarefaction wave, all present in Figures $`4a`$-$`d`$, have moved “off stage”, that is, past the boundaries at $`z=Z`$ and $`0`$. We note that passage of the forward waves through $`z=Z`$ does not produce any substantial numerical reflections at the boundary located there. Nor are there any spurious reflections generated by the backward rarefaction wave incident upon the opposite boundary at $`z=0`$. This illustrates the efficacy of the non-reflecting boundary conditions (§ 3.1).

#### 3.3.2 Nonmagnetic Shock with Inflow Conditions

The next test we present is an inflow problem. We again consider nonmagnetic, adiabatic ($`\gamma =5/3`$) flow. Initially the neutral gas within the entire region \[0, $`Z`$\] is at rest with uniform density $`\rho _{\mathrm{n0}}`$ and temperature $`T_{\mathrm{n0}}`$. At $`t=0`$ a high-velocity shock with $`v_{\mathrm{n}z}=50C_{\mathrm{n0}}`$ is incident on the boundary at $`z=0`$; this corresponds to a Mach number $`M=v_\mathrm{s}/\gamma ^{1/2}C_{\mathrm{n0}}=51.7`$, where $`v_\mathrm{s}`$ is the shock speed. The inflowing shocked material has density $`\rho _\mathrm{n}=3.99\rho _{\mathrm{n0}}`$, and temperature $`T_\mathrm{n}=835T_{\mathrm{n0}}`$.

Figures $`6a`$-$`d`$ depict the profiles of the density, velocity, pressure, and temperature of the neutrals, normalized to the uniform initial state, at a time $`t>0`$. As before, the circles represent the solution from our multifluid numerical code. The solid line is the exact Riemann solution. Comparison indicates that there is good agreement between the numerical and exact solutions.

#### 3.3.3 Counterpropagating Rarefaction Waves in a Magnetic, Flux-Frozen Fluid

This test is used to verify that our numerical code has the correct MHD behavior in the limit of flux-freezing, that is, when the magnetic field is frozen into the bulk fluid. Effective freezing of the magnetic flux into the neutrals can occur when the characteristic minimum flow time $`\tau _{\mathrm{flow},\mathrm{min}}`$ ($`=\mathrm{\Delta }z/v_{\mathrm{n}z}`$) is $`\tau _{\mathrm{n},\mathrm{tot}}^\mathrm{d}`$, the drag time for neutrals colliding with all charged particles attached to the magnetic field. In this situation, information about the magnetic field is quickly communicated to the neutrals, through collisions with the plasma, on a time scale much smaller than the natural time scale associated with the flow. When this happens, the magnetic field, plasma, and the neutrals can rapidly respond to one another, and thus move together as a single combined fluid. In our code, the condition $`\tau _{\mathrm{flow},\mathrm{min}}/\tau _{\mathrm{n},\mathrm{tot}}^\mathrm{d}1`$ is readily satisfied by using a state with a large enough abundance of charged particles (since $`\tau _{\mathrm{n},\mathrm{f}}x_\mathrm{f}^1`$; e.g., see eqs. \[2c\] and \[10a\] of Ciolek & Mouschovias 1993), and/or choosing a sufficiently large value $`Z`$ for the size of the physical region being modeled (recall that $`\mathrm{\Delta }z=Z/N`$).

The model we choose for this test is the counterpropagating MHD rarefaction flow described in § 4 of Ryu & Jones (1995). This test consists of two oppositely-directed, adiabatic rarefaction waves in a magnetized fluid. The initial state has a nonzero, discontinuous velocity profile with $`v_{\mathrm{n}z}(z<Z/2)=C_{\mathrm{n0}}`$, and $`v_{\mathrm{n}z}(z>Z/2)=C_{\mathrm{n0}}`$. However, all of the other physical quantities are uniform in the domain \[$`0,Z`$\]: $`\rho _\mathrm{n}(z)=\rho _{\mathrm{n0}}`$, $`T_\mathrm{n}(z)=T_{\mathrm{n0}}`$, and $`B_x(z)=B_0`$. The Alfvén speed in the undisturbed neutrals, $`v_{\mathrm{A},\mathrm{n0}}`$, is equal to $`C_{\mathrm{n0}}`$ for this model. Creation/destruction of charged species was not considered in the analysis of Ryu & Jones. Hence, to be able to make an accurate comparison to with their published model results, in this test we switched off the source and sink terms that appear in the continuity equation for each charged species.

Figures $`7a`$-$`d`$ exhibit the state of this test model at a time $`t>0`$, when the forward- and backward-traveling rarefaction waves have opened a cavity in the density and magnetic field where mass and magnetic flux have flowed away from the center. The curves with circles show (a) the density, (b) velocity, (c) pressure, and (d) the magnetic field profiles generated by our numerical code at this particular time. For comparison, also shown in each plot is the solution (solid lines) of Ryu & Jones (1995, their Fig. 3b), which they obtained both numerically and by using a specially-developed MHD Riemann solver. These figures show that our multifluid code reproduces the limiting ideal MHD flow. That the flow in this model is indeed effectively flux-frozen can be seen by comparing Figures $`7a`$ and $`7d`$, which show that the mass and magnetic flux have behaved in the same way, preserving the local mass-to-flux ratio $`(\rho _\mathrm{n}/B_x)/(\rho _{\mathrm{n0}}/B_0)=1`$ throughout the evolution. There is, however, a numerical overshoot at the half-way point, $`z=Z/2`$, which is the location of the contact discontinuity within the initial state. A similar overshoot was also found to occur in the numerical simulation of Ryu & Jones.

Finally, we note that there are no spurious oscillations at either boundary in Figures $`7a`$-$`d`$. This confirms that the non-reflecting condition imposed at the artificial boundaries in our code is also properly dealing with the passage of MHD flows.

#### 3.3.4 Chernoff’s Solution for Adiabatic C Shocks

We also test the effectiveness of our code in dealing with realistic, non-ideal MHD situations, such as when the magnetic field is not necessarily frozen into the neutrals, due to relatively small abundances of charged particles. To examine this regime, we modeled the formation of a C shock in a weakly-ionized plasma. As noted earlier, this problem for three-component plasmas (neutrals, ions, and electrons) has previously been considered in great detail for both steady (e.g., D80; DRD83) and also nonsteady (e.g., Tóth 1994; Smith & Mac Low 1997) flows.

For the case of time-independent flow, a detailed phase plane analysis of ‘energy-conserving’ (i.e., adiabatic) MHD shocks was presented by Chernoff (1987) for three-component plasmas. We use his results as the benchmark for this particular test.

To compare our results to those of Chernoff, we used a model containing relative abundances of dust grains that are orders of magnitude smaller than the “typical values” normally adopted (§ 3.1). This reduced the effects of the dust grains, making our test model effectively a three-component plasma. We also neglected the source terms in the continuity equations, since chemistry was also neglected by Chernoff. The source terms for momentum transfer due to collisions between ions and neutrals were taken from DRD83 and Ciolek & Mouschovias (1993).

For initial conditions, we used an adiabatic inflow model with a Mach number $`M=9.75`$ incident on a uniform initial state. The values for the density, temperature, and magnetic field strength of the inflow at $`z=0`$ were determined by the jump conditions for a perpendicular, ideal (i.e., frozen-in magnetic flux) MHD J shock (e.g., Priest 1984). The Alfvén speed in the neutrals for the pre-shocked initial state is $`6.17C_{\mathrm{n0}}`$.

As the flow expands into the computational region for $`t>0`$, the initial J shock rapidly evolves into a C shock. The transient, short-term evolution during the transition from a J-type to C-type shock is similar to that which was reported by Tóth (1994), and Smith & Mac Low (1997), including an early, rapid diffusion of the ions and magnetic field outward from the initially constricting J shock (with $`v_{\mathrm{i}z}v_{\mathrm{n}z}`$ in this phase of the evolution — for more discussion of transient phenomena in a representative model, see § 3.3.5 below).

A steady C shock is soon established, with a velocity $`v_\mathrm{s}=12.6C_{\mathrm{n0}}`$. The Alfvén Mach number $`M_{\mathrm{nA}}`$ ($`=v_\mathrm{s}/v_{\mathrm{A},\mathrm{n0}}`$) for this shock is 2.04; this result is consistent with the analysis of Chernoff (1987), who found that under these conditions (‘energy conserving’ flow), C shocks should exist for $`M_{\mathrm{nA}}<2.8`$. Figure $`8a`$ shows the profile of the neutral velocity (circles) and the ions (filled diamonds) at a time after the steady state is attained. The smooth, supersonic-to-supersonic transition in the flow, characteristic of a C shock, is apparent in this figure, and so is the magnetic precursor in the ion fluid that precedes the neutral flow at the head of the shock.

The earlier studies of MHD shocks in partially ionized fluids were usually carried out in the frame of reference moving with the shock. Figure $`8b`$ displays the shock-frame velocities of the neutrals ($`v_{\mathrm{n}z,\mathrm{shk}}=v_{\mathrm{n}z}v_\mathrm{s}`$; circles), and ions ($`v_{\mathrm{i}z,\mathrm{shk}}=v_{\mathrm{i}z}v_\mathrm{s}`$; filled diamonds), for the same time as in Figure $`8a`$. As in previous examples, the velocities are normalized to $`C_{\mathrm{n0}}`$.

We may also quantitatively compare our numerical results to the detailed predictions of the phase-plane analysis of shocks in three-component plasmas by Chernoff (1987). Figure $`8c`$ displays the values (circles) of the phase variable $`r=|v_{\mathrm{n}z,\mathrm{shk}}|/v_\mathrm{s}`$ as a function of $`q=|v_{\mathrm{i}z,\mathrm{shk}}|/v_\mathrm{s}`$, for the numerical results presented in Figure $`8b`$. The dashed lines in Figure $`8c`$ are the analytical values for $`r(q)`$ calculated by using equation (35) of Chernoff (1987), for a shock with an Alfvén Mach number $`M_{\mathrm{nA}}=2.04`$. The essential equality of these two curves in this figure indicates that the results of our multifluid numerical code are in accord with the analytical solution.

#### 3.3.5 Illustrative Solution with Grain Dynamics

We now present an exemplary model that demonstrates some of the unique features of our time-dependent, multifluid numerical MHD code. As described above in § 2, this includes the ability to follow non-steady, transient phenomena, as well being able to follow simultaneously the dynamics of several different fluid species.

The illustrative model we present is an inflow J shock. The flow for $`z<0`$, $`t<0`$, is assumed to have a frozen-in magnetic field (so that the plasma, magnetic field, and neutrals move together), with the velocity and other quantities at $`z=0`$ determined by the jump conditions for a perpendicular shock. The resulting flow inside the computational domain $`[0,Z]`$ is then determined by solving the full set of governing equations, including the continuity equations for all 9 fluids (with source terms), and also the analytical cooling function $`\mathrm{\Lambda }_\mathrm{n}`$ (§ 3.3).

For specificity, the pre-shock neutral gas has a uniform number density $`n_{\mathrm{n0}}=\rho _{\mathrm{n0}}/m_\mathrm{n}=1.03\times 10^8\text{cm}\text{-3}`$, and temperature $`T_{\mathrm{n0}}=10`$ K. The magnetic field of the uniform initial state in this model is $`B_0=600\mu \mathrm{G}`$. The initial abundances of the various particle species are $`x_{\mathrm{i0}}=2.04\times 10^{10}`$, $`x_{\mathrm{e0}}=3.26\times 10^{12}`$, $`x_{\mathrm{g}^{},0}=4.13\times 10^{13}`$, $`x_{\mathrm{g}^0,0}=2.34\times 10^{12}`$, $`x_{\mathrm{g}^+,0}=3.42\times 10^{13}`$, $`x_{\mathrm{sg}^{},0}=9.94\times 10^{10}`$, $`x_{\mathrm{sg}^0,0}=9.82\times 10^8`$, and $`x_{\mathrm{sg}^+,0}=7.93\times 10^{10}`$. The density and magnetic field strength of this state are such that the Hall parameter for the small charged grains is $`\mathrm{\Gamma }_{\mathrm{sg}}=\mathrm{\Omega }_{\mathrm{sg}}\tau _{\mathrm{sg},\mathrm{n}}^\mathrm{d}=25.41`$. This means that the magnetic field in this model is effectively frozen into the plasma of ions and also the small charged grains. The Hall parameter of the large grains is $`\mathrm{\Gamma }_\mathrm{g}=\mathrm{\Omega }_\mathrm{g}\tau _{\mathrm{g},\mathrm{n}}^\mathrm{d}=2.54\times 10^31`$, so they are not attached to the magnetic field; this means that their dynamics are primarily governed by collisions with the neutrals. Finally, the inflowing J shock has $`M=40`$, and $`M_{\mathrm{nA}}=115`$, corresponding to a speed $`v_\mathrm{s}=9.8\text{km s}\text{-1}`$.

We first examine the model flow on a length scale $`10^{14}`$ cm. On this length scale, and for these parameters, the shock travels as a typical J shock. Figures $`9a`$-$`9c`$ shows the profile of various physical quantities at $`t=5.5\times 10^8`$s. Over these length and time scales, the coupling between the neutrals and the other species is rapid enough, compared to the characteristic flow time of the shock, that the overall motion of the fluids is essentially as a single, combined fluid.

In contrast, we now consider the same model on the much smaller scale $`10^{11}`$ cm. This scale is now small enough that the effect of drift of the magnetic field and the plasma with respect to the neutrals is discernible. So too, is the drift or lag of the large dust grains with respect to the neutrals because of inertial effects. This is demonstrated in Figures $`10a`$ and $`10c`$, which show the velocities of the neutrals (circles), ions and small charged grains (filled diamonds), and the large dust grains (crosses) at times $`t=1.32\times 10^3`$s and $`t=3.68\times 10^3`$s, respectively. As can be seen from these figures, these times are so early on that very little of the inflowing neutral gas has penetrated into the physical/computational region. However, there has been significant readjustment of the ions and the small charged grains, due to plasma drift. The diffusion of the charged particles is now resolved; the plasma and magnetic field have become ‘unlatched’ from the inflowing neutrals, as evident by the higher-velocity precursor in these figures.

Figure $`10b`$ shows the profiles of the neutral density (circles) and the magnetic field (filled diamonds), normalized to the uniform reference values, at $`t=1.32\times 10^3`$s; Figure $`10d`$ shows the same thing at $`t=3.69\times 10^3`$s. Again, there is little penetration of the neutrals into the computational region during this early phase of the evolution. Because the magnetic gradient is much steeper at earlier times (being greatest at time $`t=0`$, when the gradient is initially the discontinuity of the J shock), the velocity of the ions and the small charged grains is much larger at earlier times. This accelerates the plasma ahead of the shock, and gives rise to the precursor. The driving gradient decreases with time as the magnetic field spreads out into the lower-density, lower-magnetic field strength undisturbed region ahead of the shock, and the velocity of the plasma correspondingly decreases from its maximum. From Figures $`10ad`$ the effective signal or disturbance speed of the precursor at these early times is found to be $`10^2\text{km s}\text{-1}`$. This value is intermediate between the ion magnetosound speed and the magnetosound speed for the small charged grains, which are $`1800\text{km s}\text{-1}`$ and $`36\text{km s}\text{-1}`$, respectively.

Notably lagging behind the neutrals in Figures $`10a`$ and $`10c`$ are the large dust grains; only a few non-zero velocity data points for the large grains (crosses) appear in the lower left-hand corners of these figures. This is because of the inertia of the large grains, which is noticable on these scales. The distance over which the large grains are accelerated due to collisions with the neutrals is $`L_{\mathrm{drag}}(\mathrm{g})=(v_{\mathrm{n}z}v_{\mathrm{g},z})\tau _{\mathrm{g},\mathrm{n}}^\mathrm{d}=\rho _\mathrm{s}a_\mathrm{g}/m_\mathrm{n}n_\mathrm{n}`$, where we have used the previously defined collision time scale, the relation $`m_\mathrm{g}=(4\pi /3)\rho _\mathrm{s}a_\mathrm{g}^3`$, and the fact that the relative drift speed in the shock is much greater than the sound speed. The neutral density in the shock front is $`4n_{\mathrm{n0}}`$, and using $`m_\mathrm{n}=2.33`$ amu for the assumed H2 gas with a 10% abundance of He by number, we therefore have $`L_{\mathrm{drag}}(\mathrm{g})0.3\times 10^{11}`$cm. Hence, the lag of the large grains behind the neutrals on these length scales is understood for this model.

The state of this model at $`t=2.63\times 10^5`$s is shown in Figures $`11a`$ and $`11b`$, and also at $`t=5.5\times 10^5`$s, in Figures $`11c`$ and $`11d`$. By these times the adjustment by plasma drift is approaching an asymptotic state, with an ion and small-charged grain mediated magnetic precursor still ahead of the shock. A particularly interesting feature in the magnetic field profile can be seen in Figures $`11b`$ and $`11d`$, which show that, because of the diffusion of the plasma with respect to the neutrals, magnetic field has “leaked forward” ahead of the shock. This decreases the magnetic field strength $`B_x`$ immediately behind the shock front relative to the far downstream (frozen-in) value $`4B_0`$.

The lag of the large dust grains (crosses) behind the neutrals (circles) is again apparent in Figures $`11b`$ and $`11d`$. The distance over which the motion of the large grains is delayed with respect to the neutrals, because of large grain inertia, is consistent with, and explained by, the value for $`L_{\mathrm{drag}}(\mathrm{g})`$.

Finally, we point out that the results presented in Figures $`9a`$-$`9c`$ are entirely consistent with the model results shown in Figures $`10a`$-$`10d`$ and $`11a`$-$`11d`$. This is so because the small-scale structure has a maximum extent that is basically the size of the magnetic precursor, which has a width $`10^{12}`$ cm. However, examination of Figures 10 and 11 reveal that on scales greater than this, the values of the various physical quantities are the same as the upstream and downstream values for the J shock. Hence, on the large scale the smaller substructure is unresolved, and the model appears as a discontinuous MHD shock, such as that presented in Figure 9.

This model demonstrates our ability to follow the formation of multifluid MHD shocks, including short-lived phenomena. Detailed studies, for a wide variety of parameters and physical conditions, will be presented in later work.

## 4 Summary

In this paper — the first of a series — we have discussed the formulation of a time-dependent, numerical MHD code designed to model the formation and evolution of perpendicular MHD shocks in dense molecular gas. Chief among the various features of our code is that it is a multifluid code, incorporating the dynamics of neutral, ion, and six different types of grain fluids. The grain fluids consist of populations with two distinct radii, each with three different charge states. In addition to including for the effects of elastic collisions with gas particles, we also account for the transfer of mass and momentum between the grain species due to capture of charged particles on grain surfaces. We have presented the equations required to model this physical system, and also described simplifying approximations adopted within the numerical code.

To establish the accuracy of our numerical code, we have presented several test models and compared their results to exact, analytical solutions. In general, we found good agreement between the numerical and analytical solutions. This included a test of the formation of a C shock in a dust-free, weakly ionized gas, whose analytical solution was determined by Chernoff (1987).

We also presented an illustrative model of a multifluid flow in which inertia of the large grains was shown to be important over a finite scale of the flow. The possible effects that differential motion of fluids in shocks may have on the state of the neutral gas and other species will be investigated in more detail elsewhere. The effect of grains (including their inertia) on the speed and propagation of signals in dense gas will also be a topic of discussion in future work. Using this numerical code as a tool to aid in further research, we intend in subsequent papers to study the relevant physics and chemistry of the formation and evolution of various types of nonlinear MHD waves and shocks in molecular clouds, protostellar cores, and disks.

We gratefully acknowledge support from the New York Center for Studies on the Origins of Life (NSCORT), and the Dept. of Physics, Applied Physics, and Astronomy at RPI, under NASA grant NAG5-7598. Thoughtful comments by T. W. Hartquist are also gratefully acknowledged.

## Appendix A Chemical Network


format: text

text:
# Boundary of two mixed Bose-Einstein condensates

## I Introduction

The experimental realization of Bose Einstein condensation in trapped dilute gases  has allowed to investigate variety of properties of quantum fluids both theoretically and experimentally. In the recent years, it has become possible to produce and explore mixtures of Bose-Einstein condensates corresponding to different internal states .

Theoretical treatment of mixtures  assures that depending on the relative strength of interactions inside each condensate and between them it is possible to observe spatial separation. Experimental realization of such systems  has given an opportunity to study both equilibrium properties and dynamics of separation. Although, structure of the boundary has been qualitatively analyzed in  an analytical distribution of densities has not been derived yet. Authors of  discussed asymptotic behavior of densities far from the boundary, and on the basis of these estimations analyzed surface tension.

In the present paper, we explore the boundary between two repulsively interacting condensates at zero temperature in two limit cases corresponding weak and strong separation of the condensates. We show that in the case of weak separation it is possible to derive equations for densities, and using iteration method solve them analytically. We show that asymptotic behavior of our solution coincides with those predicted in . For the sake of completness, we also provide the solutions of these equations in the case of strong separation. Importance of the solution lies in the possibility to explore quantitatively different types of excitations on the boundary. The structure of the boundary in both cases allows to consider one-particle excitations as well as surface waves associated with the boundary. Obtained expressions for dispersion relation of surface waves can be used to explore that phenomenon experimentally.

The Hamiltonian describing the mixture of two weakly interacting Bose-Einstein condensates can be written in the form

$$\begin{array}{c}\hfill H=\underset{i=1,2}{}𝑑𝐫\mathrm{\Psi }_i^+\left[\frac{\mathrm{}^2^2}{2m_i}+V_i(𝐫)+\frac{u_i}{2}\mathrm{\Psi }_i^+\mathrm{\Psi }_i\right]\mathrm{\Psi }_i\\ \hfill +u_{12}𝑑𝐫\mathrm{\Psi }_1^+\mathrm{\Psi }_2^+\mathrm{\Psi }_1\mathrm{\Psi }_2\end{array}$$
(1)
Here, $`u_i=4\pi \mathrm{}^2a_i/m_i>0`$ characterizes the interaction inside each condensate, $`u_{12}=2\pi \mathrm{}^2a_{12}(m_1+m_2)/(m_1m_2)>0`$ – the intercondensate interaction, $`m_i`$ – mass of a particle of each condensate, $`a_i`$, $`a_{12}`$ – corresponding scattering lengths, $`V_i(𝐫)`$ – external trapping potentials. Theoretical treatment of the mixtures  has shown that separation takes place when $`u_{12}/\sqrt{u_1u_2}>1`$.

Starting with the Hamiltonian (1) we get Gross-Pitaevskii equations for condensate wave functions:

$$\begin{array}{c}\hfill i\mathrm{}\frac{\mathrm{\Psi }_1}{t}=\left(\frac{\mathrm{}^2^2}{2m_1}+V_1(𝐫)+u_1|\mathrm{\Psi }_1|^2+u_{12}|\mathrm{\Psi }_2|^2\right)\mathrm{\Psi }_1,\\ \hfill i\mathrm{}\frac{\mathrm{\Psi }_2}{t}=\left(\frac{\mathrm{}^2^2}{2m_2}+V_2(𝐫)+u_2|\mathrm{\Psi }_2|^2+u_{12}|\mathrm{\Psi }_1|^2\right)\mathrm{\Psi }_2\end{array}$$
(2)
As we are interested in studying of stationary solutions of these equations, then assuming as usual $`\mathrm{\Psi }_j\mathrm{exp}(i\mu _jt)`$, where $`\mu _j`$ – chemical potentials of the condensates, we obtain two coupled nonlinear equations for densities of gases $`n_i(𝐫)=|\mathrm{\Psi }_i(𝐫)|^2`$:

$$\begin{array}{c}\mu _1=\frac{\mathrm{}^2}{2m_1}\frac{^2\sqrt{n_1}}{\sqrt{n_1}}+V_1(𝐫)+u_1n_1+u_{12}n_2,\hfill \\ \mu _2=\frac{\mathrm{}^2}{2m_2}\frac{^2\sqrt{n_2}}{\sqrt{n_2}}+V_2(𝐫)+u_2n_2+u_{12}n_1\hfill \end{array}$$
(3)
These equations are essentially nonlinear, so to find solutions we need to make some simplifications. We assume that $`V_1(𝐫)=V_2(𝐫)`$ and consider the case when the size of a boundary between condensates much less than characteristic length of the trap. In the case of a parabolic trap potential and Thomas-Fermi regime, it means that $`dR_{TF}`$, where $`d`$ is the size of the boundary, and $`R_{TF}`$ is the Thomas-Fermi radius of the atomic cloud. Physically, it helps to avoid the effect of the potential on the form of the boundary. To simplify calculations further, we also suppose that separation takes place in one dimension, $`z`$. This allows us to write the system of equations (3) in the following form:

$$\begin{array}{c}\mu _1=\frac{\mathrm{}^2}{2m_1\sqrt{n_1}}\frac{d^2}{dz^2}\sqrt{n_1}+u_1n_1+u_{12}n_2,\hfill \\ \mu _2=\frac{\mathrm{}^2}{2m_2\sqrt{n_2}}\frac{d^2}{dz^2}\sqrt{n_2}+u_2n_2+u_{12}n_1\hfill \end{array}$$
(4)
Although, there is no explicit trapping potential in these equations, we have to impose external conditions on the solutions, that implicitly take it into account. Let us assume that separation takes place along the $`z`$ direction, and the condensate with the label “$`1`$” is to the right, and the one with the label “$`2`$” is to the left of the boundary. Then, asymptotically we require:

$$\begin{array}{c}n_1(z+\mathrm{})n_{10},n_1(z\mathrm{})0,\hfill \\ n_2(z\mathrm{})n_{20},n_2(z+\mathrm{})0,\hfill \end{array}$$
(5)
where $`n_{10}`$, $`n_{20}`$ – equilibrium densities of condensates far from the boundary. Substituting these conditions in equations (4) for densities, and using the obvious condition of equilibrium we obtain

$$\begin{array}{c}\mu _1=u_1n_{10},\mu _2=u_2n_{20},\hfill \\ P_1=u_1n_{10}^2/2=P_2=u_2n_{20}^2/2\hfill \end{array}$$
(6)
Here, we used well-known expression for the pressure of a homogeneous weakly-interacting Bose gas. That condition connects equilibrium densities of condensates far from the boundary.

To reduce the number of parameters in equations (4) we notice that it is possible to exclude difference of masses by the change:

$$\begin{array}{c}u_1^{}=u_1m_1/m_2,u_2^{}=u_2m_2/m_1,\hfill \\ n_1^{}=n_1\sqrt{m_2/m_1},n_2^{}=n_2\sqrt{m_1/m_2},\hfill \\ \mu _1^{}=u_1^{}n_{10}^{},\mu _2^{}=u_2^{}n_{20}^{},\hfill \\ m^{}=\sqrt{m_1m_2}\hfill \end{array}$$
(7)
We get the same equations (4) and conditions (5), (6) but for quantities with stars and with the same mass $`m^{}`$. To simplify notations, in the following discussion we omit stars. The generalization for different masses can be easily done by following the above rules.

We can solve equations for densities (4) with conditions (5),(6) analytically in two limit cases for weak separation when $`\mathrm{\Delta }=u_{12}/\sqrt{u_1u_2}11`$, and strong separation when $`\mathrm{\Delta }1`$. Notice, that both parameters $`u_{12}`$ and $`\mathrm{\Delta }`$ are not affected by the above procedure of mass difference excluding.

## II Weak separation

Let us consider weak separation when condition $`\mathrm{\Delta }=u_{12}/\sqrt{u_1u_2}11`$ is satisfied. We expect that in the simplest case when $`u_1=u_2`$ the whole density of a gas $`n(z)=n_1(z)+n_2(z)`$ is approximately constant. That is why, to find the solution it is natural instead of $`n_1`$, $`n_2`$ to introduce other quantities and solve equations (4) using small parameter $`\mathrm{\Delta }`$. Consider the functions:

$$\begin{array}{c}\rho =(u_1/u_2)^{1/4}n_1+(u_2/u_1)^{1/4}n_2,\hfill \\ g=\left[(u_1/u_2)^{1/4}n_1(u_2/u_1)^{1/4}n_2\right]/\rho \hfill \end{array}$$
(8)
Conditions (5), (6) give us simple asymptotic behavior of these functions:

$$\begin{array}{c}\rho (z\pm \mathrm{})\rho _0=\sqrt{n_{10}n_{20}},\hfill \\ g(z\pm \mathrm{})\pm 1\hfill \end{array}$$
(9)
Densities of condensates “1” and “2” are easily obtained if functions $`\rho `$ and $`g`$ are known:

$$\begin{array}{c}n_1=(u_2/u_1)^{1/4}\rho [1+g]/2,\hfill \\ n_2=(u_1/u_2)^{1/4}\rho [1g]/2\hfill \end{array}$$
(10)
It is straightforward to derive equations for $`\rho `$ and $`g`$:

$$\begin{array}{c}\hfill \frac{\sqrt{\rho }^{\prime \prime }}{\sqrt{\rho }}\frac{\rho }{\rho _0}\left[1+\alpha g+\frac{\mathrm{\Delta }}{2}(1g^2)\right]\\ \hfill =\frac{g^{\mathrm{\hspace{0.17em}2}}}{4(1g^2)}1\alpha g,\\ \hfill \frac{g^{\prime \prime }}{1g^2}+\frac{2\sqrt{\rho }^{}g^{}}{\sqrt{\rho }(1g^2)}+\frac{gg^{\mathrm{\hspace{0.17em}2}}}{(1g^2)^2}+2\alpha \\ \hfill =\frac{\rho }{\rho _0}\left[2\alpha +\mathrm{\Delta }(\alpha g)\right]\end{array}$$
(11)
Here, $`f^{}=\xi _0\frac{df}{dz}`$, and we also introduced $`\alpha =\frac{\sqrt{u_1}\sqrt{u_2}}{\sqrt{u_1}+\sqrt{u_2}}`$ and $`1/\xi _0^2=m(u_1u_2)^{1/4}(\sqrt{u_1}+\sqrt{u_2})\rho _0/\mathrm{}^2`$.

To find asymptotic solutions of equations (11) in the case of $`\mathrm{\Delta }1`$ we use iteration method. Let us suppose that terms with derivatives of $`\rho `$ are much smaller than the others. We will justify this assumption in the end of calculations.

Neglecting terms with derivatives of $`\rho `$, from equations (8) we get

$$\begin{array}{c}\rho /\rho _0=1\frac{g^{\mathrm{\hspace{0.17em}2}}+2\mathrm{\Delta }(1g^2)^2}{4(1g^2)(1+\alpha g)},\hfill \\ \frac{g^{\prime \prime }}{1g^2}+\frac{2g+\alpha (1+g^2)}{2(1g^2)(1+\alpha g)}g^{\mathrm{\hspace{0.17em}2}}+\frac{\mathrm{\Delta }(1\alpha ^2)g}{1+\alpha g}=0\hfill \end{array}$$
(12)
The equation for $`g`$ can be solved by substitution $`g^{}=f`$, so taking into account conditions (9) we get that $`g`$ is the solution of the equation

$$g^{}=\sqrt{\mathrm{\Delta }(1\alpha ^2)}\frac{(1g^2)}{\sqrt{1+\alpha g}}$$
(13)
At this point it is clear that assumptions we made were correct. Namely, we obtain that $`\sqrt{\rho }^{}\mathrm{\Delta }^{3/2}`$, $`\sqrt{\rho }^{\prime \prime }\mathrm{\Delta }^2`$, $`g^{}\sqrt{\mathrm{\Delta }}`$, so neglected terms by $`\mathrm{\Delta }`$ times less than the others.

Now we can write down solutions for both $`\rho `$ and $`g`$ in the parametric form

$$\begin{array}{c}\hfill 1\frac{\rho }{\rho _0}=\frac{\mathrm{\Delta }}{4}\frac{1g^2}{(1+\alpha g)^2}[3\alpha ^2+2\alpha g],\\ \hfill \frac{zz_0}{\xi _0/\sqrt{\mathrm{\Delta }(1\alpha ^2)}}=\frac{\sqrt{1+\alpha }}{2}\mathrm{ln}|\frac{\sqrt{1+\alpha g}+\sqrt{1+\alpha }}{\sqrt{1+\alpha g}\sqrt{1+\alpha }}|\\ \hfill \frac{\sqrt{1\alpha }}{2}\mathrm{ln}|\frac{\sqrt{1+\alpha g}+\sqrt{1\alpha }}{\sqrt{1+\alpha g}\sqrt{1\alpha }}|\end{array}$$
(14)
Here, the second equation is the solution of (13). In these expressions an arbitrary constant $`z_0`$ defines the position of the boundary, and in the case of finite geometry with given average number of particles in condensates can be obtained from the condition of equal pressures (6). There is also physically obvious symmetry $`zz_0z_0z`$, $`gg`$, $`\alpha \alpha `$.

Finally, the densities of condensates have the following form:

$$\begin{array}{c}n_1=\frac{n_{10}}{2}\left(1\frac{\mathrm{\Delta }}{4}\frac{1g^2}{(1+\alpha g)^2}[3\alpha ^2+2\alpha g]\right)\left[1+g\right],\hfill \\ n_2=\frac{n_{20}}{2}\left(1\frac{\mathrm{\Delta }}{4}\frac{1g^2}{(1+\alpha g)^2}[3\alpha ^2+2\alpha g]\right)\left[1g\right]\hfill \end{array}$$
(15)
With the parametric equation for $`g`$, these densities are the main result of the paper. On fig With the use of the desribed method, it is possible to derive expressions for densities in the next orders of the small parameter $`\mathrm{\Delta }`$ in the form of asymptotic series. As follows from the above estimations of neglected terms in equations (11) the next order is proportional to $`\mathrm{\Delta }^2`$. The typical dependence of densities on the distance from the boundary is shown in Fig. 1 and Fig. 2.

It is interesting to notice that the whole density $`n=n_1+n_2`$ at some values of parameters has a well on the boundary as shown in Fig. 1. In the case of weak separation we see that it is broad and shallow. This result is a consequence of the interparticle interactions. The interparticle potential for each condensate acts in some sense as a wall, so we expect that probability of a particle to be close to the boundary decreases, which means lower density. The existence of such a well allows to consider the possibility of confining of a particle of another sort in the vicinity of the boundary, which is discussed below.

Although, we can not further simplify solutions (14), it is interesting to derive asymptotic behavior of these functions in particular limits:

$$\begin{array}{c}1g(z+\mathrm{})\mathrm{exp}[\frac{2z\sqrt{\mathrm{\Delta }}}{\xi _2}],\hfill \\ 1+g(z\mathrm{})\mathrm{exp}[\frac{2z\sqrt{\mathrm{\Delta }}}{\xi _1}],\hfill \\ g(zz_0)\frac{(zz_0)\sqrt{\mathrm{\Delta }(1\alpha ^2)}}{\xi _0}\hfill \end{array}$$
(16)
Here, $`\xi _i=\mathrm{}/\sqrt{2m_i\mu _i}`$ – correlation lengths of condensates, defined by chemical potential $`\mu _1`$ and $`\mu _2`$, and we used (7) to include mass difference. The asymptotic behavior far from the boundary of each condensate is defined by its correlation length. As we see, the size of the boundary can be approximated as $`d(\xi _1+\xi _2)/\sqrt{\mathrm{\Delta }}`$ and for $`\mathrm{\Delta }1`$ appears to be much larger than correlation lengths.

Solutions have the simplest form when $`\alpha =0`$:

$$\begin{array}{c}n_1(z)=\frac{n_0}{2}\left(1\frac{3\mathrm{\Delta }}{4\mathrm{cosh}^2[\frac{z\sqrt{\mathrm{\Delta }}}{\xi _0}]}\right)\left(1+\mathrm{tanh}[\frac{z\sqrt{\mathrm{\Delta }}}{\xi _0}]\right),\hfill \\ n_2(z)=\frac{n_0}{2}\left(1\frac{3\mathrm{\Delta }}{4\mathrm{cosh}^2[\frac{z\sqrt{\mathrm{\Delta }}}{\xi _0}]}\right)\left(1\mathrm{tanh}[\frac{z\sqrt{\mathrm{\Delta }}}{\xi _0}]\right)\hfill \end{array}$$
(17)
Here, we choose $`z_0=0`$.

## III Strong separation

To analyze the case of strong separation $`\mathrm{\Delta }=u_{12}/\sqrt{u_1u_2}1`$ (we use the same notation but for another quantity) we start from density equations (4). In this case we expect that density on the boundary will be approximately zero because interparticle interactions make it almost impossible for one condensate to penetrate inside the other. To estimate the density of condensates on the boundary we can use the fact that second derivatives of wave functions should be approximately zero there. That makes system (4) a set of two linear equations with solution:

$$\begin{array}{c}n_{1B}=\frac{n_{10}}{\mathrm{\Delta }+1}\frac{n_{10}}{\mathrm{\Delta }}n_{10},\hfill \\ n_{2B}=\frac{n_{20}}{\mathrm{\Delta }+1}\frac{n_{20}}{\mathrm{\Delta }}n_{20}\hfill \end{array}$$
(18)
This allows us in the zero approximation to use simple conditions for the densities $`n_1(z0)=n_2(z0)=0`$. Then equations (4) have simple form:

$$\begin{array}{c}\mu _1=\frac{\mathrm{}^2}{2m_1\sqrt{n_1}}\frac{d^2}{dz^2}\sqrt{n_1}+u_1n_1,\text{ for }z0,\hfill \\ \mu _2=\frac{\mathrm{}^2}{2m_2\sqrt{n_2}}\frac{d^2}{dz^2}\sqrt{n_2}+u_2n_2,\text{ for }z0\hfill \end{array}$$
(19)
Solutions are easily obtained:

$$\begin{array}{c}n_1(z0)=n_{10}\mathrm{tanh}^2[\frac{z}{\sqrt{2}\xi _1}],\hfill \\ n_2(z0)=n_{20}\mathrm{tanh}^2[\frac{z}{\sqrt{2}\xi _2}]\hfill \end{array}$$
(20)
Here, $`n_{10}`$ and $`n_{20}`$ are connected by the condition of equal pressures (6), and we choose the position of the boundary at $`z=0`$.

The size of the boundary in this case is approximately $`d2\sqrt{2}(\xi _1+\xi _2)`$. The dependence of the densities on the distance from the boundary is shown in Fig. 3. As in the previous section, we see that there is again a well in the whole density but in the case of strong separation it becomes narrower and deeper in comparison with the one for weak separation.

## IV One-particle excitations on the boundary

The existence of a well in the whole density allows to consider confining of a particle of another sort in the vicinity of the boundary. As a general property of a quantum-mechanical motion in a one-dimensional well there always exists such a confined state. As an example, we consider the simplest case when $`\alpha =0`$ and a particle of another sort interacts with both condensates repulsively with the same constant $`\lambda `$. The Schrödinger equation for the wave function of a particle with the mass $`M`$ has the form:

$$\left[\frac{\mathrm{}^2}{2M}\frac{d^2}{dz^2}+\lambda n(z)\right]\varphi =E\varphi $$
(21)
We can solve equation (21) for weak and strong separation cases simultaneously. It has universal form

$$\frac{d^2\varphi }{dz^2}+\frac{2M}{\mathrm{}^2}\left(ϵ+\frac{U_0}{\mathrm{cosh}^2[\beta z]}\right)\varphi =0,$$
(22)
where $`ϵ=E\lambda n_0`$, and $`U_0=3\mathrm{\Delta }\lambda n_0/4`$, $`\beta =\sqrt{\mathrm{\Delta }}/\xi _0`$ for weak; and $`U_0=\lambda n_0`$, $`\beta =1/(\sqrt{2}\xi _0)`$ for strong separation. The spectrum of energy $`ϵ`$ is well-known:

$$\begin{array}{c}ϵ_j=\mu \frac{\mathrm{\Delta }m}{4M}\left[(1+2j)+\sqrt{1+3\frac{M\lambda }{mu}}\right]^2\text{“weak”},\hfill \\ ϵ_j=\mu \frac{m}{8M}\left[(1+2j)+\sqrt{1+8\frac{M\lambda }{mu}}\right]^2\text{“strong”}\hfill \end{array}$$
(23)
where $`\mu =un_0`$ – chemical potential, $`j=0,1,\mathrm{}`$ and the condition that an expression in $`[\mathrm{}]`$ is positive defines the upper limit for $`j`$. There is always at least one state with $`j=0`$.

## V Surface waves on the boundary

There is also another type of excitations associated with the boundary. As we see condensate density distributions give rise to nonzero surface tension which was previously analyzed in , and we use the expression for surface tension derived there:

$$\sigma =\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑z\underset{i=1,2}{}\frac{\mathrm{}^2}{2m_i}\left(\frac{d\sqrt{n_i}}{dz}\right)^2$$
(24)
Substituting expressions (15) for densities in the case of weak separation and taking into account only first nonzero order in $`\mathrm{\Delta }`$ we obtain for surface tension

$$\sigma _w=\frac{P\sqrt{\mathrm{\Delta }}}{6}\left[\begin{array}{c}\hfill (\xi _1+\xi _2)\frac{2\alpha ^21+\sqrt{1\alpha ^2}}{\alpha ^2}\\ \hfill (\xi _1\xi _2)\frac{1\sqrt{1\alpha ^2}}{\alpha }\end{array}\right],$$
(25)
where $`\alpha =\left(m_1\sqrt{u_1}m_2\sqrt{u_2}\right)/\left(m_1\sqrt{u_1}+m_2\sqrt{u_2}\right)`$ and we used (7) to include mass difference; $`\xi _1`$, $`\xi _2`$ – correlation lengths of condensates; $`P`$ – the pressure given by (6). When $`\alpha =0`$, $`\sigma _w=P\sqrt{\mathrm{\Delta }}(\xi _1+\xi _2)/4`$

In the case of strong separation we use expressions (20) to get surface tension

$$\sigma _s=\frac{P\sqrt{2}}{3}(\xi _1+\xi _2)$$
(26)
Let us notice that expressions (25) and (26) differ from the ones obtained in . Although, in qualitative sense our expressions coincide with those of , using our method of solution we can get a general expression applicable for variety of parameters, and for example retrieve the correct numerical factor for the case considered in . As follows from the general expression for surface tension (24), we need to know the behavior of densities not only far but also in the vicinity of the boundary. That is why, estimate character of expressions for densities in  could only give qualitative answer for surface tension.

For velocities smaller than the speed of sound we can consider gas as incompressible, so it is possible to write down the usual hydrodynamic equations and find dispersion relation of surface waves on the boundary associated with surface tension. Suppose that we have condensates in a “box” and the position of the boundary is defined by the distances $`L_1`$ and $`L_2`$ from the “box” walls, where labels “1”, “2” correspond to the side with the same condensate, and the box sizes for other directions are $`L_x`$, $`L_y`$. Taking into account the fact that velocities are zero on the walls we get the dispersion relation

$$\omega ^2(k)=\frac{\sigma k^3}{\rho _1\mathrm{coth}(kL_1)+\rho _2\mathrm{coth}(kL_2)},$$
(27)
where $`\sigma `$ is the surface tension; $`\rho _{10}=m_1n_{10}`$, $`\rho _{20}=m_2n_{20}`$ – mass densities; wave vector along the boundary $`k=\sqrt{(\pi n_x/L_x)^2+(\pi n_y/L_y)^2}`$, $`n_x,n_y=0,1\mathrm{}`$, not $`n_x=n_y=0`$; $`L_1`$, $`L_2`$ – sizes of condensates in the direction perpendicular to the boundary. We assume that $`L_1`$, $`L_2`$ are much larger than the size of the boundary $`d`$.

For dispersion relation we consider two limit cases corresponding to wavelengths $`\lambda L_i`$ ($`kL_i1`$) and $`\lambda L_i`$ ($`kL_i1`$). Let us notice that the second case is possible only if $`L_iL_{x,y}`$.

In the first case we obtain

$$\omega (k)=\left(\frac{\sigma }{\rho _{10}+\rho _{20}}\right)^{1/2}k^{3/2}$$
(28)
In the long wavelength limit we get

$$\omega (k)=\left(\frac{\sigma L_1L_2}{\rho _{10}L_2+\rho _{20}L_1}\right)^{1/2}k^2$$
(29)
As we see, for weak separation $`\omega \mathrm{\Delta }^{1/4}`$, and surface waves are relatively “soft” in that case, which enables to consider them as a dissipative channel for other condensate excitations.

## VI Conclusion

We performed the analysis of the boundary of two Bose-Einstein condensates interacting repulsively in the limit cases corresponding weak and strong separation at zero temperature.

For weak separation we obtained solutions (15) of two coupled nonlinear Gross-Pitaevskii equations using small parameter $`\mathrm{\Delta }`$. The solutions show that the penetration depth of the condensate “$`i`$” inside the other is estimated as $`\xi _i/\sqrt{\mathrm{\Delta }}`$, so that the size of the boundary $`d(\xi _1+\xi _2)/\sqrt{\mathrm{\Delta }}`$ is much larger than correlation lengths, which was obtained experimentally . There is also a well in the full density profile, which is the consequence of the wave function behavior near the boundary. In general, the proposed method of obtaining density distributions for the weak separation case if desired can be expanded to obtain expressions for the next orders of $`\mathrm{\Delta }`$.

We also considered the case of strong condensate separation but restricted ourselves to the zero order approximation. In that case the size of the boundary is $`d(\xi _1+\xi _2)`$, and the whole density of gases goes approximately to zero on the boundary.

The existence of the well in the density profile at some parameters allowed us to consider one-particle excitations on the boundary. Using the expressions for densities we found excitation spectrum of a particle in the simplest case when constants of interaction of the particle with condensates are the same, and distribution of densities on the boundary corresponds to $`\alpha =0`$. The generalization to other cases can be easily done with the use of expressions (15),(20) for density distributions. Let us notice that the existence of a potential well depends on the interaction constants of a particle with condensates as well as the relation between interaction constants of the condensates. Observation of such confined states can be possible only if temperature are smaller than the potential well.

It was shown that there also exist collective excitations associated with surface tension. The expressions for surface tension were obtained for both weak and strong separation. The dispersion relation for surface waves was analysed in the case where condensates fill finite volumes. The dispersion relation has different forms in two cases corresponding to short- and long-wavelength limit. In the case of weak separation “soft” surface modes can be a dissipative channel of other condensate excitations.

Author gratefully acknowledges valuable discussions with Yu.M. Kagan and other members of the Theoretical Division of RRC “Kurchatov Institute”, Russia.