text:
# Magneto-Centrifugal Launching of Jets from Accretion Disks. II. Inner Disk-Driven Winds

## 1. Introduction

### 1.1. The Quest for Large-Scale Wind Structure

In a previous paper (Krasnopolsky, Li, & Blandford 1999; hereafter Paper I), we described a modification of the MHD code and its application to the problem of launching a magneto-centrifugal wind from a Keplerian disk. The disk was idealized as a boundary condition at the equatorial plane, on which open magnetic fields are firmly anchored at all radii inside the simulation box. Mass was injected supersonically but sub-Alfvénically onto the open field lines at a prescribed rate and accelerated magneto-centrifugally along the field lines to produce a high speed wind. Our approach followed that of Ustyugova et al. (1995) and Ouyed & Pudritz (1997). It permitted a decoupling of the dynamics of the wind from that of the disk and allowed for a thorough investigation of the structure of the wind. Previous numerical studies of magneto-centrifugal winds were limited to a relatively small region not far from the launching surface. In this paper, we seek to extend the wind solutions to large, observable distances, and to relate the observable flow quantities on large scales to the launching conditions.

To study the large-scale structure of a disk-driven magneto-centrifugal wind numerically, a technical problem must first be overcome. This problem is illustrated in Fig. 1 (taken from Paper I), where a disk-wind solution representative of those found in existing time-dependent simulations (e.g., Romanova et al. 1997; Ouyed & Pudritz 1997; Ustyugova et al. 1999) is shown. Note that a large fraction of the wind remains sub–fast-magnetosonic in the computational realm. Information from the edge of the simulation box can propagate upstream in a sub–fast-magnetosonic region, creating an undesirable coupling between the box boundaries and the wind-launching. This coupling, we believe, is at the root of the sensitive dependence of wind solutions on the size and shape of the simulation box found previously by us and others (see §2.1 and Ustyugova et al. 1999 for a discussion). In some cases, merely changing the shape of the box can destabilize a wind and induce chaotic flow fragmentation. This dependence must be eliminated before one can study the large scale structure of the magneto-centrifugal wind with confidence.

One way to overcome the above difficulty is to restrict the wind launching to the inner region of an accretion disk. In the specific context of star formation, there are reasons to believe that the bulk of the outflow is indeed driven from an inner disk. There are two popular magneto-centrifugal models for the production of optical jets, such as the spectacular HH 30 jet studied in detail with the HST (Burrows et al. 1996; Ray et al. 1996). The models differ on where the wind-driving open magnetic fields are anchored. If the open field lines are anchored in a narrow region of disk near the corotation radius of the stellar magnetosphere, then the wind is called an “X-wind” (Shu et al., 2000). The X-wind serves the fundamental purpose of removing angular momentum from the central star and keeping it rotating at a rate well below the breakup rate, as observed. If, on the other hand, the field lines are anchored over a wider region of the disk, then a “disk-wind” is produced (eg., Königl & Pudritz 2000). The disk-wind may be primarily responsible for driving the mass accretion through the portion of the disk where the wind is launched. Wind launching should be easier from the inner part of protostellar disks, where the temperature is high enough ($`T>10^3\text{K}`$) that thermal ionization of alkali metals ensures that the disk material and wind-launching magnetic field are well coupled. In the cooler outer part, the bulk of disk material is too weakly ionized to participate in the wind launching, although in the surface layer and at large radii (say $`10^2\text{AU}`$; Wardle & Königl 1993) protostellar X-rays and Galactic cosmic rays (if not excluded by the wind) may provide a level of ionization sufficient for good coupling. We will not consider in this paper the possible outer disk wind, which is probably less relevant for the fast jets and winds of hundreds of $`\text{km}\text{s}^1`$ that we are interested in. Such a wind can modeled with our simulation setup (see §2.1) if desired since the size of the launching region is fully adjustable.

The large-scale structure of magneto-centrifugal winds is a prerequisite for quantitatively modeling a variety of outflow phenomena observed around YSOs, including optical jets, neutral atomic winds, and bipolar molecular outflows. It is the observational data on these phenonema that will constrain the conditions at the wind launching surface, which are poorly known. Particularly tight constraints are expected from high spectral as well as high spatial resolution observations of forbidden line emission from T Tauri winds. A recent example of such observations is that of Bacciotti et al. (2000). They constructed synthetic two-dimensional images of the DG Tau wind based on HST/STIS data in several forbidden lines (including \[OI\]$`\lambda \lambda `$6363,6300, \[SII\]$`\lambda `$6731, and \[NII\]$`\lambda `$6583) and Hα in several radial velocity intervals, with a spatial resolution of order $`0.1\mathrm{}`$ (or about $`15\text{AU}`$ at the distance of DG Tau). To take full advantage of such observations, one needs a reliable model that can follow the flow from the (small) launching surface to observable distances. As a first step towards such a model, we restrict our attention to those magneto-centrifugal winds that can reach a steady state, even though most observed jets are time-dependent, possibly episodic (Reipurth & Raga, 1999). We postpone an investigation of the intrinsically time-dependent winds to a future paper.

Although our intended application is primarily to YSOs, we note that magneto-centrifugal winds driven from a limited region of accretion disks may occur in other astrophysical systems, e.g., cataclysmic variables (CVs) and possibly planetary nebulae. A related jet launching mechanism, through torsional Alfvén wave, was first discussed in Uchida & Shibata (1985), and refined in many subsequent works (see, e.g., Kato, Kudo, & Shibata 2002 and references therein).

### 1.2. Uncertain Jet Formation Efficiency

Asymptotic structure of steady magneto-centrifugal winds has been analyzed by Heyvaerts & Norman (1989) and Shu et al. (1995), among others. The focus of such analyses was on jet formation, through the hoop-stress associated with toroidal magnetic field. While there is a general agreement that a dense, narrow jet can indeed form in the axial region of a magneto-centrifugal wind (see, however, Okamoto  for a different view), it is not clear how the jet properties, such as their density, mass and momentum fluxes, depend on the launching conditions. Indeed, the term “jet” was used by different authors to mean different things. Bogovalov & Tsinganos (1999), for example, refer to the axial region with a more or less uniform density and magnetic field strength and vertical streamlines as the jet, whereas Shu et al. (1995) suggested that the observed jet is simply the dense part of an underlying anisotropic wind above, roughly speaking, some fiducial density that shows up more prominently than the more tenuous equatorial part in the emission of forbidden lines such as \[OI\] and \[SII\] (Shang, Shu, & Glassgold, 1998). We adopt the latter meaning of “jet”, which is more observationally oriented. For simplicity, the fiducial density at the jet boundary is chosen to be $`n_\text{H}=10^4\text{cm}^3`$, comparable to the values derived by Bacciotti, Eisloffel, & Ray (1999) for the famous HH 30 jet111We note that the density derived by Bacciotti et al. (1999) for the HH 30 jet is not constant along the jet, but decreases from a value $`10^5\text{cm}^3`$ close to the exciting source to a value $`10^4\text{cm}^3`$ further out.. An important question is: what fraction of the mass flux of the wind is collimated into the fiducial jet? We will address this question using numerically determined, large-scale wind solutions.

The fraction of mass flux residing in the jet measures, in some sense, the efficiency of jet formation. It can be constrained by observations. Take the HH 30 system as an example. For this object, one can eliminate a jet formation efficiency as low as $`10\%`$ from the following considerations. Bacciotti et al. (1999) showed that the mass flux in the main jet is $`1.7\times 10^9\text{M}_{\mathrm{}}\text{yr}^1`$ ($`1.8\times 10^9\text{M}_{\mathrm{}}\text{yr}^1`$ in the counter jet) close to the exciting source. If only $`10\%`$ of the mass flux resides in the jet, the total mass flux of the wind from both sides would be $`3.5\times 10^8\text{M}_{\mathrm{}}\text{yr}^1`$. If a third (Shu et al., 1995) or less (e.g., Ferreira & Pelletier 1995) of the material accreted through the disk gets ejected in the wind, then the disk mass accretion rate must be $`10^7\text{M}_{\mathrm{}}\text{yr}^1`$ or higher. Such a large accretion rate would be difficult to accommodate by the very weak IR excess from ISO observations (Stapelfeldt & Moneti, 1998), which points to a relatively inactive disk, and by the low disk mass ($`10^3\text{M}_{\mathrm{}}`$) estimated from scattered light and dust emission (Stapelfeldt & Padgett, 2001), although this estimate is uncertain. Increasing the jet mass flux fraction to $`50\%`$ would lower the disk accretion rate to $`2\times 10^8\text{M}_{\mathrm{}}\text{yr}^1`$ (a more typical value for T Tauri disks; Gullbring et al. 1998), if a third of the accreted material is ejected. If a smaller fraction of the accreted material is ejected, then the required accretion rate would be higher. It therefore appears that, at least for HH 30, the jet formation efficiency must be fairly high, of order $`50\%`$ or more. It is not clear whether such a high efficiency can be achieved naturally in a magneto-centrifugal wind.

The rest of the paper is organized as follows. In §2, we present our formulation of the problem of launching magneto-centrifugal winds from only inner disk regions. It is followed by a discussion of a “standard” wind solution which illustrates several generic features of jet formation and large-scale structure of steady magneto-centrifugal winds in §3. In §4, we examine the effects of the mass flux distribution at the launching surface on the jet/wind properties. We conclude and discuss future work in §5.

## 2. Formulation of the Inner Disk-Driven Wind Problem

### 2.1. Implementation of a Finite Launching Surface

As mentioned earlier, we limit the launching of magneto-centrifugal winds from the inner, relatively hot regions of protostellar disks. Numerically, this restriction enables us to extend the wind solution to large distances without having to worry about the undesirable coupling between the simulation box and the launching region, which occurs in the sub fast-magnetosonic region in the lower-right corner of the solution shown in Fig. 1. The reason for that region of flow to remain sub fast-magnetosonic is simple: the wind coming off the outer part of the disk encounters the edge of the simulation box too soon. In other words, the wind in the region simply does not have enough room to get accelerated to the fast magnetosonic speed. This situation remains true as long as the wind is driven off all of the (equatorial) disk plane (as assumed in all previous time-dependent disk-wind simulations), regardless of the box size. This troublesome feature can be eliminated by restricting the wind launching to only the inner region of a disk (within a radius denoted by $`R_L`$), as sketched in Fig. 2. To fill all available space above (and below) the equatorial plane, we demand that the last field line lie exactly on the equatorial plane222Of course, if the wind is confined externally, either by a flaring disk or an ambient medium, the last field line must be modified accordingly. External confinement of magneto-centrifugal winds is an important topic that we plan to explore numerically in the future. It is expected to play a more dominant role in the outflows of embedded sources than in the winds of optically revealed T Tauri stars. Note that the dense material inside the disk beyond the wind launching region prevents the “last” field lines of both sides of the disk from reconnecting. Also, the interface between the wind and disk material outside $`R_L`$ may become Kelvin-Helmholtz unstable. Such additional complications are ignored in this paper.. Wind plasma sliding along the last (horizontal) field line, anchored at the radius $`R_L`$, will become super fast-magnetosonic in the computational realm, provided that the size of the simulation box is sufficiently large. Once the whole fast magnetosonic surface is completely enclosed inside the simulation box, the size and shape of the box should have minimal effects on the structure of the wind, especially near the launching surface, since information cannot propagate upstream in a super fast-magnetosonic region. In this way, we should be able to study the wind structure up to arbitrarily large distances from the source region, limited only by computer time.

There is, however, an additional consideration. Ustyugova et al. (1999) have shown that artificial dependence of simulation results on computational box shape can appear even in a super–fast-magnetosonic flow, if the Mach cones at the boundaries intersect the simulation region. The opening angle $`\phi `$ of the Mach cone around a streamline is given (for a cold flow) by $`\mathrm{tan}^2\phi =1/(M_f^21)`$, where $`M_f`$ is the fast Mach number. As the flow becomes more super–fast-magnetosonic, the cone gets narrower, more forward-directed, and has less chance of intersecting the simulation region. Our simulations, designed to produce clearly super-fast flows, fulfill this requirement everywhere except for a small region around the axis. As we discuss below, the axial region receives a special treatment, and does not appear to create any numerical problems.

### 2.2. Equations

For completeness and to define symbols, we list below the time-dependent MHD wind equations that we solve numerically:

$`{\displaystyle \frac{\rho }{t}}+\left(\rho 𝐯\right)`$ $`=`$ $`0`$ (1)
$`\rho {\displaystyle \frac{𝐯}{t}}+\rho \left(𝐯\right)𝐯`$ $`=`$ $`p+\rho \mathrm{\Phi }_g+𝐣\times 𝐁/c`$ (2)
$`{\displaystyle \frac{𝐁}{t}}`$ $`=`$ $`\times \left(𝐯\times 𝐁\right)=\times 𝓔`$ (3)
$`{\displaystyle \frac{u}{t}}+\left(u𝐯\right)`$ $`=`$ $`p𝐯`$ (4)
$`p`$ $`=`$ $`(\gamma 1)u`$ (5)
where

$`\rho `$ $`=`$ matter density
$`𝐯`$ $`=`$ velocity flow field
$`𝐁`$ $`=`$ magnetic field
$`𝐣`$ $`=`$ $`(c/4\pi )\times 𝐁=\text{current density}`$
$`𝓔`$ $``$ $`𝐯\times 𝐁=c𝐄\text{where}𝐄\text{is the electric field}`$
$`\mathrm{\Phi }_g`$ $`=`$ gravitational potential
$`p`$ $`=`$ thermal pressure
$`u`$ $`=`$ internal energy density (per unit volume)
$`\gamma `$ $`=`$ adiabatic index
We have adopted a cylindrical coordinate system $`(R,\varphi ,z)`$, with the spherical radius given by $`r=\sqrt{R^2+z^2}`$. The gravity field is smoothed inside a sphere of radius $`r_g`$, through

$$\mathrm{\Phi }_g=\frac{\mathrm{\Omega }_0^2}{2}(r^23r_{g}^{}{}_{}{}^{2}),$$
(6)
where $`\mathrm{\Omega }_0`$ and $`r_g`$ are fixed parameters, whose product $`v_{\text{Ko}}=\mathrm{\Omega }_0r_g`$ gives the scale for the rotational speed $`v_\text{K}\sqrt{R\left(\mathrm{\Phi }_g/R\right)_{z=0}}=\mathrm{\Omega }_0R`$ on the disk inside $`r_g`$. Outside $`r_g`$, the gravity is not smoothed, and the potential has the usual form

$$\mathrm{\Phi }_g=\frac{\mathrm{\Omega }_0^2r_{g}^{}{}_{}{}^{3}}{r}.$$
(7)
It yields a Keplerian speed of

$$v_\text{K}=v_{\text{Ko}}\left(\frac{r_g}{R}\right)^{1/2}$$
(8)
on the disk. The parameter $`\mathrm{\Omega }_0`$ is simply the angular speed at the inner edge of the (unsoftened) Keplerian disk. We label magnetic field lines with the magnetic flux function $`\psi `$ defined by $`\psi /R=RB_z`$. Except for a few specific differences noted here and below, the notations and numerical methods we use are the same as those in Paper I, to which we refer for more details.

The governing equations are solved numerically using the code, subject to a set of appropriate initial and boundary conditions. We use a parallelized version of the code (Clarke, Norman, & Fiedler 1994), whose main algorithms are based on Stone & Norman (1992a, b) and Hawley & Stone (1995). The magnetic field is evolved by the well-known constrained transport method (Evans & Hawley, 1988), which preserves the condition $`𝐁=0`$ at all times as long as this condition is satisfied everywhere initially.

### 2.3. Boundary Conditions

On the four sides of our 2D simulation box (see Fig. 2), we need to impose boundary conditions. Numerically, these conditions are imposed by assigning values to flow quantities in a few ghost zones just outside the active computational grid. The ghost zone quantities include the density $`\rho `$, internal energy $`u`$, the three components of velocity $`𝐯`$, and the three components of the field $`𝓔𝐯\times 𝐁`$.

The conditions on the axis $`R=0`$ and the outer boundaries $`R=R_{\mathrm{max}}`$ and $`z=z_{\mathrm{max}}`$ are relatively straightforward: on the axis we demand reflection symmetry, whereas on the outer boundaries we allow for outflow. These conditions are standard options for the code. The conditions on the equatorial plane $`z=0`$ are more difficult to specify and implement. They represent the main technical improvement that distinguishes this work from others.

For our problem, the $`z=0`$ surface is divided into two distinct regions on which different boundary conditions are imposed (see Fig. 2). Inside the radius $`R_L`$, which divides these two regions, magnetic field lines are anchored and outflowing materials are injected. The boundary conditions there are the same as those used in Paper I. Specifically, the quantities $`\rho `$, $`B_z`$ and $`v_z`$ are specified as known functions of the radius $`R`$ at all times, and the two horizontal components of the $`𝓔`$ field satisfy $`_\varphi =0`$ and $`_R=v_\text{K}B_z`$. To implement these boundary conditions for $`𝓔`$, the horizontal components of $`𝐯`$ in the ghost zones are chosen so that $`v_R/v_z=B_R/B_z`$, and $`v_\varphi =v_\text{K}+B_\varphi v_p/B_p=v_\text{K}+B_\varphi v_z/B_z`$. In the ghost zones below $`z=0`$, we force the components of $`𝓔`$ to have values: $`_\varphi (z)=_\varphi (z)`$, $`_R(z)=2v_\text{K}(R)B_z(R)_\varphi (z)`$, and $`_z(z)=_z(z)`$. These conditions allow the foot point of a wind-launching open field line to be pinned firmly on the launching surface and at the same time leave the field line free to bend both radially and azimuthally in response to the stresses exerted by the wind material.

Our boundary conditions on the launching surface inside $`R_L`$ are similar to those used by Ustyugova et al. (1995, 1999) and Romanova et al. (1997), except that they also considered the case of subsonic wind injection. In the subsonic case, the wind must pass smoothly through a sonic point close to the launching surface, which fixes the mass loading rate and reduces the number of the quantities that must be specified at the $`z=0`$ surface by one. The sonic point is sensitive to the thermal and magnetic structure of the disk atmosphere, which is not well understood. We will postpone a detailed treatment of the sonic transition to a future investigation. Ouyed & Pudritz (1997) appear to have fixed all flow quantities in the ghost zones, which is different from our treatment. Fixing all quantities in the ghost zones risks creating discontinuities in the radial and toroidal components of the magnetic field between the ghost and active zones, since in the active zones these two components cannot be prescribed; they are determined by the stresses in the wind.

Outside the radius $`R_L`$, the equatorial $`z=0`$ surface is occupied by the last flux surface along which the wind material emanating from the radius $`R_L`$ flows. We demand that the vertical field component $`B_z=0`$ on this portion of the $`z=0`$ surface, and $`B_z(z)=B_z(z)`$ in the ghost zones below $`z=0`$. These conditions on the magnetic field are enforced by the conditions on the field $`𝓔`$, which are $`_\varphi (z)=_\varphi (z)`$, $`_R(z)=_R(z)`$, and $`_z(z)=_z(z)`$. Equation (3) then ensures that the initial value of $`B_z`$ at $`z=0`$ will not change. For the hydrodynamic quantities, we demand that $`\rho `$, $`u`$, $`v_R`$ and $`v_\varphi `$ be continuous across the $`z=0`$ surface and assign the value $`v_z=0`$ to the $`z=0`$ surface and the ghost zones below.

To complete our discussion of boundary conditions on the $`z=0`$ surface, we note the following three additional features. First, to fill out space, the field lines anchored in the wind launching region must become more and more horizontal as the edge of the region (at $`R=R_L`$) is approached. This condition is satisfied by forcing the functions defining $`v_z(R)`$ and $`B_z(R)`$ for $`(z=0,R<R_L)`$ to zero at $`R=R_L`$, keeping their ratio finite. Second, to prevent inflow in the region surrounding the rotation axis where the magneto-centrifugal mechanism fails, we inject a cold, light, fast axial flow capable of escaping the gravitational pull of the central compact object. Plausible physical justifications for this fast injection are given in Paper I. Finally, during the wind acceleration, there can appear temporary negative values of $`v_z`$ in active zones next to the launching disk. For this case only, we allow the boundary condition at $`z=0`$ to absorb the flow, by setting $`v_z`$ to zero at $`z=0`$, and allowing the other hydrodynamic quantities to have the same value in the ghost and active zones, but not altering the conditions on the $`𝓔`$ field. This procedure avoids some numerical defects, as explained in Paper I. The need for this special treatment disappears once a steady state near the disk surface is reached.

### 2.4. Initial Magnetic Field Distribution

Following Ouyed & Pudritz (1997), we choose as the initial (purely poloidal) magnetic field distribution above the disk a current-free field compatible with the function $`B_z(R)`$ defined at $`z=0`$. Such an initial field has the advantage of perturbing the ambient medium as little as possible. (This gas will eventually be pushed outside the computational domain by the fast moving wind launched from the disk.) One can solve the force free equation $`^2\psi =0`$ for the flux function $`\psi `$ through separation of variables, with the following set of boundary conditions: $`\psi (R=0,z)=0`$, $`\psi (R,z=\mathrm{})=0`$, $`\psi (R=\mathrm{},z)=0`$, and $`\psi (R,z=0)=\psi _D(R)`$. The function $`\psi _D(R)`$ is prescribed in the launching region inside the radius $`R_L`$; beyond $`R_L`$, it has a constant value $`\psi _L`$. A formal solution that satisfies these conditions is

$`\psi (R,z)`$ $`=`$ $`R{\displaystyle _0^{\mathrm{}}}𝑑ke^{kz}J_1(kR)g(k)`$ (9)
$`g(k)`$ $`=`$ $`k{\displaystyle _0^{\mathrm{}}}𝑑RJ_1(kR)\psi _D(R)`$ (10)
$`=`$ $`k{\displaystyle _0^{R_L}}𝑑RJ_1(kR)\psi _D(R)+J_0(kR_L)\psi _L`$ (11)
where $`J_0`$ and $`J_1`$ are Bessel function of zeroth and first order. Doing these integrations at each grid point can be expensive, especially in regions where the integrals are slow to converge, due to the rapidly oscillatory integrand. A more convenient and surer procedure to find the solution consists of first evaluating with sufficient approximation the integrals at the outer boundaries $`R=R_{\mathrm{max}}`$ and $`z=z_{\mathrm{max}}`$, and then using these values to solve an elliptic problem with Dirichlet boundary conditions, together with the known values of $`\psi _D`$ at the disk and axis. Since the size of the simulation box is much larger than the size of the launching region, one can expand the integrals in powers of $`k`$ at distances far from the launching region (i.e., $`rR_L`$). The first two terms of expansion, which proved sufficient to calculate our boundary values, are

$$\psi (R,z)=\psi _L\left(1\frac{z}{r}\right)+\frac{3CzR^2}{r^5},$$
(12)
where $`C=_0^{R_L}R\psi _D(R)𝑑R\psi _LR_{L}^{}{}_{}{}^{2}/2`$. We solved the elliptic equation by a relaxation method, taking as an initial guess the simple distribution $`\psi (R,z)=\psi _D(r)(1z/r)`$.

With the governing equations, numerical method, initial and boundary conditions explained, we are ready to examine in detail a representative wind solution that has reached a steady state.

## 3. An Illustrative Wind Solution

### 3.1. Simulation Setup

We first establish a “standard” run against which other simulations will be compared. All simulations are carried out with dimensionless quantities, but to facilitate comparison with high resolution observations that are becoming available, especially from the HST, we will present results with dimensional units. To fix units, we adopt one solar mass for the central object $`M_{}`$, and $`R_0=0.1\text{AU}`$ for the inner radius of the wind-launching region, which is roughly the corotation radius of the stellar magnetosphere (inside which the disk is either very sub-Keplerian or non-existent; Königl 1991; Shu et al. 1995). The gravitational softening radius $`r_g`$ is set equal to $`R_0`$. The central mass and inner radius set a scale for the velocity, which is the Keplerian speed at $`R_0`$, $`v_{\text{Ko}}=\sqrt{GM_{}/R_0}=94\text{km}\text{s}^1`$. The scales for the magnetic field strength, $`B_0`$, and the density, $`\rho _0`$, are related by the requirement that $`B_0^2/(4\pi \rho _0)=v_{\text{Ko}}^{}{}_{}{}^{2}`$. Both of these scales are fixed by demanding that the total mass flux of the outflow $`\dot{M}`$ (from each side of the disk) be $`\dot{M}_0=10^8\text{M}_{\mathrm{}}\text{yr}^1`$, a value representative for classical T Tauri stars (e.g., Edwards, Ray, & Mundt 1993).

For the standard run, we choose an outer radius for the wind launching region, $`R_L=1\text{AU}`$, ten times the inner radius $`R_0`$. On the equatorial plane $`z=0`$ inside $`R_L`$, one needs to specify the distributions of the density $`\rho `$ at the base of the wind, and the vertical components of both the injection velocity ($`v_z`$) and the magnetic field strength ($`B_z`$). We consider a single-component distribution for the magnetic field

$$B_z(R)=B_{z,\mathrm{max}}\left[1+\left(R/R_0\right)^2\right]^{e_b/2}\times f(R),$$
(13)
where the parameter $`B_{z,\mathrm{max}}`$ sets the scale for the field strength and the exponent $`e_b`$ the distribution. They are taken to be $`4B_0`$ and $`1.5`$ respectively for the standard run (where $`B_0`$ is the magnetic field unit defined above). The function $`f(R)`$ is set to unity between the origin and a transitional radius $`R_T`$, and

$$f(R)=\left(1+2\frac{R^2R_T^2}{R_{L}^{}{}_{}{}^{2}R_T^2}\right)\left(\frac{R_{L}^{}{}_{}{}^{2}R^2}{R_{L}^{}{}_{}{}^{2}R_T^2}\right)^2$$
(14)
between $`R_T`$ and the outer radius of the launching region $`R_L`$. It is used to force both $`B_z`$ and $`v_z`$ to vanish together at the outer edge of the launching region, as required by the boundary conditions (see §2.3). We choose $`R_T=0.8\text{AU}`$.

For the injection speed at the launching surface, we consider a two-component distribution, with

$$v_z(R)=\frac{v_{z0}}{1+\left(R/R_0\right)^2}$$
(15)
in the fast injection axial region inside $`R_0`$, and

$$v_z(R)=\alpha v_\text{K}\times f(R)$$
(16)
in the magneto-centrifugal wind launching region between $`R_0`$ and $`R_L`$. The parameter $`v_{z0}`$ sets the maximum injection speed on the axis, and is taken to be $`160\text{km}\text{s}^1`$ for the standard run. The parameter $`\alpha `$ denotes the ratio of the injection speed to the Keplerian speed, and is set to a small value of 0.1 for the standard run.

For the density distribution at the base of the wind $`z=0`$, we adopt a small, constant value of $`\rho _{\text{in}}=0.1\rho _0`$ (where $`\rho _0`$ is the density unit) in the fast injection region inside the radius $`R_0`$. On the Keplerian disk between $`R_0`$ and $`R_L`$, we prescribe a power-law distribution

$$\rho (R)=\rho _{\text{out}}\left(\frac{R_0}{R}\right)^{e_\rho },$$
(17)
up to the outer edge of the launching region $`R_L`$. The parameters $`\rho _{\text{out}}`$ and $`e_\rho `$ set the density scale and its decline rate with radius. They are taken to be $`\rho _{\text{out}}=10\rho _0`$ and $`e_\rho =1`$ for the standard run. Note the factor of $`10^2`$ contrast between the densities of the fast axial injection ($`\rho _{\text{in}}`$) and the inner edge of the magneto-centrifugal wind ($`\rho _{\text{out}}`$).

The launching conditions specified above yield the following units for the field strength and density: $`B_0=1.1\text{G}`$ and $`\rho _0=1.2\times 10^{15}\text{g}\text{cm}^3`$. These conditions are displayed in Fig. 3, together with a distribution of the accumulative mass flux injected into the wind and the magnetic flux enclosed within a given radius on the disk. It is clear that the fraction of mass flux in the fast injection region inside $`R=R_0`$ is small ($`1.5\%`$), despite the large injection speed in this region; the large speed is more than offset by the low density specified. The large velocity shear near $`R=R_0`$ does not appear to create any instability, possibly because of the strong (poloidal) magnetic field present in the region.

Above the equatorial surface $`z=0`$, we start with a purely poloidal, potential magnetic field, constructed using the method outlined in §2.4. The computational domain is initially filled with a low density medium, whose density distribution is prescribed either analytically or taken from other, similar runs that have reached a steady state. Judicious choice of the initial ambient density distribution can speed up the flow convergence to the steady state enormously, although the final wind solution does not depend on the choice, since the ambient medium will eventually be swept out of the simulation box completely. To isolate the magneto-centrifugal wind acceleration and collimation, the focus of our study, from thermal effects, we keep the flow cold at all times by arbitrarily setting a small value of $`0.02\text{km}\text{s}^1`$ for the sound speed.

Our main goal is to study the structure of magneto-centrifugal winds on large, observable scales. For this purpose, we adopt a simulation box of $`100\text{AU}\times 100\text{AU}`$. It is $`100`$ times the size of the wind launching region and $`1000`$ times the inner Keplerian disk radius $`R_0`$. Such a scale is directly accessible to high resolution observations using the HST and ground-based adaptive optics (Dougados et al., 2000) for nearby optically revealed, T Tauri winds. Outflows from deeply embedded sources can be probed on a similar scale using water masers. For example, the water masers of the Class 0 source S106 FIR appear to lie on a $`3\text{AU}`$ (width) $`\times `$ $`4\text{AU}`$ (length) U-shaped surface some $`25\text{AU}`$ away from the central source, possibly driven by a well-collimated “micro jet” Furuya et al. 1999). Numerically, we cover the active computation domain with $`190\times 210`$ zones. In the vertical direction, there are 40 uniform grid points covering the region from 0 to $`0.8\text{AU}`$, and 150 non-uniform zones covering the region from 0.8 to $`100\text{AU}`$, with the zone size increasing by a small factor of 1.035023 between adjacent zones when moving outwards. In the horizontal direction, the first 60 grid points from 0 to $`1.2\text{AU}`$ (thus covering the entire launching region) are uniform, followed by 150 non-uniform grid points with approximately the same ratio between zones as before (1.034987). Larger boxes are in principle possible, but limited by computer time.

### 3.2. Global Structure, Flow Collimation, and Jet Formation

We follow the evolution of the material injected into the computation domain numerically with our modified code, keeping the boundary conditions shown in Fig. 3 fixed at all times. The injected wind material pushes aside the low-density ambient medium and propagates towards the outer edges of the simulation box. After all of the ambient medium is swept out of the box, the wind settles quickly into a steady state, with no appreciable time variations. The overall appearance of the steady-state wind is given in Fig. 4, on two scales. In Fig. 4a, we show on the $`10\text{AU}`$ scale a meridian view of nine streamlines (also field lines) that divide the wind into ten zones of equal mass flux. Superposed on the streamlines are the isodensity contours (in shades), with values of density decreasing exponentially outwards. A salient feature is the fast magnetosonic surface, shown in dashed lines in Fig. 4a. As mentioned earlier, the fast surface closes inside the simulation box, except near the axis where the light, fast axial flow resides. The tenuous flow in this narrow axial region remains dominated by the (largely poloidal) magnetic field, which can, in principle, provide some stability to the wind as a whole against kink instabilities (Shu et al. 1995). The axial injection occupies a decreasing fraction of the total volume on larger scales. It becomes nearly invisible in panel (b) of Fig. 4, where streamlines and isodensity contours are plotted on the $`100\text{AU}`$ scale. Here, most of the space is filled with the flow driven from the disk magneto-centrifugally, whose acceleration and collimation properties we want to study in detail.

We first concentrate on the collimation properties of the standard solution. Self-collimation is commonly accepted as a hallmark of the magneto-centrifugal winds. The degree of collimation is still a matter of debate, however (see, e.g., Okamoto 1999). Our self-consistently determined solution up to large distances can shed light on this important issue. From panel (a) of Fig. 4, we find that the standard wind is not well collimated on the small scales close to the launching surface, although hints of gradual collimation are evident in the shapes of both streamlines and density contours: the streamlines bend slightly towards the rotation axis, which forces the density contours to elongate along the axis. The elongation of density contours shows up more clearly on the larger scale in panel (b) of Fig. 4, with the contours appearing nearly cylindrical in the region within $`15\mathrm{°}`$ of the rotation axis. Asymptotically cylindrical density stratification was predicted by Shu et al. (1995) analytically. The cylindrical shape does not extend all the way to the equator, however, at least on scales up to $`100\text{AU}`$. The density contours in the equatorial part of the wind bulge outwards. The bulging is similar to that found numerically by Sakurai (1987; see also Ouyed & Pudritz 1997). The relative shape of the bulge appears to change little from one isodensity contour to the next, which is probably related to the logarithmically slow collimation of streamlines in the equatorial region at large distances (see Bogovalov & Tsinganos 1999 for a recent discussion). It will be slow to disappear, if ever, on even larger scales.

The shape of isodensity contours on large scales is important to compute because it is intimately related to the observed shape of optical jets. A density stratification bulging out prominently near the base as in Fig. 4 may be in conflict with high resolution observations of some HH jets. In particular, the spectacular HST images of HH 30 (Burrows et al. 1996; Ray et al. 1996) show that this jet remains nearly cylindrical all the way to the disk surface; if anything, it appears to be slightly pinched in the equatorial region, keeping in mind, however, that contamination of the jet emission from scattered stellar light is a concern near the disk. Since the emissivity is sensitive to density, the optical appearance of a jet should follow to a large extent that of isodensity contours (Shang et al., 1998), although the thermal structure affecting the emission is uncertain. The requirement that the density contours be nearly cylindrical all the way to the disk surface (assuming the likely situation of a much less rapid variation in temperature than in density) imposes a constraint on possible combinations of magnetic field and mass flux distributions on the launching surface. The standard run does not appear to satisfy this constraint.

Another generic feature predicted by the asymptotic analysis of Shu et al. (1995) is the strong stratification in density transverse to the rotation axis. This feature shows up clearly in the density contours. We stress that the stratification is a most important characteristics of the magneto-centrifugal wind, even though individual streamlines appear hardly bent at all. We estimate that, at a height of $`90\text{AU}`$, the density distribution decreases rapidly outside a small core of $`2\text{AU}`$, roughly as a power-law of the cylindrical radius $`R^{1.3}`$. The decline is slower than the asymptotic estimate of $`R^2`$ of Shu et al. (1995), obtained assuming a constant flow speed on all field lines; the difference is partly accounted for by the lower speed achieved along the field lines in the more equatorial part of the standard solution (Matzner & McKee, 1999). Nevertheless, the fact that the density decline is significantly faster than $`R^1`$ means that the central high-density “jet” will stand out in maps of integrated emission along the line of sight. The steep density stratification, coupled with a moderate anisotropy in flow velocity, gives rise to a strong concentration of energy and momentum fluxes in the direction along the axis, with implications on the dynamics of bipolar molecular outflows (Li & Shu, 1996), which are probably formed when the highly anisotropic winds run into the ambient medium. The decrease in density away from the axis is qualitatively consistent with the HST observations of the wind of DG Tau by Bacciotti et al. (2000).

Even though a dense, nearly cylindrical “jet” is naturally produced around the axis, at least well above the disk, the fraction of the wind mass flux that gets collimated into the jet is a concern. To fix ideas, we adopt a fiducial number density of $`n_\text{H}=10^4\text{cm}^3`$ to delineate the outer boundary of the central “jet”. From Fig. 4b, we find that only $`20\%`$ of the wind mass flux resides in the fiducial jet at the height of $`z=10^2\text{AU}`$; the majority remains in the wide-angle component. Such a jet formation efficiency may be too low, at least for the famous HH 30 system, as discussed in the introduction. This example illustrates that there is no guarantee that a disk-driven magneto-centrifugal wind would automatically (a) concentrate most of its mass flux into an axial jet or (b) produce a density stratification cylindrical all the way to the equator.

The relatively low efficiency of jet formation and undesirable shape of the density contours of the standard solution are related to the particular distribution of mass flux chosen. It is dominated by the mass injected at the outer part of the launching region, since the mass flux per octave in radius is

$$\frac{d\dot{M}}{d\mathrm{ln}R}R^{1/2},$$
(18)
between the inner radius of the Keplerian disk $`R_0`$ and the radius $`R_T`$ (beyond which the injection speed is modified by equation ). Mass from the outer part is loaded onto more horizontally inclined field lines, which are harder to collimate. The above deficiencies are rectified to some extent in the wind solutions to be discussed in the next section. Moderately collimated winds like the standard solution may be applicable to the outflows from young high-mass stars, which could be more dominated by the wide-angle component than their low-mass counterparts (Richer et al., 2000).

We note that the standard solution (and other solutions to be discussed below) does not recollimate towards the rotation axis at large distances, and does not show any oscillations along the axis. Both behaviors are seen in some self-similar (Chan & Henriksen 1980; Blandford & Payne 1982; Achtenberg, Blandford, & Goldreich 1983) and width-averaged jet solutions (e.g., Spruit, Foglizzo, & Stehle 1997). They are most likely a result of the self-similarity assumption or the averaging over the jet cross-section, not the general properties of magneto-centrifugal winds.

### 3.3. Flow Acceleration and Kinetic-Energy Domination

We now concentrate on the acceleration, rather than collimation, of the standard solution. Magneto-centrifugal acceleration along field lines is illustrated in panel (a) of Fig. 5, where the poloidal velocity is plotted as a function of (log) spherical radius along the four streamlines that divide the wind into five zones of equal mass flux. As expected, the flow along each field line crosses the fast point smoothly, getting accelerated to a speed $`23`$ times the Keplerian speed at the foot point of the field line. Note that for this particular solution, the flow moves faster along the streamlines in the polar region than in the equatorial region, by a maximum factor of $`4`$. This pole-equator velocity contrast is qualitatively similar to that deduced for the wind of DG Tau from HST observations (Bacciotti et al., 2000). These authors inferred that the high velocity component of forbidden line emission of DG Tau wind occurs near the axis, and is spatially bracketed by the lower velocity component. Whether wind solutions like our standard run can explain the density distribution and velocity field observed quantitatively remains to be seen.

At large distances from the launching region, most of the energy extracted magnetically from the disk, as measured by the Poynting flux, is expected to be converted into the kinetic energy of the ordered plasma flow. This kinetic domination is indeed seen in the standard solution, and is illustrated in Fig. 5b, where contours of constant fast (magnetosonic Mach) number, $`M_\text{f}=v_p/v_\text{f}`$ (where $`v_p`$ and $`v_\text{f}=[(B_p^2+B_\varphi ^2)/(4\pi \rho )]^{1/2}`$ are the poloidal and fast magnetosonic speed) are plotted. Note that the fast number along each field line reaches a maximum value between $`1.6`$ and $`3.0`$ at the $`10^2\text{AU}`$ distance from the origin. Since the ratio of Poynting and kinetic energy fluxes is approximately $`2/M_{\text{f}}^{}{}_{}{}^{2}`$ at large distances (Spruit, 1996), some $`60\text{}80\%`$ of the total energy is kinetic. Similarly high (low) kinetic (magnetic) energy fraction has been found numerically by Ouyed & Pudritz (1997) and Ustyugova et al. (1999). This efficient conversion is of interest because the remaining, low magnetic energy makes potential destruction of toroidal magnetic field by kink instabilities (Spruit 1996; Begelman 1998) less disruptive to the wind; the wind has become basically ballistic by the time the outer edges of the simulation box are reached, particularly in the region not far from the polar axis where the fast number is the highest. The transition of the flow from magnetic to kinetic domination is an esthetically pleasing feature of the magneto-centrifugal wind. The asymptotic kinetic domination is in contrast with that of the “classical” 1-D Weber-Davis (1967) wind solution, where the kinetic energy flux is at most a third of the total (Spruit, 1996). It also differs from the ultra-relativistic case where kinetic domination is difficult, if not impossible, to achieve in an ideal, steady MHD wind that fills the entire $`4\pi `$ steradians (Chiueh, Li, & Begelman 1998; Okamoto 2002).

We should emphasize that even though the fast number $`M_\text{f}`$ is greater than unity, it is not greater by much; it has values of only a few on the $`100\text{AU}`$ scale, as seen from Fig. 5b. In other words, the magnetic field, even though its energy fraction is small, is not completely negligible on large, observable scales. In particular, the relatively modest values of $`M_\text{f}`$ make the effects of magnetic cushion important in the interaction of the magnetized wind with its ambient medium (or with itself if the wind is variable) through shocks, especially if the shocks are oblique.

## 4. Jets with Better Collimation and Higher Formation Efficiency

In the standard wind solution discussed in the last section, only a relatively small fraction of the wind mass flux resides in the fiducial “jet”, and the isodensity contours bulge out in the equatorial region. Both features are undesirable for modeling jets like HH 30. In this section, we seek to rectify the situation, by concentrating more mass loading at smaller disk radii, through a steeper power-law distribution for the density at the launching surface while keeping the magnetic field distribution fixed. The density distribution is controlled by the exponent $`e_\rho `$ in the equation (17). The standard run has a relatively flat density distribution specified by $`e_\rho =1`$. We have varied this exponent over a range of values, and found that the collimation properties of the wind can change substantially from the standard run. We choose a case with $`e_\rho =3`$ to illustrate the changes.

For the $`e_\rho =3`$ run, we keep the distributions of the magnetic field and injection speed on the disk the same as in the standard run, and set the density at the inner edge of the Keplerian disk $`\rho _{\text{out}}=10\rho _0`$ as before. Since the density now declines more rapidly with radius, the (dimensionless) mass flux ejected from the Keplerian disk is reduced. We correspondingly reduce the density $`\rho _{\text{in}}`$ at the base of the fast injection region by a factor of 4 (to $`\rho _{\text{in}}=0.025\rho _0`$), so that its mass flux fraction remains small ($`2\%`$ of the total). The vast majority of the wind material is launched magneto-centrifugally from the surface of the Keplerian disk, according to the distribution

$$\frac{d\dot{M}}{d\mathrm{ln}R}R^{3/2}.$$
(19)
Nearly half of the mass flux is concentrated at the inner edge of the Keplerian disk, between $`1R_0`$ and $`1.5R_0`$ (or $`0.10`$ and $`0.15\text{AU}`$). The concentration is shown clearly in panel (a) of Fig. 6, where nine streamlines (also field lines) dividing the wind into ten zones of equal mass flux are plotted, after a steady state is reached. Most of the streamlines emanate from close to the inner Keplerian disk, creating an impression of an “X-wind”, although most of the wind-launching magnetic flux lies outside the region where the mass flux is concentrated (see Fig. 3b).

The density contours in the $`e_\rho =3`$ case started out more collimated in the launching region than those in the standard run. On the $`10^2\text{AU}`$ scale shown in panel (b) of Fig. 6, they become nearly cylindrical all the way to the equator; if anything, the contours pinch slightly inwards in the equatorial region, in contrast with the bulging observed in the standard solution. Thus, a better collimated jet is produced.

The efficiency of jet formation is also improved. From Fig. 6b, we estimate that $`45\%`$ of the wind mass flux is enclosed, at the height of $`100\text{AU}`$, within the density contour of $`10^4\text{cm}^3`$, which marks the boundary of the fiducial jet. In contrast, the jet fraction is only $`20\%`$ for the standard run. Therefore, compared with the standard solution, the wind is more dominated by the axial jet, although there is still more mass flux in the wide-angle component outside the jet than inside the jet itself. For the chosen size of the Keplerian disk (one decade in radius, from $`0.1`$ to $`1.0\text{AU}`$) and magnetic field distribution ($`B_z`$ roughly proportional to $`R^{3/2}`$), we find it difficult to collimate more than half of the mass flux into the jet, at least on the $`10^2\text{AU}`$ scale. The jet fraction may be increased by external collimation, which can act in addition to the intrinsic self-collimation. This will be particularly true during the earliest, Class 0 stage of star formation, when a massive envelope (André, Ward-Thompson, & Barsony, 2000) is present, which can confine the wind. The combination of external and intrinsic collimation may produce a more jet-dominated magneto-centrifugal wind; such a wind may be needed to explain the kinematics of some bipolar molecular outflows (Lee et al., 2000). We will address the issue of external collimation in the future.

We note that the cylindrically-stratified density decreases rapidly away from the rotational axis, roughly as $`R^{2.2}`$ at a height of $`90\text{AU}`$ in this example. The decline is steeper than in the standard case, where $`\rho R^{1.3}`$ roughly. It is not far from $`R^2`$, perhaps signaling that the wind has reached an asymptotic state close to that analyzed by Shu et al. (1995). We note that the asymptotic wind speed in this particular example is almost independent of the polar angle, showing little signature of the Keplerian rotation at the launching surface. It illustrates the obvious point that a disk-wind does not have to be slower in the equatorial region; it depends on the distribution of mass loading for a given magnetic field configuration.

As a final remark, we stress that the mass loading in the axial region where the magneto-centrifugal mechanism fails is somewhat artificial. Nevertheless, it provides a relatively non-intrusive inner boundary to the magneto-centrifugally driven outflow. It has the added advantage of possibly representing a fast stellar coronal wind emanating along open stellar field lines; such a wind may have a predominantly poloidal magnetic field, which could add to the stability of the entire flow (Shu et al., 1995). We have done several tests to ensure the non-intrusive nature of the fast injection, including increasing by a factor of 4 the density in the injection region for the two examples discussed in this and last sections. The flow structure in the magneto-centrifugal region outside the axial injection remains little affected, especially at large distances. Even though the axial injection appears to have minimal effects on the region of interest, it controls the timestep of the simulation through the Courant condition, because of its low density and large magnetic field strength (and thus high Alfvén speed). To explore parameter space in a reasonable amount of time, as we do in a companion paper (J. Anderson et al. 2003, in preparation), one can in principle adopt a relatively heavy axial injection.

## 5. Summary and Future Work

We have developed a -based numerical code capable of following the acceleration and collimation of magneto-centrifugal winds from their launching surface to large, observable distances. The code is made possible by treating the wind-driving accretion disk as a boundary, which allows for a clear separation of the dynamics of the disk and the wind. By limiting the launching region to the inner part of the disk, we were able to follow the outflow up to, and well beyond, the fast magnetosonic surface for most of the streamlines. Close to the rotation axis where the magneto-centrifugal mechanism fails, a light, fast jet is injected, which could stay sub-fast magnetosonic in the computational domain; it occupies an increasingly small fraction of space at large distances, and provides a non-intrusive inner boundary to the super-fast magneto-centrifugal wind that we are interested in.

We find that, as expected, most of the magnetic energy extracted from the disk can be converted into the bulk flow kinetic energy at large distances, producing an essentially ballistic wind. Its fast magnetosonic Mach number remains moderate, however, which has implications for wind interaction. Our unique simulation setup allows us to investigate quantitatively the large-scale wind structure, particularly jet formation. In agreement with the asymptotic analysis of Shu et al. (1995), we find nearly cylindrical stratification in the wind density close to the axis at relatively large heights. Closer to the equator, the density contours can either bulge outwards or pinch inwards, depending on the launching conditions. The fraction of the wind mass flux residing in the dense axial jet can vary substantially, again depending on the conditions on the disk, particularly the distribution of mass loading rate for a given magnetic field configuration. It appears possible to have, at least on the $`10^2\text{AU}`$ scales around T Tauri stars, both winds dominated by the wide-angle component and more jet-dominated winds. The exact demarcation of these two regimes on the $`10^2\text{AU}`$ and larger scales is unclear. It will require a time-consuming, systematic exploration of parameter space.

Our robust determination of the large-scale structure of magneto-centrifugal jets and winds opens up several new possibilities. First, it avoids the artificial instabilities seen in some of previous simulations where a major fraction of the wind leaves the simulation box with a sub-fast speed. However, non-axisymmetric instabilities of physical origin, such as the kink instability, are possible. They will be examined in a subsequent paper in this series with 3D calculations (R. Krasnopolsky et al. 2003; Paper III) using axisymmetric solutions as a starting point. Second, the size of our wind launching region is adjustable. It allows us to study both the nominal “disk-wind” (e.g. Königl & Pudritz 2000) and the “X-wind” (Shu et al., 2000), which differ mainly in the width of the launching region. The effects of the size of the launching region and other parameters, such as the mass loading rate, on the large-scale wind structure will be explored in a companion paper (J. Anderson et al. 2003, in preparation). Finally, by adopting a time-dependent launching condition, we will be able to study the formation of internal shocks, which may be used to interpret the knots commonly observed in YSO jets.

We thank J. Anderson for helpful discussion and C. Gammie for computational facilities used in the production runs. The work is supported in part by NASA grants NAG 5-7007, 5-9180, 5-12102 and by NSF grant AST 00-93091.


format: text

text:
# Anisotropic dynamics in a shaken granular dimer gas experiment

## Abstract

The dynamics, velocity fluctuations, and particle-plate interactions for a 2D granular gas of shaken, non-spherical particles are studied experimentally. The experiment consists of a horizontal plate that is vertically oscillated to drive the dynamics of macroscopic dimers, spherical pairs that are loosely connected by a rod that couple the interaction each of the spheres has with the shaking plate. The extended nature of the particles results in more than one energy-momentum transfer between the plate and each dimer per shaking cycle. This complex interaction results in anisotropic behavior for the dimer that is a function of the shaking parameters.

preprint: APS/123-QED

Recently, a novel driven granular gas experiment comprised of two layers of separate species of granular material has demonstrated robust velocity fluctuations that are nearly Gaussian over a wide range of parameters Baxter and Olafsen (2003). In the experiment, a horizontal plate vertically drives a first layer of one species of particles whose collisions, in turn, fluidize a second layer comprised of a different species of particles. While the velocity fluctuations in the first layer, driven by the plate, remain strongly non-Gaussian and vary with the driving parameters, the second layer, thermalized by collisions with the first layer, demonstrate nearly Gaussian velocity statistics over a wide range of shaking parameters. Similar to this mechanically fluidized bed, a gas fluidized bed system used to thermalize the motion of a larger sphere also demonstrates nearly Maxwell-Boltzmann statistics for the indirectly driven particle Ojha et al. (2004). Clearly, the way in which energy is injected into each layer in the mechanically fluidized experiment is an important detail that must be examined more closely to understand how each layer could simultaneously demonstrate drastically different velocity fluctuations.

In this report, we detail the injection of energy between the shaking plate and the species that comprise the lower layer of the two layer experiment: granular dimers, spherical pairs loosely connected by a short rod. Experiments to probe the dynamics of driven granular gases, collections of large numbers of dissipative macroscopic particles whose motion is maintained through various forms of external driving, have focused largely on species of identical spheres Olafsen and Urbach (1999a); Losert et al. (1999); Blair and Kudrolli (2001); Aranson et al. (2000); Painter and Behringer (2000); Rouyer and Menon (2000). However, identical spheres represent only one special species of macroscopic granular particles. We demonstrate that many of the differences in the dynamic behavior from that of identical spheres can be understood through anisotropies that directly relate to the geometry of the extended particles. In particular, the manner in which the particles are thermalized by the shaking plate in the vertical direction is intimately connected to how the horizontal dynamics are manifested in the plane of the 2D granular gas. A prior simulation has investigated developing a kinetic theory for such gases in a freely cooling case Huthman et al. (1999).

The experiment consists of an aluminum plate of radius 14.6 cm that is vertically oscillated by an electromagnetic shaker. The plate is level and flat, and the acceleration of the plate, $`\mathrm{\Gamma }=A(2\pi \nu )^2`$, shaking with peak amplitude, A, at frequency, $`\nu `$, is uniform, having a spatial variance across the surface that is less than 1$`\%`$. The particles are comprised of two hollow spheres that are each 3.2 mm in diameter and are loosely connected by a thin rod that allows a spacing between the spheres of zero to 1.6 mm. The total length of the dimer pair can therefore vary from 6.4 mm to 8.0 mm. Except at highest densities where excluded volume effects compress the dimers along the interconnecting rod, the motion of the bouncing on the plate is observed to keep the dimers fully extended.

The dimers are constructed by cutting pairs from chains similar to those used in prior investigations of the motion of chains Ben-Naim et al. (2001); Prentis and Sisan (2002); Belmonte et al. (2001). The loose connection of the rod between the two spheres of the dimer allows for an additional degree of movement where the spheres may rotate about the ends of the rod as wheels on an axle, letting the dimer roll in a direction perpendicular to the connecting rod.

The average mass of a single dimer is approximately 185 mg but can vary up to $`\pm `$ 7 mg from dimer to dimer. The aspect ratio of the dimer, which is the ratio of the length (as measured along the direction of the rod connecting the spheres) to the width, varies from 2 to 2.5 depending upon whether the dimer is fully compressed or fully extended. A coefficient of restitution is determined experimentally for the dimers to be approximately 0.3 coe . For comparison, the identical technique yields a value of 0.5 for one of the single hollow balls from a dimer. A delrin ball is measured with the same technique to have a value of 0.9.

The rod connecting the two spheres is a natural geometric axis for the purpose of describing the motion of the particles. When viewed from above, the two dimensional motion of the dimer in the horizontal plane can be considered in terms of two angles, as shown in Figure 1. Because of the extended nature of the particle, the dynamics can be decomposed into the motion of the center of mass, $`\theta `$ and $`\mathrm{\Delta }𝐫`$, and the motion about the center of mass, $`\varphi `$.

The motion of a single dimer on a vertically shaken horizontal plate has been measured using high-speed digital photography dal . By overlaying two sequential images, a composite picture similar to Figure 1 can be produced. The locations of the two spheres that comprise the dimer can be identified using software analysisidl . Identifying the location of the spheres that constitute the dimer in each frame allows the center of mass and the orientation of the dimer to be determined. If the orientation of the dimer relative to a fixed frame (such as that of the camera) is Di and Di+1 in frame i and i+1, respectively, then the angular displacement about the center of mass $`\varphi `$ is determined from the following relationship:

$$sin(\varphi )=\frac{𝐃_i\times 𝐃_{i+1}}{𝐃_i𝐃_{i+1}}.$$
(1)
In a similar fashion, the direction of the motion of the center of mass can be determined by examining the motion of the center of mass, $`\mathrm{\Delta }`$r, from one frame to the next relative to the director of the dimer, Di:

$$cos(\theta )=\frac{𝐃_i\mathrm{\Delta }𝐫}{𝐃_i\mathrm{\Delta }𝐫}.$$
(2)
Probability distributions of $`cos(\theta )`$ for a single dimer on a plate for different shaking frequencies at a peak plate acceleration of 2g are shown in Figure 2. The particle was tracked from above with a high speed camera and the results demonstrate a preference for the center of mass of the dimer to move along the direction of the rod connecting the two spheres at low frequencies. The reader should note the change in this bias by the decrease in the probability at $`cos(\theta )=\pm 1`$ as the shaking frequency is increased. There is an interesting exchange between the motion of the center of mass and about the center of mass as the frequency is increased. At 50 Hz, when the motion about the center of mass is biased toward the ends, the center of mass motion demonstrates Gaussian velocity statistics. As the frequency of shaking is increased to 70 Hz and the motion about the center of mass becomes more uniform, the velocity statistics for the center of mass become slightly non-Gaussian. At 90 Hz, where the motion about the center of mass is now biased for motion perpendicular to the director, the velocity statistics for the center of mass motion is strongly non-Gaussian. To understand this anisotropy in the dimer motion, it is necessary to examine how the single dimer interacts with the vertically shaking plate.

Figure 3 is an example of the signal acquired from the accelerometer to measure the peak acceleration of the shaking plate in this experiment. The data shown is for one typical oscillation of the plate with a peak acceleration of 2g at a frequency of 50 Hz. The signal is comprised of two parts. The sinusoidal wave is the measured acceleration of the plate. The two “noisy” portions of the signal (one near the trough and one near the peak of the oscillation) are the pings caused by each of the balls of the dimer colliding with the plate. The trough (negative acceleration) ping corresponds to a ball of the dimer colliding with the plate as it accelerates downward (see also Figure 4(a)). The peak ping corresponds to the other ball of the dimer colliding with the plate as it accelerates upward (Figure 4(c)). The bias in Figure 2 can be understood in light of this dimer-plate interaction in the vertical direction in terms of the individual collisions of each sphere of the dimer.

The motion of a single sphere interacting with a shaking plate has been studied previously Mehta and Luck (1990). In carefully examining the motion of the single sphere throughout the shaking cycle, regions of the system’s phase space were observed to either increase or decrease the sphere’s momentum. These regions were deemed transmitting and absorbing, respectively Mehta and Luck (1990), and were responsible for the general “chattering” up and down in the amplitude of the sphere’s motion over several shaking cycles. The analysis for the motion of the dimer may be interpreted for two spheres colliding at nearly opposite phases of the plate cycle.

When the one ball of the dimer collides with the plate during its downward acceleration, the dimer loses momentum (as would a single sphere in that phase of the shaking cycle Mehta and Luck (1990)) in the direction along the rod (Figure 4(a)). The first ball of the dimer recoils from this “soft” collision with the plate with a reduced velocity as the other ball of the dimer, under the influence of gravity, begins to fall toward the plate (Figure 4(b)). When the plate is accelerating upward, the dimer gains a momentum kick from the second sphere colliding with the plate (Figure 4(b)), as would a single sphere under the same circumstances Mehta and Luck (1990). That momentum kick is transmitted up the loose connection of the rod of the dimer. The result is a net momentum transfer along the rod of the dimer. This net interaction with the shaking plate results in a bias of the particle to move in a direction along the geometric axis of the dimer. The biasing effect grows at fixed acceleration as the frequency is decreased because the amplitude of shaking increases, increasing the vertical angle the dimer makes relative to the horizontal plane and emphasizing the bias between the two collisions. As the amplitude of shaking decreases, this angular rotation decreases and the effect is minimized but never completely disappears (see Figure 2). It should also be noted that this net momentum transfer led to a second important difference from the dynamics of monomers Olafsen and Urbach (1999a, b). In the monomer case, clustering could be avoided within the system by shaking at an acceleration as little as 1.25g for some values of A and $`\nu `$. Because of the more complicated momentum transfer between the dimers and the plate, accelerations of $`4g`$ were necessary to produce a nearly uniform density gas of dimers, i.e. to avoid clustering.

For more than one dimer on the plate, the accelerometer signal becomes more complicated and it is difficult to determine which particle is interacting with the plate. However, in a future paper we will demonstrate a method of describing the net momentum transfer for several dimers interacting with the plate to show that collisions occur more uniformly throughout the shaking cycle than in the monomer case. Our goal in the current work is to describe the velocity statistics of an inelastic granular gas comprised of many dimers (800 - 3000) as measured via high speed photography from above the plate for motion in the horizontal plane. The results for some parameters are shown in Figure 5. In nearly all regimes, the velocity distributions demonstrate non-Gaussian behavior similar to that of the monomer experiment of Olafsen and Urbach Olafsen and Urbach (1999a, b). It is interesting to note, however, that at the lower density, the velocity distribution function is nearly Gaussian. This can also be understood in light of the complex interaction the particles have with the plate.

The density of the dimer gas, or coverage of the plate, is calculated in terms of the equivalent number of dimer balls that would be needed to have a single layer on the shaking plate. This number is 7320 spheres or 3660 dimers (assuming full compression of all dimers along their connecting rods for one full layer). The density is then calculated as a fraction of the number of dimers in the cell normalized by 3660.

At low densities, the dimers do not undergo as many collisons with the other particles in the gas and the horizontal motion is dominated by local surface roughness of the plate, resulting in the nearly Gaussian distributions in the horizontal direction. However, the dimers still interact, so the motion cannot completely be thought of as that of a single isolated dimer as was the case in Figure 2. Indeed, at higher densities, the horizontal motion is dominated by dimer-dimer collisions, which in general can be thought of belonging to one of three types. Because each dimer has two interactions with the shaking plate as previously described, each dimer can be thought of as a pair of loosely coupled spheres, one for which there was a net momentum gain, $`p_+`$, and one for which there was a net momentum loss, $`p_{}`$, in its most recent interaction with the shaking plate.

Collisions between any two dimers $`i`$ and $`j`$ are therefore interactions of their constituent spheres that have most recently both received positive momentum kicks, $`<p_+^ip_+^j>`$, both received negative momentum kicks, $`<p_{}^ip_{}^j>`$, or two spheres that have received one of each type of kick from interacting with the plate, $`<p_+^ip_{}^j>`$ or $`<p_{}^ip_+^j>`$. The way in which energy and momentum are transferred from the vertical to the horizontal is therefore much more complicated to model. In general, as the density of the dimer gas increases, the velocity statistics become more non-Gaussian. However, because of the additional complexity of dimer collisions, it is not surprising that the overall shape of the deviations from Guassian statistics do not mimic the monomer case, where the shape of the entire distribution could be “tuned” from nearly Gaussian to nearly exponential depending on the shaking parameters and how “two dimensional” the system was constrained to be Olafsen and Urbach (1999b). The other general trend observed in the data is that the deviation from Gaussian behavior increases as the root mean square velocity of the inelastic gas decreases.

In conclusion, in this paper we detail similarities and differences in the observed dynamic behavior of a driven granular gas composed of dimers (linked spheres) from that of granular gases with monomer species. The dynamic differences are directly related to the geometric anisotropies of the constituent particles and the additional complexity with which these particles interact with the shaking plate. The results clearly underscore the need to better understand how energy and momentum are injected in such systems in order to study the dynamics of driven granular gases.

###### Acknowledgements.

This work was supported by a grant from the Petroleum Research Fund of the American Chemical Society and a grant from the General Research Fund of the University of Kansas. One of us (JA) was also partially supported by an Undergraduate Research Award (UGRA) from KU. The authors would like to thank G. William Baxter for several helpful conversations and feedback in the preparation of this manuscript.