text:
# Gap Formation by Planets in Turbulent Protostellar Disks

## 1 Introduction

Understanding the process of planet formation in nascent solar systems is a long-standing goal of astrophysical theory. The traditional picture is an orderly one in which planets slowly build up their mass first by accreting rocky planetesimals, then (if sufficiently massive) gas from the surrounding protostellar disk. The recent multiple detections of extrasolar planets close to the central star, a configuration once thought to be highly improbable, pose a stern challenge to this relatively simple picture of planet formation. In the wake of these discoveries, contemporary theories emphasize the importance of the dynamical interaction between a developing planet and the ambient gaseous disk.

A possible explanation for the presence of giant planets near their central star is that they have migrated inwards after having formed farther out in the disk. This requires angular momentum to be removed from the planet by an imbalance in the positive and negative torques exerted upon it (Goldreich & Tremaine, 1980; Ward, 1986, 1997). Two mechanisms have been identified, known as Type I and Type II migration. Type I migration is essentially a linear process in which the planet excites trailing density waves at the Lindblad resonances. The waves carry angular momentum outward through the disk. There are torques of opposite signs exerted at both the inner and outer resonances, but for a given $`m`$ (azimuthal wavenumber), the outer resonance lies closer to the planet and the coupling is stronger (Ward, 1997). Therefore, if the interaction is dominated by the Lindblad torques, then more angular momentum is lost than acquired by the planet, which would then drift inward. The role of corotation torques, however, is still under investigation, and needs to be understood if Type I migration is to be put on a firm theoretical footing (Tanaka, Takeuchi, & Ward, 2002). In Type II migration, nonlinear gravitational interactions between the disk and the planet clear a gap at the orbital radius of the planet (Lin & Papaloizou, 1979, 1993; Goldreich & Tremaine, 1980). Disk material then no longer accretes across the planet’s orbital radius, and the planet is carried inward with the disk itself on an accretion time scale. The rate of migration thus depends on whether a gap forms or not, and whether there is, in fact, an turbulent stress to drive accretion in the disk. Unless for some reason disk material remains abundant in the immediate vicinity of the planet, one important consequence of gap formation is likely to be that planet formation will be severely limited.

Traditionally, turbulence has been modeled as a Stokesian viscosity in accretion disk studies, and the criteria for a planet to open a gap are not the same in viscous and inviscid models. Lin & Papaloizou (1993) and Bryden et al. (1999) discuss the conditions pertinent to each class of model. In the case of an inviscid accretion disk, the competing forces are the local fluid pressure, which resists the creation of a gap, and the planet’s gravitational tidal force, which promotes gap formation. The tidal criterion for gap formation is that a gap can form if the planet’s Roche (or Hill) radius $`R_R`$ exceeds the pressure scale height $`H`$, which we define as the ratio of isothermal sound $`c_s`$ to local orbital frequency $`\mathrm{\Omega }`$. For a planet of mass $`M_p`$ orbiting a central star of mass $`M_{}`$ at radius $`R_p`$, the Roche limit is $`R_R=(M_p/3M_{})^{1/3}R_p`$. Thus a gap forms if

$$q\frac{M_p}{M_{}}3\left(\frac{H}{R_p}\right)^3.$$
(1)
Denoting the right hand side of equation (1) as $`q_t`$, the tidal criterion is simply $`qq_t`$.

In viscous disk models, there is another condition that must be satisfied before a gap can open. The is set by the requirement that the rate at which planetary torques clear the gap must exceed the rate at which viscous torques replenish it. The condition is (Lin & Papaloizou, 1993)

$$qq_v40\alpha \left(\frac{H}{R_p}\right)^2,$$
(2)
where $`\alpha `$ is defined as the ratio of the viscosity to $`c_sH`$. This is essentially the Shakura & Sunyaev (1973) parameterization

These gap conditions have been extensively investigated via two-dimensional numerical simulations of $`\alpha `$ disks (Bryden et al., 1999; Kley, 1999; Lubow, Seibert, & Artymowicz, 1999). An important limitation of this approach is that the most interesting consequences of the planet-disk interaction depend upon the nature of the stress, and viscosity is not synonymous with turbulence. The origin of the anomalous disk stress and turbulent transport is no longer unapproachable. If $`\alpha `$ represents a stress arising from internal disk turbulence, then the source of this turbulence is magnetic (Balbus & Hawley, 1991). As discussed by Balbus & Hawley (1998) purely hydrodynamic sources of turbulence are problematic at best, and none has been observed and sustained in a fully three-dimensional simulation of a Keplerian disk. The problem is illustrated by studies of convective turbulence: when this form of turbulence is imposed, it drives inward, not outward, angular momentum transport (Stone & Balbus, 1996). This is no fluke of the peculiar properties of convection, it is what happens quite generally if incompressible turbulence is driven in a rotationally stable disk. If nonlinear instabilities were a possibility, they would be present in convective and shearing sheet simulations (Hawley, Balbus, & Winters, 1999). The shearing sheet is scale-free and any instabilities would accordingly not require extremely large Reynolds numbers to appear. (There is nothing special about the behavior of very small scale disturbances that would not also be seen at larger scales.) This is why Hawley, Balbus, & Winters (1999) easily found nonlinear planar Couette flow instabilities. Their absence from any numerical simulation of a Keplerian disk casts doubt on their existence. Trailing density waves transport angular momentum outwards of course, but their excitation and maintenance must be explained, and this is most readily done only if the disk itself is strongly self-gravitating. Even in this case, the transport properties are generally nonlocal, and not well-suited to an $`\alpha `$ formalism description (Balbus & Papaloizou, 1999). The nature of, and difficulties with, an $`\alpha `$ disk associated with direct planetary torques have recently been discussed by Goodman & Rafikov (2001) (see also Lin & Papaloizou (1996)).

In this paper we perform numerical experiments of planets in protostellar disks in which the internal turbulence and its stress arise self-consistently from the magnetorotational instability. These simulations examine how turbulence affects the formation of gaps in the disk, and, equally importantly, how the planet affects the turbulence. In this first study we will assume ideal MHD, and not address the issue of the ionization state of the disk, or the degree of magnetically coupling in the gas. These issues are far from resolved (Fromang, Terquem, & Balbus, 2002). It is of interest to examine the “simplest” case first.

The plan of this paper is as follows: in §2 we describe the numerical model used for the simulations. In §3 we present several two dimensional hydrodynamic simulations. These models serve as code tests, they provide a comparison with previous numerical work, and they are the controls for the MHD simulations which follow in §4. We present our conclusions in §5.

## 2 Numerical Model

The governing equations for the evolution of a magnetized gaseous disk in the ideal MHD limit are the equation of continuity,

$`{\displaystyle \frac{\rho }{t}}+\left(\rho 𝐯\right)`$ $`=`$ $`0,`$ (3)
the equation of motion,

$`\rho {\displaystyle \frac{𝐯}{t}}+(\rho 𝐯)𝐯`$ $`=`$ $`\left(P+{\displaystyle \frac{B^2}{8\pi }}\right)\rho \mathrm{\Phi }+\left({\displaystyle \frac{𝐁}{4\pi }}\right)𝐁,`$ (4)
and the induction equation

$`{\displaystyle \frac{𝐁}{t}}`$ $`=`$ $`\times \left(𝐯\times 𝐁\right).`$ (5)
In these equations $`\rho `$ is the disk mass density, $`𝐯`$ is the fluid velocity, $`P`$ is the pressure, $`𝐁`$ is the magnetic field vector, and $`\mathrm{\Phi }`$ is the gravitational potential. These equations are solved in a cylindrical coordinate system, $`(R,\varphi ,z)`$, in the “cylindrical disk” limit, in which vertical gravitational forces are ignored. Time and length units are fixed by setting $`GM_{}=1`$, which makes the Keplerian angular velocity, $`\mathrm{\Omega }\left[=(GM_{}/R^3)^{1/2}\right]=1`$ at $`R=1`$.

The equation of state is “locally isothermal,” with a prescribed temperature profile, $`T(R)R^1`$. This is the same equation of state used in many of the previous planet-disk simulations, e.g., Kley (1999); Bryden et al. (1999); Nelson et al. (2000). Specifically,

$$P=\rho c_s^2=\rho \left(\frac{H}{R}\right)^2\frac{GM_{}}{R},$$
(6)
where $`H/R`$ is a fixed input parameter of a given simulation. In this model, the Mach number of the flow ($`v_\varphi /c_s`$) is constant with radius. The initial disk is given a density profile $`\rho R`$ with $`\rho =1`$ at $`R=1`$. This ensures that the initial pressure gradient is zero. The equations that are explicitly integrated also contain an artificial viscosity of the form described by Stone & Norman (1992a). This is nonzero only in regions of local compression and acts only on the momentum, since the temperature at any given radius is fixed. Its principal function is to spread shock transitions over a few grid zones.

The computational grid consists of a cylindrical annulus centered on the central star, with the radial coordinate running between $`0.25<R<3.75`$. The use of an annulus avoids the singularity associated with $`R=0`$. The angle $`\varphi `$ runs over the full range from 0 to $`2\pi `$. The vertical grid size, $`L_z`$, is set equal to twice the scale height $`H`$ at the radius of the planet, $`R_p`$, i.e., for $`H/R=0.05`$, $`L_z=0.2`$. The three dimensional grid resolution is $`128\times 128\times 32`$ in $`(R,\varphi ,z)`$. We work in the ideal MHD limit of a single, infinitely conducting medium.

The gravitational potential for the central star is $`\mathrm{\Phi }=GM_{}/R`$. The planet is assumed to remain on a fixed circular orbit with a fixed orbital velocity at radius $`R_p=2`$. The planet’s orbital frequency is

$$\mathrm{\Omega }_p^2=\frac{G}{R_p^3}(M_{}+M_p),$$

and the planet’s angular location $`\varphi _p`$, is a function of time,

$$\varphi _p=\varphi _0+\mathrm{\Omega }_pt.$$

The gravitational potential of the planet is that of a softened point particle, specifically,

$$\mathrm{\Phi }_p=GM_p/(d^2+b^2)^{1/2},$$
(7)
where $`d`$ is the distance from the planet to a location within the grid, $`d^2=R^2+R_p^22RR_p\mathrm{cos}(\varphi \varphi _p)`$, and $`b`$ is the softening parameter set to be a fraction of the size of the planet’s Roche radius. In the simulations below, $`b`$ was set to be equal to one-fifth of the size of the Roche radius for a planet with $`q=0.001`$, corresponding to $`b=0.0277`$. The radial and angular gravitational force per unit volume from a planet at a point $`(R,\varphi )`$ are

$`f_R`$ $`=`$ $`\rho GM_p\left[{\displaystyle \frac{(RR_p\mathrm{cos}(\varphi \varphi _p))}{(R^2+R_p^22RR_p\mathrm{cos}(\varphi \varphi _p)+b^2)^{\frac{3}{2}}}}{\displaystyle \frac{\mathrm{cos}(\varphi \varphi _p)}{R_p^2}}\right]`$ (8)
$`f_\varphi `$ $`=`$ $`\rho GM_p\left[{\displaystyle \frac{(R_p\mathrm{sin}(\varphi \varphi _p))}{(R^2+R_p^22RR_p\mathrm{cos}(\varphi \varphi _p)+b^2)^{\frac{3}{2}}}}+{\displaystyle \frac{\mathrm{sin}(\varphi \varphi _p)}{R_p^2}}\right].`$ (9)
Since the origin of the coordinate system is non-inertial, the potential includes the indirect term (which accounts for the acceleration of the computational frame). Note that the planet is represented only as a potential; there is no mass loss from the grid at its location, and we do not attempt to model accretion onto the planet in detail.

We have run two types of simulation: inviscid hydrodynamic models in two dimensions $`(R,\varphi )`$, and three-dimensional turbulent MHD simulations. Because neither the MRI nor field amplification can operate without vertical structure, simulations of MRI-induced turbulence must be three-dimensional. However, some simplification is still required to make the problem tractable. The cylindrical disk approximation makes it practical to compute disks that have small $`H/R`$ values without the need for an unduly large number of grid zones. This cylindrical limit was previously employed by Hawley (2001) in simulations of MHD turbulence in cold Keplerian disks.

For the MHD simulations the initial magnetic field strength is set by the input parameter,

$$\beta =8\pi P/B^2,$$
(10)
which specifies the ratio of the gas pressure to the magnetic pressure. To minimize difficulties with the magnetic radial boundaries, the initial magnetic field is confined to the region between $`R=0.45`$ and $`R=3.55`$. The subsequent evolution does allow field to leave through the radial boundaries, however.

The equations were solved using time-explicit Eulerian finite differencing. The numerical algorithm is that used by the ZEUS code for hydrodynamics (Stone & Norman, 1992a) and MHD (Stone & Norman, 1992b; Hawley & Stone, 1995). The imposed boundary conditions are periodic along the azimuthal and vertical directions. The radial boundary conditions are somewhat more complex and are set in ghost zones exterior to the inner and outer radial limits of the disk. The density $`\rho `$ in the ghost zones is held fixed in time at its initial value. The azimuthal velocity $`v_\varphi `$ in the ghost zones is held constant and equal to the appropriate Keplerian value. Both the $`v_r`$ and $`v_z`$ velocity components in the ghost zones are set equal to the last physical values computed at the adjacent boundary. In the 3D MHD cases, the magnetic ghost radial boundaries are that $`B_\varphi `$ and $`B_z`$ vanish, and $`B_r`$ is set so that it is divergence-free in the ghost zones. Although straightforward to implement, these boundary conditions are at least partially reflecting, and therefore less than ideal.

The simulations produce a significant amount of both spatial and temporal data. To facilitate the analysis, two-dimensional data were averaged azimuthally, and three-dimensional quantities were averaged both vertically and azimuthally. The definition for the average of some quantity $`X`$ is

$$X=\frac{XR𝑑\varphi 𝑑z}{R𝑑\varphi 𝑑z}.$$
(11)
These averages were computed at regular time intervals for each radial grid zone, and the results used to construct space-time diagrams that characterize the overall evolution.

## 3 Hydrodynamic Tests

We begin with a series of purely hydrodynamic simulations to test the numerical model and to establish a baseline for comparison with the results from MHD simulations. Two-dimensional $`(R,\varphi )`$ hydrodynamic simulations have been carried out by a number of other researchers, e.g, Bryden et al. (1999); Lubow, Seibert, & Artymowicz (1999); Kley (1999); Nelson et al. (2000). Kley, D’Angelo, & Henning (2001) performed three-dimensional simulations and found that their results largely agree with those from two dimensional simulations so long as the scale height of the disk is comparable to or less than the Roche radius $`R_\mathrm{R}`$ of the planet, which is the tidal condition for gap formation. For $`H>R_\mathrm{R}`$ the accretion rate onto the planet can differ between two and three dimensional studies.

The disk temperature profile is set by the selected scale height ratio. The standard value, $`H/R=0.05`$, corresponds to a constant Mach number of 20 throughout the disk. The initial components of the velocity are $`v_r=0`$ and $`v_\varphi =(GM_{}/R)^{1/2}`$. Three different planet masses are used, with mass ratios $`q=2\times 10^4`$, to $`1\times 10^3`$, and $`5\times 10^3`$. Other parameters for the simulations may be found in Table 1. The three main simulations are labeled HSP, HMP, HLP for Hydrodynamic Small, Medium, and Large mass Planet respectively. The ratio of $`q`$ to the tidal limit $`q_t`$ (eqn. 1) is given for each problem. For $`H/R=0.05`$ the small mass planet has $`q<q_t`$, the medium mass is exceeds $`q_t`$, and the large mass amply satisfies the tidal criterion. For the HMPH(ot) simulation the temperature is increased sufficiently that the medium mass planet no longer meets the tidal condition for gap formation. In the HSPC(old) simulation the temperature is reduced so that the small mass planet can satisfy the tidal condition. Thus the ensemble of simulations bracket $`q_t`$ both through different planet masses and different values of $`H/R`$.

Figure 1 shows the distribution of density after 100 planetary orbits for the three simulations with the same $`H/R`$. A qualitative comparison of these plots illustrates how the width of the gap produced changes with the mass ratio: HSP fails to produce a gap, while HLP produces an extensive one. Although the small mass planet does not create a gap, the density is nevertheless reduced within a region near its orbit. In the HSP run the reduced density region is spanned by approximately 10 radial grid zones while the gaps in the HMP and HLP runs are spanned by 30 and 45 radial grid zones, respectively. Also visible are the spiral density waves that the planets generate. The amplitude of these waves increases with planetary mass. Note too the density enhancements at the trailing and leading Lagrange points, L4 and L5, particularly in the HMP run.

The density contours immediately surrounding the planet appear elongated. This is an artifact arising from the aspect ratio of the $`R`$ and $`\varphi `$ grid at the radius of the planet. The grid is uniformly spaced in both $`R`$ and $`\varphi `$ and (at the resolution used) $`Rd\varphi >dR`$ at the planet’s orbital radius. The zones are long thin rectangles, and this results in an elongated planet shape. To verify this, a test calculation was run with the number of angular zones increased by a factor of 3.5 giving $`R_pd\varphi =dR`$ at the planet’s radius. With this grid, the resulting density contours around planet are symmetric. Ideally, one would prefer to have square zones at all radii. This would be particularly important in investigating the accretion of material onto the planet. The best approach might be a nested grid around the planet (Lubow, Seibert, & Artymowicz, 1999; Ciecielag, Plewa, & Rózyczka, 2000; D’Angelo et al, 2002). In the present experiments, however, we are primarily concerned with the radial structure in the disk, not the accretion onto the planet. The grid we use is a compromise between radial resolution and a practical amount of computational time.

The top panel of Figure 2 shows the radial density profile (normalized by the initial profile) after 120 orbits for runs HSP, HMP, and HLP. In both the middle and high mass cases a gap forms. The bottom panel of Figure 2 shows the time history of the density evolution at the radial location of the planet. The slope of the line is a measure of the rate of gap clearing; not surprisingly this is a function of planet mass. The clearing rate may be characterized by $`\tau _{\mathrm{gap}}`$, defined as the amount of time required for the density at the planet’s radius to decrease by a factor of 10. The value of $`\tau _{\mathrm{gap}}`$ is obtained by fitting the density versus time at the planet’s radius with a function of the form $`\rho /\rho _o10^{t/\tau _{\mathrm{gap}}}`$, where $`t`$ is the orbital time. The value of $`\tau _{\mathrm{gap}}`$ is listed in Table 1.

The three simulations above test the effect of varying the planet mass, and two further hydrodynamical simulations test the tidal condition with varying disk temperature. In run HMPH the middle mass planet was embedded in a disk that was too hot for it to produce a gap, and in run HSPC the small mass planet was placed in a disk that was cool enough to produce a gap. Figure 3 is the density profile and time-evolution for these two runs. Comparing this data to Figure 2 one can see that cooling the disk increases the rate at which the small mass planet clears its Roche radius, reducing $`\tau _{\mathrm{gap}}`$ to 140 orbits. Heating the disk greatly increases $`\tau _{\mathrm{gap}}`$ for the middle mass planet, and eliminates the sharp gap boundaries in the density distribution.

To summarize, these hydrodynamic simulations are in agreement with the results of previous work and demonstrate the role of the tidal condition in gap formation in stress-free disks. The Roche radius determines the width of the gap formed, and the temperature in the disk determines whether that radius can be efficiently cleared by the planet’s gravity. These disk simulations will serve as the baseline comparisons for the series of disks with stress from MHD turbulence.

## 4 MHD Turbulence and Gap Formation

### 4.1 MHD turbulent disks

To establish a turbulent background flow, we begin with a planet-free MHD simulation. Previous work with cylindrical representations of MHD turbulence (Hawley, 2001; Steinacker & Papaloizou, 2002; Nelson, & Papaloizou, 2002a) provides some guidance for choosing a suitable initial field configuration. An initial vertical magnetic field produces more powerful turbulence than does an initial toroidal field, but also tends to produce a series of evacuated gaps in the disk, even without the presence of a planet. The extent to which this is physical or an artifact of the cylindrical scheme is not yet well-understood. We have therefore elected to work primarily with toroidal fields, though this is likely to underestimate the strength of the turbulence actually present in a magnetized disk.

Our fiducial simulation starts with a uniform toroidal field with $`\beta =4`$, and is run for 140 orbits at the radius at which the planet will be inserted. The relatively strong toroidal field ensures that the stress will be reasonably large. Turbulence develops promptly and the resulting stresses drive accretion. The component of greatest importance in the stress tensor is the $`R\varphi `$ component. For simplicity, we shall refer to this component as “the stress.” The Reynolds stress needs to be defined carefully, since the concept of azimuthal velocity fluctuation makes sense only if the mean azimuthal velocity is known. But the latter must, in general, be itself determined from the simulated flow. We follow Hawley (2000), and define the Reynolds stress as

$$\rho v_R\delta v_\varphi =\frac{1}{R}\left(\frac{(\rho v_R)(\rho Rv_\varphi )}{\rho }\frac{\rho v_R\rho Rv_\varphi }{\rho }\right).$$
(12)
The total turbulent stress is a sum of the Reynolds and ($`R\varphi `$) Maxwell stresses.

It is useful to have a dimensionless form for the stress. This, of course, is role of the classical $`\alpha `$ parameter, here defined as

$$\alpha \frac{1}{P_{\mathrm{initial}}}\rho v_R\delta v_\varphi \frac{B_RB_\varphi }{4\pi },$$
(13)
the stress divided by the initial pressure. The $`\alpha `$ parameter is particularly useful for comparisons with the prediction of the viscous gap formation criterion (2), and with the viscous hydrodynamic simulations, e.g., Takeuchi, Miyama, & Lin (1996), Bryden et al. (1999), Kley (1999), and Nelson et al. (2000). The function $`\alpha (R)`$ can be obtained by averaging the stress in accordance with equation (11). Figure 4 shows the $`\alpha (R)`$ profile at $`t=40`$ orbits for the baseline toroidal field simulation. The Maxwell and Reynolds stresses are also plotted separately. The Maxwell stress is clearly dominant, consistent with earlier studies of MRI-driven accretion disk turbulence (Hawley, Balbus, & Winters, 1999). The diminishing of $`\alpha `$ at the radial edges is almost certainly a consequence of the initial and boundary conditions, and the true radial dependence (assuming there is one) of $`\alpha `$ cannot be determined with certainty. One can, however, extract a reasonable global average, which is $`\alpha =0.02`$. The magnitude scale of $`\alpha `$, its fluctuation amplitudes, and its general radial variation are all similar to those found in other toroidal field cylindrical disk studies (Hawley, 2001; Steinacker & Papaloizou, 2002). Translating the globally averaged $`\alpha `$ into an effective viscosity and using it in the viscous condition (2) for gap formation, one would conclude that only the large mass planet ($`q=5\times 10^3`$) should be able to clear a gap.

The initial condition for the planet simulations corresponds to orbit 40 of the unforced toroidal field run. For comparison, the unforced run is evolved a further 100 orbits. The top panel of figure 5 shows the radial profile of density normalized by the initial profile at $`t=40`$ (dotted line) and 140 orbits (solid line). The density in the outer disk gradually decreases with time, while the density in the inner disk increases. The density at $`R=2`$, the radius at which the planet will be inserted, remains nearly constant.

The bottom panel of figure 5 shows the radial profile of the Maxwell stress normalized by the initial pressure. As the simulation evolves, the Maxwell stress in the outer disk decreases. This may be due in part to the reduction in mass and field energy in this region of the disk caused by losses through the outer boundary. As a practical matter, the time over which the simulation can be evolved is limited to $`200`$ orbits, by which point significant global evolutionary changes occur.

While the toroidal field model appears to be best suited to the present study, it is desirable to investigate other magnetic geometries, such as a vertical field. But vertical magnetic fields have proven to be so vigorously unstable (Hawley, 2001; Steinacker & Papaloizou, 2002), that large $`\alpha `$ values and substantial evolution occur over a relatively short time interval. In a cylindrical disk, an initial uniform vertical field can only be lost through the boundary, and is otherwise always present to drive the linear “channel mode” of the MRI (Hawley & Balbus, 1992). Vertical fields also tend to open up gaps in cylindrical disks even without a planet present: the Maxwell stress increases where the density decreases, and underdense regions become yet more underdense.

One way around these difficulties is to choose an initial magnetic field strength that is so small that the fastest growing unstable wavelength is unresolved. The field strength is set so that fastest growing wavelength, $`\lambda _{max}`$, is set equal to the vertical size of one grid zone, $`\mathrm{\Delta }z`$, i.e. (Balbus & Hawley, 1998),

$$B_z(R)=(15/16)^{1/2}\mathrm{\Omega }(R)(4\pi \rho )^{1/2}\frac{\mathrm{\Delta }z}{2\pi }.$$
(14)
This constraint artificially “cripples” the MRI, reducing the linear growth rate of the resolved wavelengths, and allowing long duration simulations. The linear growth rate is proportional to $`v_A/N\mathrm{\Omega }/N`$, where $`v_A`$ is the Alfvén speed, and $`N`$ is the number of grid zones of the resolved unstable mode. In the outer portion of the disk this becomes very low. Even with this crippled field, however, there is still some tendency to produce underdense gaps. The unforced simulation develops radially-alternating high and low density regions throughout the disk by orbit 140. The crippled vertical field model should be regarded not so much as a good model for a disk (it is, after all, intentionally underresolved), but as a numerical experiment that provides an alternative turbulent state for comparison with the toroidal field results.

### 4.2 Planets in turbulent disks

The parameters for the MHD simulations with a planet are listed in Table 2. The runs are labeled M for MHD, NP for no planet, or SP, MP, LP for small, medium and large mass planet, and a Z indicates a run with the initial vertical field. Except for the inclusion of magnetic fields, these planet simulations had the same parameters as the hydrodynamic tests described in §3. The value of the viscous gap parameter $`q_v`$ used in Table 2 assumes $`\alpha =0.02`$.

Figure 6 shows the density after 100 planetary orbits for runs MSP, MMP, and MLP, and figure 8 is the radial distribution of the Maxwell stress at late time for these simulations. Figure 8 illustrates that the level of turbulent stress is not independent of the planet; this has consequences. One such consequence is illustrated by MSP, a small mass planet that fails to satisfy either the tidal or the viscous gap opening condition. In MSP the density increases within the Roche radius of the planet, forming an “anti-gap.” The Maxwell stress has decreased near the planet and gas accretes into this region faster than it leaves. This density rise further decreases the local Alfvén velocity of the fluid, reducing the wavelength of the fastest growing MRI mode, further reducing the stress. Outside of the planet’s Roche radius the low mass planet produces a weak spiral wave, but otherwise has a minimal effect on the evolution of the accretion disk.

The top panel of figure 7 is a comparison of the radial density profiles for the three embedded planet toroidal MHD simulations, similar to the top panel of figure 2. The corresponding plots for each one of the planetary masses are easy to separate by the depth of the gap produced. The low mass planet simulation produces a density increase at the planet’s radius, while the other planets’ gaps are shallower than those seen in the hydro runs, and their outer edges were not as well-defined. Note that the overall slope in the radial density profile is consistent with the presence of ongoing accretion. The bottom panel of figure 7 displays the time history of the density, similar to the bottom panel of figure 2. Once again, the three different planetary masses can be separated by the rate at which the density decreases. The small mass planet shows a steady increase in density at the planet’s orbit. In the other two cases, the rate of density decrease is much less than with hydrodynamics alone; accretion hampers gap formation.

The middle mass planet (run MMP) satisfies the tidal condition but not the viscous condition. Exceeding the tidal condition ensures that the planet will exert a significant influence on the the accretion disk whether or not it is ultimately able to clear a gap. As is clear from figure 6 the planet was able to produce a noticeable density reduction around its Roche radius. The planet also creates a strong spiral wave in the disk. Most surprisingly, after about 50 planetary orbits, the magnetic stress inside $`R_p`$ decreases to very low levels (see fig. 8). Again, not only does the stress in the disk affect the planet’s ability to form a gap, but, at least in these reduced growth rate simulations, the planet can influence the turbulent stress within the disk.

The large mass planet (run MLP) easily satisfies both gap criteria, and its effect on the disk is the most extreme. A large density gap develops rapidly, although the formation time $`\tau _{\mathrm{gap}}`$ is increased over the purely hydrodynamic model. A strong global spiral wave is generated, from the inner to the outer radial boundaries. And, as with run MMP, the planet modifies the disk by lowering the Maxwell stress, here in both the inner and outer disk.

In these three runs the planet was placed into an already turbulent magnetized disk. In the next experiment, we approach the problem from another direction: we add a magnetic field to a purely hydrodynamical disk in which a gap has formed. This new simulation, labeled HMMP, is initialized by adding a $`\beta =4`$ toroidal field to the disk from simulation HMP after 100 planetary orbits. The result is a disk with just the same parameters as in simulation MMP. The ordering of the appearance of turbulence and the appearance of a planet does not matter. In MMP, a shallow gap forms even though the mass ratio $`q`$ does not exceed the viscous criterion. Adding a magnetic field to HMP has an immediate effect, and accretion begins to fill in the much deeper gap produced without stress.

In figure 9, the top panel shows the radial density profile at the moment the magnetic field is added to the disk (dotted line), and after 100 orbits of MHD evolution (solid line). The bottom panel depicts the time history of the density at the planet’s radius. The initial hydrodynamic phase exhibits an exponential decrease in density as the gap deepens. However, once the toroidal magnetic field is included and turbulence develops, the density at the planet’s radius begins to increase, returning to about the same level that it had at the end of run MMP. At the end of the simulation the gap has become shallower, wider, and much less well-defined than in the hydrodynamic simulation. The turbulent stresses drive accretion across the gap, increasing the mass in the inner disk at the expense of the outer disk. Apparently a planet that satisfies the tidal criterion, but not the viscous criterion, can form only a partial gap.

We have also examined the effect of varying $`H/R`$ by altering the temperature of the underlying disk. To keep the average $`\alpha `$ value comparable to the previous simulations, the initial magnetic field strength must also be adjusted. The ratio of the fastest growing wavelength to the disk radius, $`\lambda _{max}/R`$, is held constant by varying $`\beta `$, since

$$\frac{\lambda _{max}}{R}\frac{1}{\sqrt{\beta }}\left(\frac{H}{R}\right).$$
(15)
Run MMPH is a hotter disk ($`H/R=0.09`$) with the medium mass planet. This planet does not satisfy either the tidal or viscous criterion for gap formation. In fact, an anti-gap, or over-density, is produced, just as in the MSP simulation.

Run MSPC is a cooler disk with $`H/R=0.02`$, $`\beta =0.64`$, and a small mass planet. It satisfies the tidal criterion, but not the viscous criterion. During the first 60 orbits of evolution, the density drops inside the planet’s Roche radius; after that time the rate of decrease declines nearly to zero. The density in the gap levels off at around 8% of its original value, similar to the results of the MMP simulation. When the tidal criterion is satisfied, the density within the Roche radius is reduced, but if the viscous condition is not met, the gap is not completely emptied.

Initially, in run MSPC there is a measurable Maxwell stress throughout the disk, but once the planet is inserted, the magnetic stress in the inner disk is reduced. A comparison between this simulation and MSP is instructive. When the planet was not able to produce a density gap in MSP, the Maxwell stress was largely unaffected. In the cooler disk, however, the small mass planet did produce a gap, and the turbulence in the inner disk was diminished. This suggests that the threshold for gap production is closely related to the threshold for affecting the disk stress. Evidently, both processes mark the onset of nonlinear disk-planet interactions.

In run MZMP, a medium mass planet in a vertical field, the disk was evolved for 100 orbits, allowing a gap to form. In contrast to run MMP, the Maxwell stress was sustained inside of the planet’s orbit, but it is significantly reduced outside of the planet’s orbit. Of course even in run MZNP which does not have a planet, it is clear that at large radii in the MHD turbulence is rather weak, so this effect cannot be attributed solely to the presence of the planet. The continued presence of Maxwell stresses in the inner portion of the disk, however, is a significant difference from what is seen in the toroidal field simulation.

Comparing the density evolution in run HMP with MZMP we find remarkably good agreement between the two simulations. Even though the turbulent stresses are initially present in the MHD simulation, the resulting density profile and time history are essentially those of the hydro simulation. The reason for this is clear. In the MHD simulation at late times, the Maxwell stress in the outer disk is drastically reduced, effectively ending inward accretion. With no fluid moving from the outer disk to the inner, nothing flows across the gap. Consequently, the rate of gap formation and radial profile of the gap are the same as in the hydro simulation.

### 4.3 Effect of planet on stress

The simulations have provided evidence that not only does the stress affect the formation of a gap around the planet, but, under some circumstances, the planet can directly influence the turbulent stress. In run MSP, the stress was reduced only within the planet’s Roche radius, while in the larger mass planet runs (MMP and MLP), throughout most of the inner disk the Maxwell stress was significantly reduced at late times. This occurs more rapidly with the large mass planet.

The question naturally arises whether the presence of the planet directly causes the reduction of the MHD stress, or whether changes in the disk and magnetic field due to the overall evolution have altered the MRI stability properties of the disk. To distinguish between these possibilities we return to the end of run MMP and remove the planet. The disk, magnetic field, and gap remain, but the spiral wave is quickly lost once the planet is gone.

Figure 10 shows the before and after radial profiles of density (top panel) and Maxwell stress (bottom panel). The solid lines correspond to the point in time when the planet was removed, and the dashed lines are from 20 orbits later. The bottom panel dramatically shows how the Maxwell stress in the inner disk recovers after the removal of the planet, increasing by a factor of 10 after 20 orbits. The density distribution, on the other hand, changes only slightly over this same time. The most noticeable change in the radial density profile is that the gap has become shallower and its edges less defined. This experiment suggests that interactions between the spiral density wave and the MHD turbulence affects the level of the Maxwell stress.

The planet’s influence on the disk as a whole is clearly manifest as a large spiral wave. Figure 11 displays vertically-averaged Maxwell stress images from planetless simulation and the two small-planet MHD simulations, MSP and MSPC. Without a planet, the magnetic stress is spread evenly throughout the disk. In the MSP run, the planet is not able to form a gap. The planet produces a spiral wave that is strong only near the planet. When the disk is cooled in the MSPC simulation, the small mass planet is able to form a gap. As can be seen from the bottom image, the planet now has a much more extended influence on the magnetic stress. The inference is that if the planet is massive enough to create a gap, it is massive enough to produce a spiral density wave that interacts with, and lowers, the turbulent Maxwell stress in the disk interior.

It should be noted that this effect may well be sensitive to the choice of a toroidal field. The Maxwell stress in the crippled vertical field run MZMP was not substantially reduced by the planet, even though it did form a gap. This point requires further study.

## 5 Conclusions

In this paper we have examined the formation of gaps by planets in protostellar disks in the presence of Maxwell and Reynolds stresses produced by MHD turbulence. The aims of this first study were to test the tidal and viscous conditions for gap formation, and to observe the effects of the planet on an MHD turbulent disk.

We began with a series of inviscid hydrodynamic simulations as controls. Under these conditions gap formation is determined by the tidal criterion: a planet will clear a gap when the local Roche radius is greater than the pressure scale height. Our hydrodynamic simulations agreed with this condition and with the many simulations of such disks already in the literature.

To establish a turbulent disk, we ran an initial toroidal field configuration for 40 planetary orbits prior to the insertion of the planet itself. The resulting disk had significant Maxwell and Reynolds stresses that vary with radius, but which gave an average $`\alpha `$ value of 0.02. A planet was then added, and the results compared with the tidal and with the viscous criteria for gap formation.

When the planet was insufficiently massive, or the disk too hot to satisfy the tidal condition, a gap did not form. In fact, the density around the planet increases. The MHD turbulence drives an overall accretion flow from the outer to the inner disk, and because the planet is able to reduce the turbulent stress within its Roche radius, there is a local buildup in density.

In several simulations, the planet satisfied the tidal criterion, but not the viscous criterion. In this case, the density within the Roche radius is still significantly reduced, but not to zero. The resulting gaps are more shallow than those seen in the pure hydro runs, and the radial density profiles at outer edges of the gap are wider.

When a planet satisfies both the viscous and tidal criteria, a gap formed, but at a reduced rate compared with the hydrodynamic runs. Accretion across the gap continues, and the edges of the gap are somewhat smeared out. The accretion rate is not large enough to prevent the density within the gap from being driven to a very low level.

In a recent paper, Nelson, & Papaloizou (2002b) present a simulation of a planet placed in an MHD turbulent disk. Their simulation considered a $`q=5\times 10^3`$ mass planet (equivalent to our high mass planet) on a cylindrical grid with a disk model that has the same radial power law dependence on $`\rho `$ and $`P`$ in the main part of the disk as the model as described in §2. They use $`H/R=0.1`$ (hotter than our models), and an initial sinusoidally-varying vertical field that produces MHD turbulence with an averaged $`\alpha `$ value of $`0.005`$ compared with our models that use a stronger toroidal field with $`\alpha 0.02`$. Nelson & Papaloizou ran a single simulation with higher resolution than our models. The focused on the interaction of the planet with a disk and a pre-existing gap, whereas our aim is to consider the gap formation process itself in the presence of MHD turbulence.

Nelson, & Papaloizou (2002b) found, as we did, that spiral waves produced by the planet reduced the stresses provided by the MHD turbulence; they observed this effect mainly near the gap where the waves were strong. They do not observe evidence for accretion across the gap, whereas some of our simulations show this. This suggests that whether or not there is accretion across the gap depends on such parameters as $`q`$, $`\alpha `$, and $`H/R`$.

Some results do not differ dramatically from those seen in viscous hydrodynamic studies. For example, Kley (1999) found that the tidal condition had a greater effect on the width of the gap and the strength of the spiral waves than variations in viscosity. By determining the rate at which gas flows into the gap, the viscosity affects mainly the accretion rate onto the planet. Nelson, & Papaloizou (2002b) compared their model directly with a viscous model and found a number of differences in detail. The gap in the turbulent disk model was deeper than the viscous model with equivalent $`\alpha `$ value.

Clearly, matters are much more interesting than the viscous models would predict. The presence of the planet itself directly influences the turbulence levels and stress within the disk. This could not be seen with a standard $`\alpha `$ viscosity simulation, for in such a model the viscosity can influence, but not be influenced by, the planet. We found that the small mass planet reduced the Maxwell stress in its immediate vicinity, allowing matter to accumulate here. More significantly, the medium and large mass planets were able to reduce the turbulence levels in the disk region interior to their orbits. Tests showed that this was not because the MRI was stabilized by changes within the disk structure. Rather, it appears to be due to a nonlinear interaction between the spiral wave induced by the planet and the turbulence. This wave has to exceed a certain strength to have this effect. In the main part of the disk the rule is when the planet is massive enough to open a gap it can also reduce the turbulent stress. It should be noted that this effect was observed only in the disks where the initial field was toroidal, and not in the simulation with a weak vertical field.

We acknowledge support under NSF grant AST-0070979, and NASA grants NAG5-10655 and NAG5-9266. Some of the simulations describe here were carried out on computational platforms supported by the National Center for Supercomputing Applications and the San Diego Supercomputer Center.


format: text

text:
# Structure Function Scaling in Compressible Super–Alfvénic MHD Turbulence

## Abstract

Supersonic turbulent flows of magnetized gas are believed to play an important role in the dynamics of star–forming clouds in galaxies. Understanding statistical properties of such flows is crucial for developing a theory of star formation. In this letter we propose a unified approach for obtaining the velocity scaling in compressible and super–Alfvénic turbulence, valid for arbitrary sonic Mach number, $`M_\mathrm{S}`$. We demonstrate with numerical simulations that the scaling can be described with the She–Lévêque formalism, where only one parameter, interpreted as the Hausdorff dimension of the most intense dissipative structures, needs to be varied as a function of $`M_\mathrm{S}`$. Our results thus provide a method for obtaining the velocity scaling in interstellar clouds once their Mach numbers have been inferred from observations.

1. Introduction. Star formation occurs in our galaxy primarily inside cold interstellar clouds. Observations indicate that gas motion inside these clouds is highly supersonic in the broad range of scales $`0.11000`$ pc. At these scales the clouds lack any characteristic structure and their scaling properties are to some extent universal. This universality is the signature of the supersonic and super–Alfvénic magneto–hydrodynamic (MHD) turbulence that governs the dynamics of the clouds Padoan et al. (1999); Ballesteros-Paredes et al. (1999); Padoan and Nordlund (1999); Klessen et al. (2000); Ostriker et al. (2001); Padoan et al. (2002a, b).

The properties of turbulent MHD flows are controlled by the rms sonic Mach number, $`M_\mathrm{S}`$, (ratio of rms flow velocity and sound velocity) and the rms Alfvénic Mach number, $`M_\mathrm{A}`$, (ratio of rms flow velocity and Alfvén velocity). In interstellar clouds these numbers are large, $`M_\mathrm{S}130`$, $`M_\mathrm{A}>1`$. The corresponding supersonic motion produces the observed strong density enhancements Padoan et al. (2001), some of which are Jeans–unstable and undergo gravitational collapse. This initial ‘turbulent’ stage of density fragmentation should be described statistically, in terms of velocity and density probability distribution functions or correlators; its properties are inherently related to the statistical properties of compressible MHD turbulence Padoan and Nordlund (2002).

Early studies suggesting the importance of astrophysical supersonic turbulence were undertaken by von Weizsäcker von Weizsäcker (1951) and later by Larson Larson (1981), who discovered that the spectrum of the velocity field in molecular clouds was steeper than the spectrum predicted by the standard Kolmogorov model. The discrepancy was attributed to the strong intermittency of supersonic turbulence. However, the analytical understanding of such turbulence remained rather poor.

Advances in computational and observational techniques have recently allowed measurements of power spectra and higher order correlators of velocity and density fields, reviving the interest in supersonic turbulence. Modern observations generate surveys covering a range of scales comparable or even exceeding the resolution of numerical simulations. The interstellar medium can now be viewed as a laboratory for investigating supersonic motion unachievable with Earth–based experiments.

The turbulent field can be characterized by the scaling of its structure functions \[see below\]. The intermittency is defined as the departure of this scaling from the Kolmogorov one. So far, only the limiting cases of $`M_\mathrm{S}1`$ and $`M_\mathrm{S}=0`$ have been investigated. An analytical description of supersonic turbulence, $`M_\mathrm{S}1`$, has been proposed in Boldyrev (2002) and confirmed numerically in Boldyrev et al. (2002a, b); Kritsuk and Norman (2003). It is based on the observation that the most intense dissipative structures in supersonic flows are two–dimensional shocks. This approach makes use of the model by She and Lévêque She and Lévêque (1994) for incompressible turbulence, $`M_\mathrm{S}=0`$, where the dissipative structures are one–dimensional vortex filaments.

In this letter we report the results of a series of numerical MHD simulations where we traced the change of the velocity structure functions over the broad range of Mach numbers $`0.3M_\mathrm{S}10`$.

We discovered that their scalings can be described by the general model where only one parameter, the dimension of the most intense dissipative structures, needs to be changed with $`M_\mathrm{S}`$. The scalings obtained in She and Lévêque (1994); Boldyrev (2002) naturally appear as the limiting cases of this unifying model. Our results suggest that the velocity scaling within interstellar clouds can be inferred from measurements of their Mach number. Observationally, this is a much easier task than measuring the structure function directly. We present the main results here; the detailed discussion will appear elsewhere.

2. Velocity Structure Function Scaling. The velocity structure functions of order $`p`$ are defined as:

$$S_p(\mathrm{})=|u(x+\mathrm{})u(x)|^p\mathrm{}^{\zeta (p)}$$
(1)
where the velocity component $`u`$ is parallel (longitudinal structure function) or perpendicular (transversal structure function) to the vector $`\mathrm{}`$ and the spatial average is over all values of the position $`x`$.

The scaling of velocity structure functions in incompressible turbulence is best described by the She–Lévêque formula She and Lévêque (1994); Dubrulle (1994):

$$\zeta (p)/\zeta (3)=p/9+2\left[1\left(2/3\right)^{p/3}\right],$$
(2)
This is a specific case of the more general formula

$$\zeta (p)/\zeta (3)=\gamma p+C(1\beta ^p)$$
(3)
with the normalization condition,

$$C=(13\gamma )/(1\beta ^3).$$
(4)
Two basic assumptions are used to derive (3), as discussed in She et al. (2001). The first assumption is the existence of a symmetric relation between the intensity of fluctuations $`F_p(\mathrm{})`$ and $`F_{p+1}(\mathrm{})`$ with respect to a translation in $`p`$:

$$F_{p+1}(\mathrm{})=A_pF_p(\mathrm{})^\beta F_{(\mathrm{})}(\mathrm{})^{1\beta },$$
(5)
where the $`p`$th order intensity of fluctuations is defined as $`F_p(\mathrm{})=S_{p+1}(\mathrm{})/S_p(\mathrm{})`$, and $`A_p`$ are constants independent of $`\mathrm{}`$ (generally found to be independent of $`p`$ as well). The functions $`F_p(\mathrm{})`$ are associated with higher intensity fluctuations with increasing values of $`p`$. The parameter $`\beta `$ is a measure of intermittency. In Kolmogorov’s 1941 model Kolmogorov (1941), $`\beta 1`$ and $`\zeta (p)=p/3`$ This is the limit of no intermittency. Lower $`\beta `$ corresponds to higher degree of intermittency, or to the increasing role played by strong persistent fluctuations (structures).

The second assumption is

$$F_{\mathrm{}}S_3^\gamma .$$
(6)
The value of the parameter $`\gamma `$ is related to very high order moments and is therefore difficult to constrain with experiments. Kolmogorov’s turbulence corresponds to $`\gamma =1/3`$.

The parameter $`C`$ has been interpreted as the Hausdorff codimension of the support of the most singular dissipative structures Dubrulle (1994). In incompressible turbulence the most dissipative structures are along vortex filaments with Hausdorff dimension $`D=1`$ and so $`C=2`$.

Boldyrev Boldyrev (2002) has proposed to extend the scaling (3) to supersonic turbulence, with the assumption that the Hausdorff dimension of the support of the most singular dissipative structures is $`D=2`$ ($`C=1`$), because dissipation of supersonic turbulence occurs mainly in sheet–like shocks. He used $`\gamma =1/9`$ as in incompressible turbulence, which implies $`\beta ^3=1/3`$, that is a higher degree of intermittency than in incompressible turbulence. With these parameters one obtains the velocity scaling proposed for supersonic turbulence in Boldyrev (2002):

$$\zeta (p)/\zeta (3)=p/9+1\left(1/3\right)^{p/3}.$$
(7)
This velocity scaling has been found to provide a very accurate prediction for numerical simulations of highly supersonic and super–Alfvénic turbulence Boldyrev et al. (2002a, b). The same scaling has also been proposed for incompressible MHD turbulence Müller and Biskamp (2000), where dissipation occurs mainly in two dimensional current sheets.

In the present work we obtain the structure function scaling in simulations of compressible super–Alfvénic MHD turbulence as a function of the rms sonic Mach number of the flow, $`M_\mathrm{S}`$, for constant value of the rms Alfvénic Mach number, $`M_\mathrm{A}`$$`1`$ (the limiting case of strong magnetic field was investigated analytically and numerically in Lithwick and Goldreich (2001); Cho et al. (2002); Cho and Lazarian (2002)). We use an isothermal equation of state and vary the value of $`M_\mathrm{S}`$ in different experiments by varying the thermal energy. As a result, the initial value of $`M_\mathrm{A}`$=10 remains unchanged from run to run. In our numerical method the total magnetic flux is conserved. However, the magnetic energy is amplified.

The initial value of the ratio of average magnetic and dynamic pressures is $`P_\mathrm{m}_{\mathrm{in}}/P_\mathrm{d}_{\mathrm{in}}=0.005`$ for all runs.

The value of the same ratio, averaged over the time interval for which the structure functions are computed, is given in Table 1. Because the average dynamic pressure is always in excess of the average magnetic pressure, $`P_\mathrm{m}/P_\mathrm{d}1`$, the dissipation is expected to occur primarily in vortex filaments (roughly corresponding to magnetic flux tubes) in subsonic runs and in shocks in supersonic runs. In the next section we argue that when we increase the value of $`M_\mathrm{S}`$, the Hausdorff dimension of the support of the most dissipative structures grows continuously from $`D1`$ to $`D2`$.

3. Results. We have solved the three dimensional compressible MHD equations, for an isothermal gas with initially uniform density and magnetic fields, in a staggered mesh with 2503 computational cells and periodic boundary conditions. Turbulence is set up as an initial large scale random and solenoidal velocity field (generated in Fourier space with power only in the range of wavenumber $`1K2`$) and maintained with an external large scale random and solenoidal force, correlated at the largest scale turn–over time. Details about the numerical method are given in Padoan and Nordlund (1999).

We have run six experiments with different values of $`M_\mathrm{S}`$, in the range $`0.4M_\mathrm{S}9.5`$ (Table 1), for approximately 10 dynamical times. The transverse structure functions of velocity are computed for the last four dynamical times of each experiment, up to the order $`p=10`$. Their power law slope, $`\zeta (p)`$, is obtained from a least square fit in the range $`2\mathrm{}96`$. The standard deviation of the slope is used as an estimate of the uncertainty of $`\zeta (p)`$.

Figure 1 shows the structure functions plotted versus the third order one for the $`M_\mathrm{S}`$=9.5 experiment. The structure functions are compensated by their least square fit slope, $`\zeta (p)`$. We study the structure functions relative to the third order to exploit the concept of extended self–similarity Benzi et al. (1993); Dubrulle (1994), that is the property of the relative scaling to extend up to a range of scales affected by dissipation (numerical dissipation in this case).

The exponents $`\zeta (p)`$, normalized to the third order one, are plotted in Figure 2, divided by their value predicted in Boldyrev (2002) and given by (7). The $`M_\mathrm{S}`$=0.4 and $`M_\mathrm{S}`$=0.8 models follow the She–Lévêque scaling of incompressible turbulence. As the value of $`M_\mathrm{S}`$ increases, the velocity scaling becomes gradually more intermittent. At $`M_\mathrm{S}`$=9.5 the scaling is very close to the analytical prediction for supersonic turbulence Boldyrev (2002), as previously shown Boldyrev et al. (2002a, b).

The validity of the assumption of hierarchical symmetry expressed by (5), can be tested with a log–log plot of $`F_{p+1}/F_2`$ versus $`F_p/F_1`$,

obtained by varying the argument $`\mathrm{}`$ in (5). If the plot is a straight line, the parameters $`A_p`$ are constants independent of $`\mathrm{}`$ and the assumed hierarchical symmetry expressed by (5) is satisfied. This is called the “$`\beta `$ test” in She et al. (2001). Our numerical simulations pass the $`\beta `$ test for all values of $`M_\mathrm{S}`$. In Figure 3 we show the two extreme cases, $`M_\mathrm{S}`$=0.4 and $`M_\mathrm{S}`$=9.5. We have used values of $`F_p(\mathrm{})`$ for all values of $`\mathrm{}`$ used in the least square fit estimate of the exponents $`\zeta (p)`$. In all the runs we find that the constants $`A_p`$ in (5) are independent not only of $`\mathrm{}`$, but also of $`p`$, as in previous numerical and laboratory studies Lévêque and She (1997); She et al. (2001); Chavarria et al. (1995). The plot in Figure 3 is in fact obtained varying both $`\mathrm{}`$ and $`p`$.

The slope of the plot in Figure 3 provides an estimate of $`\beta `$. For the models with $`M_\mathrm{S}`$ from 0.4 to 9.5 we obtain $`\beta =0.82`$, 0.81, 0.75, 0.74, 0.77 and 0.77. We instead find no evidence of a systematic variation in the parameter $`\gamma `$ with $`M_\mathrm{S}`$ and our data is consistent with $`\gamma =1/9`$.

The $`\beta `$ test gives support to the basic assumption of the She–Lévêque formalism, but does not necessarily yield the value of $`\beta `$ providing the best fit to the structure function scaling $`\zeta (p)/\zeta (3)`$. Values of $`\beta `$ computed from a best fit to $`\zeta (p)/\zeta (3)`$ from our numerical experiments are slightly different from the $`\beta `$ test values (see Table 1). We use these best–fit $`\beta `$ values to compute the parameter $`D`$.

¿From the definition of the parameter $`C`$ in (4), the values of $`\beta `$ estimated above and $`\gamma =1/9`$ we can compute $`D=3C`$ as a function of $`M_\mathrm{S}`$. The values of $`\beta `$ and $`D`$ are given in Table 1. The value of $`D`$ grows with $`M_\mathrm{S}`$, from $`D1.0`$ in the subsonic runs to $`D2.0`$ in the most supersonic runs.

4. Conclusions. We propose a unifying description of the velocity structure function scaling of compressible and super–Alfvénic MHD turbulence for an arbitrary value of the rms sonic Mach number of the flow, $`M_\mathrm{S}`$. Only one model parameter, interpreted as the Hausdorff dimension of the support of the most dissipative structures in the She–Lévêque formalism, needs to be varied as a function of $`M_\mathrm{S}`$. Besides its general importance, our result implies that the velocity structure function scaling in interstellar clouds can be predicted because the value of $`M_\mathrm{S}`$ in the clouds is easily estimated from observations. This may help improve our understanding of the process of cloud fragmentation that leads to the birth of stars.

###### Acknowledgements.

RJ was supported by NSF grant AST-0206031. ÅN was supported by a grant from the Danish Natural Science Research Council, and in part by the Danish National Research Foundation, through its establishment of the Theoretical Astrophysics Center. SB was supported by the NSF Center for magnetic self–organization in astrophysical and laboratory plasmas at the University of Chicago.


text:
# Crossover from weak localization to weak antilocalization in a disordered microbridge

## Abstract

We calculate the weak localization correction in the double crossover to broken time-reversal and spin-rotational symmetry for a disordered microbridge or a short disordered wire using a scattering-matrix approach. Whereas the correction has universal limiting values in the three basic symmetry classes, the functional form of the magnetoconductance is affected by eventual non-homogeneities in the microbridge.

Interference of time-reversed paths causes a small negative quantum correction to the conductance of a disordered metal termed the weak localization.Anderson79 ; Gorkov79 ; bergmann84 ; LeeRamakrishnan  This correction is suppressed by a time-reversal symmetry breaking magnetic field, whereas in the presence of strong spin-orbit scattering, the sign of the correction is reversed.hikami80  In that case, the interference correction is known as weak antilocalization.

In a wire geometry at zero temperature, the weak localization correction takes a particularly simple and universal form,BeenakkerReview

$$\delta G=\frac{2e^2(\beta 2)}{3\beta h},$$
(1)
where the symmetry parameter $`\beta `$ denotes the appropriate symmetry class: In the presence of an applied magnetic field, $`\beta =2`$, and without a magnetic field, $`\beta =4`$ or 1 with or without strong spin-orbit scattering, respectively. Equation (1) was obtained using random-matrix theory,mello88 ; mello91 ; macedo92  and diagrammatic perturbation theory,LeeRamakrishnan ; mello91  and is valid if the length $`L`$ of the wire is much smaller than the localization length $`\xi `$ and the dephasing length $`L_\varphi `$, but much larger than the mean free path $`l`$. The validity of Eq. (1) extends to the case when sample parameters are non-homogeneous, e.g., for wires of varying cross section, mean free path, or electron density.

For wires with weak spin-orbit scattering, a crossover between weak localization and weak antilocalization takes place when the spin-orbit scattering length $`l_{\mathrm{so}}`$ becomes comparable to $`L`$ or $`L_\varphi `$ (whichever is smaller). Experimentally, this crossover regime has been well studied in wires with length $`LL_\varphi `$.kurdak92 ; Moon ; Gougam  In this regime, weak (anti)localization takes the form of a small correction to the conductivity of the wire, rather than of a correction to the conductance. Theoretically, the weak localization to weak antilocalization crossover in the regime $`LL_\varphi `$ has been considered in Refs. aa, ; santhanam84b, ; santhanam87b,  using diagrammatic perturbation theory. The opposite regime $`LL_\varphi `$, where the universal correction (1) to the conductance $`G`$ can be observed, would be relevant for relatively short high-purity metal wires,Pothier  or disordered microbridges.

The goal of this paper is threefold: (i) to generalize the random-matrix methods for quantum wires to the crossover between weak localization and weak antilocalization, thus extending the equivalence of the two methods to the interpolation between the three symmetry classes, (ii) to find an explicit expression for $`\delta G`$ for $`LL_\varphi `$, and (iii) to extend the theory for the crossover regime to the case of non-homogeneous wires, for which the electron density, impurity concentration, or cross section varies along the sample. In this case, both the crossover scale and the functional form of $`\delta G`$ in the crossover are affected by non-homogeneities. The fact that the crossover scale, characterized by the spin-orbit length $`l_{\mathrm{so}}`$ and the magnetic length $`l_H`$, is non-universal is well known, both for homogeneous and for non-homogeneous microbridges.BeenakkerVanHouten  Our finding that the functional form of the crossover is affected by the non-homogeneity is markedly different from crossovers between the three basic symmetry classes in quantum dots, where the functional forms are universal and given by random-matrix theory.BeenakkerReview  For homogeneous wires, $`\delta G`$ is a universal function of $`L/l_{\mathrm{so}}`$ and $`L/l_H`$.

The main assumption underlying our calculations is that the wire width $`WL`$, i.e.  quasi one-dimensionality. We also assume that the wire is well in the diffusive regime, $`lL,l_{\mathrm{so}},l_H\xi `$, where $`l`$ is the elastic mean free path, and, for a non-homogeneous microbridge, that the number of propagating channels at the Fermi level $`N`$ has only one minimum along the wire (excluding the possibility of a “cavity”). We first discuss our calculations for homogeneous wires; the case of non-homogeneous samples is discussed at the end of this paper.

Starting point of our calculation is a random-matrix model similar to that used by Dorokhov.Dorokhov88  A disordered wire with $`N`$ propagating channels at the Fermi level is modeled by $`N`$ one-dimensional channels and periodically inserted scatterers that scatter within and between the channels. The electronic wavefunction is represented by a $`2N`$-component vector of spinors. The $`2N`$ components of the wavefunction refer to the transverse channel and to the left/right mover index. Linearizing the kinetic energy in each of the channels, the Hamiltonian $`H`$ takes the form of a differential operator with respect to the coordinate $`x`$ along the wire and a $`2N`$-dimensional quaternion matrix with respect to the channel and left/right mover indices and spinor degree of freedom,

$$H=i\sigma _0\tau _3𝟙_{}\frac{}{𝕩}+\underset{𝕛}{}𝕍_𝕛\delta (𝕩𝕛𝕒),$$
(2)
with $`\sigma _0`$ the $`2\times 2`$ unit matrix for the spinor degree of freedom, $`\tau _3`$ the Pauli matrix in left-mover/right-mover grading, $`𝟙_{}`$ the $`N\times N`$ unit matrix in the channel grading, $`V_j`$ a Hermitian $`2N\times 2N`$ quaternion matrix representing the $`j`$th scatterer along the wire, and $`a`$ the distance between scatterers. A quaternion is a $`2\times 2`$ matrix acting in the spinor grading with special rules for transposition and complex conjugation Mehta : The “dual” $`X^\mathrm{R}`$ of a quaternion matrix is $`X^\mathrm{R}=\sigma _2X^\mathrm{T}\sigma _2`$; the quaternion complex conjugate is defined as $`X^{}=(X^{})^\mathrm{R}`$. We have chosen units such that the Fermi velocity is one. A model similar to Eq. (2) has been used in Ref. brouwer02a,  to study weak localization in unconventional superconducting wires.

The ensemble-averaged conductance $`G`$ of the wire is given by the Landauer formula,

$$G=\frac{e^2}{h}g,g=\text{tr}(1r^{}r),$$
(3)
where $`r`$ is the $`N\times N`$ quaternion reflection matrix of the wire. To calculate $`r`$, we start from a wire of zero length and add slices of length $`a`$ at the wire’s ends. The scattering matrix of the $`j`$th scatterer is

$$S_j=\left(\begin{array}{cc}t_j& r_j^{}\\ r_j& t_j^{}\end{array}\right)=\frac{2iV_j}{2i+V_j}.$$
(4)
Hence, if a scatterer is added at the lead end of the wire, the new reflection matrix of the wire is calculated according to the composition rule

$$rr_j+t_j^{}r(1r_j^{}r)^1t_j.$$
(5)
(A similar composition rule, involving both transmission and reflection matrices of the disordered wire, applies if a scatterer is added at the far end of the wire.BeenakkerReview )

In left-mover/right-mover grading, the potential $`V_j`$ is parameterized as

$$V=\left(\begin{array}{cc}v_{LL}& v_{LR}\\ v_{RL}& v_{RR}\end{array}\right),$$
(6)
where $`v_{LL}`$, $`v_{LR}`$, $`v_{RL}`$, and $`v_{RR}`$ are $`N\times N`$ quaternion matrices,

$`v_{LL}(\alpha _f,\eta _f)=v_{RR}^{}(\alpha _f,\eta _f)`$
$`=`$ $`\sqrt{{\displaystyle \frac{a}{l_fN}}}\left[(u_f^0+\eta _fx_f)\sigma _0+i\alpha _f{\displaystyle \underset{\mu =1}{\overset{3}{}}}u_f^\mu \sigma _\mu \right],`$
$`v_{LR}(\alpha _b,\eta _b)=v_{RL}^{}(\alpha _b,\eta _b)`$
$`=`$ $`\sqrt{{\displaystyle \frac{a}{l(N+1)}}}\left[(u_b^0+\eta _bx_b)\sigma _0+i\alpha _b{\displaystyle \underset{\mu =1}{\overset{3}{}}}u_b^\mu \sigma _\mu \right].`$
In Eq. (7), $`u_f^0`$ and $`x_f`$ are random Hermitian $`N\times N`$ matrices, $`u_f^\mu `$, $`\mu =1,2,3`$, is a random anti-Hermitian matrix, $`u_b^0`$ is a random symmetric matrix, and $`u_b^\mu `$, $`\mu =1,2,3`$ and $`x_b`$ are random antisymmetric matrices. All of these random matrices have independent and Gaussian distributions with zero mean and unit variance. (Variances are specified for the off-diagonal elements; diagonal elements have double variance for symmetric matrices and are zero for antisymmetric matrices.) The parameters $`\alpha _b`$ and $`\alpha _f`$ describe the strength of the breaking of spin-rotational symmetry; The parameters $`\eta _b`$ and $`\eta _f`$ describe the strength of the breaking of time-reversal symmetry. Finally, $`l_f`$ is the elastic mean free path for forward scattering and $`l`$ is the transport mean free path.

To find the conductance of the wire we calculate the change of $`g`$ if one scatterer is added to the wire. To this end, we expand the scattering matrix $`S_j`$ of Eq. (4) in powers of $`V_j`$, use the composition rule (5), and calculate the Gaussian average over the potential $`V_j`$. In the limit $`al`$ of weak disorder we thus find

$$2Nl\frac{}{L}g=g^2h_0+3h_1.$$
(8)
We abbreviated

$`h_0`$ $`=`$ $`\text{tr}(1r^{}r)(1r^{}r^\mathrm{R}),`$ (9a)
$`h_1`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mu =1}{\overset{3}{}}}\text{tr}(1r^{}r)\sigma _\mu (1r^{}r^\mathrm{R})\sigma _\mu ,`$ (9b)
and omitted terms that vanish in the diffusive regime $`lL,l_{\mathrm{so}},l_HNl`$. The subscripts $`0`$ and $`1`$ refer to singlet and triplet contributions, respectively.

To leading order in $`N`$, Eq. (8) can be solved without the interference corrections $`h_0`$ and $`h_1`$, with the result

$$g=\frac{2Nl}{L}+O(1),$$
(10)
corresponding to the Drude law for the conductance. The $`O(1)`$ correction in Eq. (10) gives the weak localization correction $`\delta g`$, which we now compute.

To find the weak localization correction, we need to calculate $`h_0`$ and $`h_1`$. Proceeding as before, we find that the $`L`$-dependence of $`h_m`$, $`m=0,1`$ is governed by the evolution equation

$$2Nl\frac{h_m}{L}=2\left(\frac{2Nl}{L}+k_m\right)h_m+\frac{8N^2l^2}{L^2},m=0,1,$$
(11)
where we abbreviated

$`k_0=\text{tr}(1r^{}r),k_1={\displaystyle \frac{1}{3}}{\displaystyle \underset{\mu =1}{\overset{3}{}}}\text{tr}(1r^{}\sigma _\mu r\sigma _\mu ).`$ (12)
Evolution equations for $`k_0`$ and $`k_1`$ are obtained similarly and read

$`2Nl{\displaystyle \frac{k_0}{L}}`$ $`=`$ $`\left({\displaystyle \frac{2Nl}{l_H}}\right)^2k_0^2,`$ (13a)
$`2Nl{\displaystyle \frac{k_1}{L}}`$ $`=`$ $`\left({\displaystyle \frac{2Nl}{l_H^{}}}\right)^2k_1^2,`$ (13b)
where the length scales $`l_H`$ and $`l_H^{}`$ are defined in terms of the parameters of the random-matrix model (7),

$`l_H^2`$ $`=`$ $`2(l^2\eta _b^2+l^1l_f^1\eta _f^2),`$ (14a)
$`l_{\mathrm{so}}^2`$ $`=`$ $`6(l^2\alpha _b^2+l^1l_f^1\alpha _f^2),`$ (14b)
$`(l_H^{})^2`$ $`=`$ $`l_H^2+\frac{4}{3}l_{\mathrm{so}}^2.`$ (14c)
Equations (11) and (13) have the solution

$`k_0`$ $`=`$ $`{\displaystyle \frac{2Nl}{l_H}}\mathrm{coth}{\displaystyle \frac{L}{l_H}},`$ (15a)
$`h_0`$ $`=`$ $`{\displaystyle \frac{2Nl}{L}}\left(1+{\displaystyle \frac{l_H}{L}}\mathrm{coth}{\displaystyle \frac{L}{l_H}}\mathrm{coth}^2{\displaystyle \frac{L}{l_H}}\right).`$ (15b)
Expressions for $`k_1`$ and $`h_1`$ are obtained from Eq. (15) after the substitution $`l_Hl_H^{}`$. Substitution of $`h_0`$ and $`h_1`$ into Eq. (8) then allows for the calculation of the weak-localization correction to the conductance,

$$\delta g=\frac{l_H}{L}\mathrm{coth}\frac{L}{l_H}\frac{l_H^2}{L^2}3\left(\frac{l_H^{}}{L}\mathrm{coth}\frac{L}{l_H^{}}\frac{(l_H^{})^2}{L^2}\right).$$
(16)
At zero magnetic field, Eq. (16) simplifies to

$$\delta g=\frac{1}{3}+\frac{9l_{\mathrm{so}}^2}{4L^2}\frac{3l_{\mathrm{so}}\sqrt{3}}{2L}\mathrm{coth}\frac{2L}{l_{\mathrm{so}}\sqrt{3}}.$$
(17)
Equation (17) reproduces the limits $`\delta G=2e^2/3h`$ without spin-orbit scattering and $`\delta G=e^2/3h`$ with strong spin-orbit scattering. Without spin-orbit scattering, Eq. (16) agrees with the weak localization correction calculated in Ref. altshuler84, . For large magnetic fields, $`Ll_H`$, Eq. (16) simplifies to

$$\delta g=\frac{1}{L}\left(l_H3(l_H^2+\frac{4}{3}l_{\mathrm{so}}^2)^{1/2}\right),$$
(18)
which has the same functional form as the weak localization obtained using diagrammatic perturbation theory.altshuler81 ; aa ; santhanam84b ; santhanam87b  Comparison of Eq. (18) and Refs. altshuler81, ; aa, ; santhanam84b, ; santhanam87b,  allows us to identify $`l_{\mathrm{so}}`$ as the spin-orbit length, and, for a channel (with width $`Wl`$) in a two-dimensional electron gas in a perpendicular magnetic field $`B`$,

$$l_H^2=3(\mathrm{}/WBe)^2.$$
(19)
The case of a cylindrical wire of radius $`Rl`$ and magnetic field perpendicular to the wire is obtained by the substitution $`W^23R^2/2`$. For $`l>W`$ (or $`l>R`$) the crossover length $`l_H`$ has a more complicated $`l`$-dependent expression.BvH

Figure 1(a) shows $`\delta g`$ as a function of the magnetic field for several values of the spin-orbit coupling. In Fig. 1(b) we show $`\delta g`$ as a function of $`l_{\mathrm{so}}^1L`$ for several values of the magnetic field.

We now turn to a description of the weak localization correction in a non-homogeneous microbridge. Examples of non-homogeneous microbridges with varying widths are shown in the inset of Fig. 2. If the wire cross section or the electron density vary with the coordinate $`x`$ along the wire, the number of propagating channels at the Fermi level $`N`$ also varies with $`x`$. We assume that $`N(x)`$ has a minimum for $`x=0`$ and that $`dN/dx>0`$ ($`dN/dx<0`$) for all $`x>0`$ ($`x<0`$). Further, $`x`$-dependence of the impurity concentration, the smoothness of the boundary, the shape of the cross section, etc., causes an $`x`$-dependence of the length scales $`l`$, $`l_H`$, and $`l_{\mathrm{so}}`$.

The reflection matrix of the wire is constructed by building the wire from thin slices, starting at the narrowest point $`x=0`$. This way, the number of channels in the slices added to both ends of the wire can increase, but not decrease. For the construction of an evolution equation for the conductance $`g`$ and for the auxiliary functions $`h_0`$, $`h_1`$, $`k_0`$, and $`k_1`$, we distinguish between two types of added slices: A thin slice that contains a scattering site but for which the number of channels remains constant, and a thin slice without scatterer in which $`N`$ increases by unity. Addition of a slice of the former type causes a small change in the reflection matrix $`r`$, which is the same as for a quantum wire of constant thickness, see Eq. (5) above. Addition of a slice for which $`N`$ increases by unity does not cause a change of the conductance $`g`$ or of the auxiliary functions $`h_0`$, $`h_1`$, $`k_0`$, or $`k_1`$, as can seen by inspecting the cases $`x>0`$ and $`x<0`$ separately: For $`x>0`$, an increase of $`N`$ does not cause additional reflection, and hence does not affect the reflection matrix $`r`$; for $`x<0`$, an increment of $`N`$ changes the dimension of the reflection matrix $`r`$ by one,

$$r\left(\begin{array}{cc}r& 0\\ 0& 1\end{array}\right),$$
(20)
but does not change the conductance $`g`$, or the functions $`h_0`$, $`h_1`$, $`k_0`$, or $`k_1`$. Combining the two types of slices, we conclude that the only effect of the $`x`$-dependence of $`N`$ and $`l`$ is indirect, through the explicit appearance of $`N`$ and $`l`$ in statistics of the scattering matrix of the added slice, see Eq. (7). In the diffusive regime, $`N(x)`$ and $`l(x)`$ only appear in the combination

$$s(x)=N(x)l(x)/N_0l_0,$$
(21)
where $`N_0`$ and $`l_0`$ are number of propagating channels and mean free path at $`x=0`$. For large $`N`$ the function $`s(x)`$ may be considered continuous, and the evolution equations become differential equations which now include explicit reference to the function $`s(x)`$. If the wire length $`L`$ is replaced by the effective length $`\overline{L}`$,

$$\overline{L}=\frac{dx}{s(x)},$$
(22)
the evolution equations for $`g`$, $`h_0`$, $`h_1`$, $`k_0`$, and $`k_1`$ keep the same form as for homogeneous wires, provided we make the substitutions $`NN_0`$, $`L\overline{L}`$, $`ll_0`$, $`l_H\overline{l_H}=l_H/s(x)`$, and $`l_{\mathrm{so}}\overline{l_{\mathrm{so}}}=l_{\mathrm{so}}/s(x)`$.

The functional form of the leading-in-$`N`$ contribution to the conductance remains unchanged, $`G=(e^2/h)(2N_0l_0/\overline{L})`$. Also, for the limiting cases of no spin-orbit scattering and strong spin-orbit scattering, the weak localization correction $`\delta G`$ is still given by the universal result Eq. (1).BeenakkerMelsen94  However, because of the $`x`$-dependence of the length scales $`\overline{l_H}`$ and $`\overline{l_{\mathrm{so}}}`$, $`\delta g`$ acquires an explicit dependence on the shape of the disordered microbridge or the non-homogeneity of the mean free path or the electron density in the crossover region between the symmetry classes. For a large magnetic field ($`l_H^1L1`$), the weak-localization correction can be found in closed form,

$`\delta g`$ $`=`$ $`{\displaystyle \frac{1}{\overline{L}}}\left(l_{H,\mathrm{eff}}3l_{H,\mathrm{eff}}^{}\right),`$ (23)
$`l_{H,\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{\overline{L}}}{\displaystyle \frac{l_H(x)dx}{s(x)^2}},l_{H,\mathrm{eff}}^{}={\displaystyle \frac{1}{\overline{L}}}{\displaystyle \frac{l_H^{}(x)dx}{s(x)^2}}.`$
Equation (23) simplifies to Eq. (18) in the case of $`s(x)`$ constant. The same result follows if Eq. (18) is interpreted as a quantum interference correction to the one-dimensional resistivity and $`l_H`$ is taken $`x`$-dependent. For weaker magnetic fields with $`l_H^1L`$ of order unity, a numerical solution of the evolution equations is required.

In Fig. 2, we show results of a numerical solution of $`\delta g`$ for the examples $`s(x)`$ constant, $`s(x)=1+4|2x/L|`$ and $`s(x)=1+4(2x/L)^2`$, $`L/2<x<L/2`$. These functional forms correspond to diffusive microbridges in a two-dimensional electron gas of the form shown in the inset of Fig. 2 with uniform impurity concentration and mean free path $`lW`$. The three sets of curves in the figure represent strong, intermediate and weak spin-orbit scattering, respectively. For the intermediate case (middle set of curves in Fig. 2), three different values of $`l_{\mathrm{so}}`$ were chosen so that the weak-localization correction $`\delta g=0`$ is equal in the three cases for zero magnetic field. The magnetic field is characterized by the ratio $`l_{H,\mathrm{eff}}^1L`$, cf. Eq. (23), in order to remove a spurious shape dependence for the large-field asymptotes. While there is no dependence on the form of the function $`s(x)`$ in the limiting cases of zero and large magnetic fields, we observe that, indeed, $`\delta g`$ depends on the precise form of the non-homogeneity for intermediate magnetic field strengths, although, with proper scaling, the difference between the results for the three cases we considered is less than 10%.

In conclusion, we have shown that the scattering matrix approach to quasi one-dimensional weak localization can be used to obtain a detailed description of the crossover between the different universality classes. We have recovered some results known from diagrammatic perturbation theory, and have discovered one aspect of the problem that has not been noticed previously: The dependence of the functional form of the crossover on non-homogeneities in the conductor.

We thank V. Ambegaokar, N. W. Ashcroft, and D. Davidovic for discussions. This work was supported by the NSF under Grant Nos. DMR 0086509 and DMR 9988576, by the Packard Foundation, by the Natural Sciences and Engineering Research Council of Canada, and by the Dutch Science Foundation NWO/FOM.


text:
# The afterglow of GRB021211: Another case of reverse shock emission

## 1 Introduction

GRB021211 is a bright, long gamma-ray burst detected by HETE II on 2002 December 11 at 11:18:34 UT. The burst duration in the 8 - 40 keV band was $`>5.7`$ seconds, the fluence was about $`1\times 10^6`$ ergs cm-2 and the peak flux was $`>8\times 10^7`$ ergs cm-2 s-1 (Crew et al. 2002). This burst was also observed by Ulysses and Konus - Wind. As observed by Ulysses, it had a duration of about 15 seconds, a 25 - 100 keV fluence was approximately $`1.8\times 10^6`$ erg cm-2, and a peak flux was about $`4.5\times 10^7`$ erg cm-2 s-1 (Hurley et al. 2002). The spectroscopic observations of the optical afterglow identified three emission lines as \[OII\] 3727, and \[OIII\] 4959, 5007 at a redshift of $`z=0.800\pm 0.001`$ (Vreeswijk et al. 2002). Assuming $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`h=0.65`$, the isotropic gamma-ray energy is $`3.2\times 10^{51}`$ ergs.

The prompt localization of GRB021211 by HETE II allowed the rapid follow-up observation of the afterglow at very early time. Several groups had detected the optical emission shortly after the gamma-ray burst (Park et al. 2002; Li et al. 2002; Wozniak et al. 2002). The observations show that the optical flux declined steeply at early time, with a power-law index of about $``$1.6, while at later time the flux decayed with a slope of about $``$1, the break time is about 12 minutes after the burst (Chornock et al. 2002).

This break from a steep initial decline to a shallow later decline is similar to the early behavior of GRB990123 (Akerlof et al. 1999). The early emission of GRB990123 is believed to be due to the reverse shock, while the later emission is ascribed to the normal forward shock (Sari & Piran 1999a). Recently Kobayashi & Zhang (2002) have shown that the re-brightening in the GRB021004 optical afterglow light curve can be explained by the reverse shock and forward shock. Here we will show that this early break around 12 minutes after the burst can be interpreted as the superposition of the reverse shock and forward shock emission.

## 2 The emission from forward shock and reverse shock

### 2.1 Forward shock

Multiwavelength follow-up of gamma-ray burst afterglows has revolutionized GRB astronomy in recent years, yielding a wealth of information about the nature of GRBs. The observed properties of GRB afterglows are broadly consistent with models based on relativistic blast waves at cosmological distances (Meszaros & Rees 1997; Wijers et al. 1997). In the standard fireball models, the huge energy released by an explosion ($`10^{52}`$ ergs) is converted into kinetic energy of a shell expanding at ultra-relativistic speed. After the main GRB event occurred, the fireball continues to propagate into the surrounding gas, driving an ultra-relativistic blast wave (forward shock) into the ambient medium. The forward shock continuously heats fresh gas and accelerates relativistic electrons to very high energy, which produce the afterglow emission through synchrotron radiation.

Sari et al. (1998) have discussed the emission features of forward shock in great details. Using their results, we have

$$\nu _{\mathrm{m},\mathrm{f}}=5.1\times 10^{15}\left(1+z\right)^{1/2}ϵ_\mathrm{B}^{1/2}ϵ_\mathrm{e}^2g^2E_{52}^{1/2}t_\mathrm{d}^{3/2}\mathrm{Hz}$$
(1)
$$\nu _{\mathrm{c},\mathrm{f}}=2.7\times 10^{12}\left(1+z\right)^{1/2}ϵ_\mathrm{B}^{3/2}E_{52}^{1/2}n^1t_\mathrm{d}^{1/2}\mathrm{Hz}$$
(2)
$$t_{\mathrm{m},\mathrm{f}}=2.9\left(1+z\right)^{1/3}ϵ_\mathrm{B}^{1/3}ϵ_\mathrm{e}^{4/3}g^{4/3}E_{52}^{1/3}\nu _{\mathrm{R},15}^{2/3}\mathrm{days}$$
(3)
$$t_{\mathrm{c},\mathrm{f}}=7.3\times 10^6\left(1+z\right)^1ϵ_\mathrm{B}^3E_{52}^1n^2\nu _{\mathrm{R},15}^2\mathrm{days}$$
(4)
$$F_{\nu ,\mathrm{max},\mathrm{f}}=110\left(1+z\right)ϵ_\mathrm{B}^{1/2}E_{52}n^{1/2}D_{28}^2\mathrm{mJy}$$
(5)
where $`\nu _{\mathrm{m},\mathrm{f}}`$ is the typical synchrotron frequency of forward shock emission, $`\nu _{\mathrm{c},\mathrm{f}}`$ is the cooling frequency, $`ϵ_\mathrm{B}`$ and $`ϵ_\mathrm{e}`$ are the fractions of the shock energy transferred to the magnetic field and electrons, $`g=(p2)/(p1)`$, $`p`$ is the index of electron energy distribution, $`E_{52}`$ is the burst energy in units of $`10^{52}`$ ergs, $`t_\mathrm{d}`$ is the observer’s time in units of 1 day, $`n`$ is the surrounding medium density in units of 1 atom cm-3, $`\nu _{\mathrm{R},15}=\nu _\mathrm{R}/10^{15}\mathrm{Hz}`$, $`D_{28}`$ is the luminosity distance in units of $`10^{28}`$ cm, $`t_{\mathrm{m},\mathrm{f}}`$ ($`t_{\mathrm{c},\mathrm{f}}`$) is the time when the frequency $`\nu _{\mathrm{m},\mathrm{f}}`$ ($`\nu _{\mathrm{c},\mathrm{f}}`$) crosses the observed optical frequency, and $`F_{\nu ,\mathrm{max},\mathrm{f}}`$ is the peak flux.

According to the fireball model, before the peak time $`t_{\mathrm{m},\mathrm{f}}`$, the observed optical flux is expected to increase as $`F_\nu t^{1/2}`$, when the typical synchrotron frequency crosses the observed optical band, the flux reaches the maximum flux $`F_{\nu ,\mathrm{max},\mathrm{f}}`$, and then when $`t>t_{\mathrm{m},\mathrm{f}}`$, the flux decays as $`F_\nu t^{3(p1)/4}`$.

### 2.2 Reverse shock

The emission of reverse shock have been discussed by Meszaros & Rees (1997) and Sari & Piran (1999b). This shock heats up the shell’s matter and accelerates its electrons, then these electrons loss energy through synchrotron radiation. The reverse shock and the forward shock are separated by a contact discontinuity, across which the pressure is equal, so the energy density in both shocked regions is the same, therefore the total energy in both shocks is comparable.

The typical synchrotron frequency $`\nu _{\mathrm{m},\mathrm{r}}`$ and cooling frequency $`\nu _{\mathrm{c},\mathrm{r}}`$ of the reverse shock at the time when it crosses the shell can be easily calculated by comparing them to those of the forward shock. Since at the shock crossing time ($`t_\mathrm{A}`$), the reverse shock and the forward shock have the same Lorentz factor and energy density (which suggests the magnetic fields in both shocked regions are the same, if we assume that the magnetic equipartition factor is same in both regions), so the cooling frequency of the reverse shock $`\nu _{\mathrm{c},\mathrm{r}}`$ is the same as that of the forward shock,

$$\nu _{\mathrm{c},\mathrm{r}}(t_\mathrm{A})\nu _{\mathrm{c},\mathrm{f}}(t_\mathrm{A})$$
(6)
The typical synchrotron frequency is proportional to the electrons random Lorentz factor squared and to the magnetic field and to the Lorentz boost. The Lorentz boost and the magnetic field are the same for both the reverse shock and forward shock, while the random Lorentz factor of reverse shock is $`\gamma _0/\gamma _\mathrm{A}`$ compared to $`\gamma _\mathrm{A}`$ of the forward shock, i.e. the effective temperature of reverse shock is much lower than that of the forward shock (by a factor $`\gamma _\mathrm{A}^2/\gamma _0`$, where $`\gamma _0`$ is the initial Lorentz factor and $`\gamma _\mathrm{A}`$ is the Lorentz factor at the crossing time.), then the typical synchrotron frequency of reverse shock at the crossing time is

$`\nu _{\mathrm{m},\mathrm{r}}(t_\mathrm{A})`$ $``$ $`{\displaystyle \frac{\gamma _0^2}{\gamma _\mathrm{A}^4}}\nu _{\mathrm{m},\mathrm{f}}(t_\mathrm{A})`$ (7)
$``$ $`3.5\times 10^{15}ϵ_\mathrm{B}^{1/2}\left({\displaystyle \frac{ϵ_\mathrm{e}}{0.1}}\right)^2g^2\left({\displaystyle \frac{\gamma _0}{300}}\right)^2n^{1/2}\mathrm{Hz}`$ (8)
The peak flux at the typical frequency is proportional to the number of electrons, the magnetic field and the Lorentz boost. The magnetic fields and the Lorentz boost are the same for both reverse and forward shock, while at the crossing time the mass of the shell is larger by a factor of $`\gamma _\mathrm{A}^2/\gamma _0`$ than that of the ambient medium swept by the forward shock (Kobayashi & Zhang 2002), so we have

$`F_{\nu ,\mathrm{max},\mathrm{r}}(t_\mathrm{A})`$ $``$ $`{\displaystyle \frac{\gamma _\mathrm{A}^2}{\gamma _0}}F_{\nu ,\mathrm{max},\mathrm{f}}`$ (9)
$``$ $`110\left(1+z\right){\displaystyle \frac{\gamma _\mathrm{A}^2}{\gamma _0}}ϵ_\mathrm{B}^{1/2}E_{52}n^{1/2}D_{28}^2\mathrm{mJy}`$ (10)
After the reverse shock has passed through the ejecta, the ejecta cools adiabatically. Sari & Piran (1999a) have shown that, for $`t>t_\mathrm{A}`$, the particle Lorentz factor evolves as $`\gamma _\mathrm{e}t^{13/48}`$, the emission frequency drops quickly with time according to $`\nu _\mathrm{e}t^{73/48}`$, and the peak flux falls like $`F_{\nu _\mathrm{e}}t^{47/48}`$. So if the optical band $`\nu _\mathrm{R}`$ is below the typical synchrotron frequency $`\nu _{\mathrm{m},\mathrm{r}}`$, then the flux decays as $`F_\nu t^{17/36}`$, while when $`\nu _\mathrm{R}>\nu _{\mathrm{m},\mathrm{r}}`$, the flux decreases as $`F_\nu t^{(21+73p)/96}`$.

## 3 Fitting the afterglow of GRB021211

Using the emission features of reverse shock and forward shock described above, we can fit the optical light curve of GRB021211. Here we take the values $`z=0.8`$, $`E=3.2\times 10^{51}`$ ergs, and $`p=2.3`$.

For the forward shock emission, the observed optical flux is

$$\frac{F_{\nu ,\mathrm{f}}(t)}{F_{\nu ,\mathrm{max},\mathrm{f}}}=\{\begin{array}{cc}\left(\frac{t}{t_{\mathrm{m},\mathrm{f}}}\right)^{1/2}\hfill & \mathrm{for}t<t_{\mathrm{m},\mathrm{f}}\hfill \\ \left(\frac{t}{t_{\mathrm{m},\mathrm{f}}}\right)^{3(p1)/4}\hfill & \mathrm{for}t_{\mathrm{m},\mathrm{f}}<t<t_{\mathrm{c},\mathrm{f}}\hfill \end{array}$$
(11)
Using eqs.(3)(5) and (11) we can give the afterglow light curve from the forward shock, as shown in Fig.1 by the dashed line. From fitting the observed data we can obtain the relation

$$ϵ_\mathrm{B}\left(\frac{ϵ_\mathrm{e}}{0.1}\right)^{8/5}n^{3/5}9.1\times 10^4$$
(12)
In addition, the observation implies that $`t_{\mathrm{m},\mathrm{f}}`$ should be less than 100 seconds (if $`t_{\mathrm{m},\mathrm{f}}>100`$ s, there will be a bump in the afterglow light curve), and $`t_{\mathrm{c},\mathrm{f}}`$ should be larger than 1 day (otherwise there will be a steepening of the light curve) , so from eqs.(3)(4) we have

$$ϵ_\mathrm{B}^{1/4}\left(\frac{ϵ_\mathrm{e}}{0.1}\right)0.14$$
(13)
$$ϵ_\mathrm{B}n^{2/3}0.04$$
(14)
For the reverse shock emission, for $`t>t_\mathrm{A}`$, we have the relations $`\nu _{\mathrm{m},\mathrm{r}}(t)=\nu _{\mathrm{m},\mathrm{r}}(t_\mathrm{A})\left(\frac{t}{t_\mathrm{A}}\right)^{73/48}`$, $`F_{\nu ,\mathrm{max},\mathrm{r}}(t)=F_{\nu ,\mathrm{max},\mathrm{r}}(t_\mathrm{A})\left(\frac{t}{t_\mathrm{A}}\right)^{47/48}`$, then the observed flux can be written as

$`F_{\nu ,\mathrm{r}}(t)`$ $`=`$ $`F_{\nu ,\mathrm{max},\mathrm{r}}(t)\left[{\displaystyle \frac{\nu }{\nu _{\mathrm{m},\mathrm{r}}(t)}}\right]^{(p1)/2}`$ (15)
$`=`$ $`F_{\nu ,\mathrm{max},\mathrm{r}}(t_\mathrm{A})\left[{\displaystyle \frac{\nu }{\nu _{\mathrm{m},\mathrm{r}}(t_\mathrm{A})}}\right]^{\frac{(p1)}{2}}\left({\displaystyle \frac{t}{t_\mathrm{A}}}\right)^{\frac{73p+21}{96}}`$ (16)
Using eqs.(8)(10) and (16) we can give the afterglow light curve from the reverse shock, as shown in Fig.1 by the dotted line. From fitting we can obtain the relation

$$\left(\frac{\gamma _0}{300}\right)^{23/10}\left(\frac{\gamma _\mathrm{A}}{\gamma _0}\right)^2ϵ_\mathrm{B}^{33/40}\left(\frac{ϵ_\mathrm{e}}{0.1}\right)^{13/10}n^{33/40}t_\mathrm{A}^29.8$$
(17)
Combining eqs. (12)(14) and (17), we get

$$\left(\frac{\gamma _0}{300}\right)\left(\frac{\gamma _\mathrm{A}}{\gamma _0}\right)^{20/23}31ϵ_\mathrm{B}^{1/184}n^{27/184}t_\mathrm{A}^{20/23}$$
(18)
$$ϵ_\mathrm{B}^{1/16}\left(\frac{ϵ_\mathrm{e}}{0.1}\right)0.077$$
(19)
Therefore eqs.(12)(13) and (19) give the constraint on the parameters $`ϵ_\mathrm{B}`$, $`ϵ_\mathrm{e}`$ and $`n`$. Fig.2 shows the relation between $`ϵ_\mathrm{B}`$, $`ϵ_\mathrm{e}`$ and $`n`$. The dotted, dash-dotted, dashed and dot-dot-dashed lines represent $`n=`$0.1, 1, 10 and 0.0016 respectively. We find that the allowed values of $`ϵ_\mathrm{B}`$ and $`ϵ_\mathrm{e}`$ lie in the region confined by two lines Lc1 (eq.(19)) and Lc2 (eq.(13)). It is obvious that $`n`$ must be larger than 0.0016, and $`ϵ_\mathrm{e}`$ must be larger than 0.0077. If we take $`n=1`$, $`ϵ_\mathrm{e}=0.07`$, then $`ϵ_\mathrm{B}=1.6\times 10^3`$. We propose that more observations are needed in order to further estimate the values of $`ϵ_\mathrm{B}`$, $`ϵ_\mathrm{e}`$ and $`n`$.

From eq.(18) we see that the initial Lorentz factor $`\gamma _0`$ depends on $`ϵ_\mathrm{B}`$ and $`n`$ very weakly, so as an approximation, and taking $`\gamma _\mathrm{A}\gamma _0`$, then we have $`\gamma _09300t_\mathrm{A}^{20/23}`$. Since the duration is about 15s and the first observation time is 65s after the burst, so the value of $`t_\mathrm{A}`$ should lie between 15s and 65s, and therefore we can get the initial Lorentz factor $`250<\gamma _0<900`$, which is consistent with the lower limit estimates base on the $`\gamma `$ \- $`\gamma `$ attenuation calculation (Fenimore et al. 1993).

## 4 Discussion and conclusion

The current afterglow observations usually detect radiation several hours after the burst, at this stage the Lorentz factor is independent of the initial Lorentz factor, thus these observations do not provide useful information on the initial extreme relativistic motion. The initial Lorentz factor is a very important quantity for constraining the GRB models since it specifies how ”clean” the fireball is. Therefore to detect the early afterglow of GRBs is very important, since it can provide the information on the initial Lorentz factor. It is fortunately that the early afterglow of GRB021211 have been observed, by fitting its optical light curve we obtain its initial Lorentz factor $`250<\gamma _0<900`$, this value seems reasonable since it is widely believed that the initial fireball Lorentz factor should be larger than 100 in order to avoid photon-photon attenuation. Of course, further constraint on the value of $`\gamma _0`$ needs more early afterglow observations.

GRB990123 is the first burst for which its optical flash was observed, the peak flux was about 1 Jy in R-band. After that many efforts have been made to try to find the optical flash from other GRBs, but only upper limits are given (Akerlof et al. 2000). Here we also note that the optical flux of GRB021211 is substantially fainter than 990123 at similar epochs. Why GRB990123 is so bright? One reason may be that GRB990123 is a very bright burst, so its reverse shock emission is also very strong. On the other hand, from fitting we note that the values of $`ϵ_\mathrm{B}`$ and $`ϵ_\mathrm{e}`$ of GRB021211 are somewhat smaller, which leads to the fact that the typical synchrotron frequency of reverse shock is well below the optical band, so the early afterglow (or optical flash) is weak. While for GRB990123 the typical synchrotron frequency of reverse shock is close to the optical band (Sari & Piran 1999a; Kobayashi & Zhang 2002).

For smaller values of $`ϵ_\mathrm{B}`$ and $`ϵ_\mathrm{e}`$, not only the typical synchrotron frequency of reverse shock is small, but also the typical synchrotron frequency of forward shock is small, so the time $`t_{\mathrm{m},\mathrm{f}}`$ when the typical frequency of forward shock crosses the optical band is also small, for GRB021211, the observations required $`t_{\mathrm{m},\mathrm{f}}100`$s. The late time afterglow for $`t>t_{\mathrm{m},\mathrm{f}}`$ is $`F_\nu =F_{\nu ,\mathrm{max},\mathrm{f}}\left(t/t_{\mathrm{m},\mathrm{f}}\right)^{3(p1)/4}`$, so for smaller value of $`t_{\mathrm{m},\mathrm{f}}`$, the observed optical flux should be much fainter than those with larger values of $`t_{\mathrm{m},\mathrm{f}}`$, so we suggest that the so-called dark bursts whose afterglow have not been observed might be due to their very small values of $`t_{\mathrm{m},\mathrm{f}}`$.

The early afterglow of GRB021211 shows that there is a early break in its optical light curve, before the break time the flux declined with a power-law index of about $``$1.6, while at later time the flux decayed with a slope of about $``$1. Although the reverse shock model predicts that the optical flux should decay with a power-law index of about $``$2, here we show that the superposition of both the forward shock and the reverse shock emission can well account for the observed light curve. Therefore we suggest that this early break might be a common feature in early optical afterglow, and before the break time the slope of flux decline may be flatter than $``$2 since it contains the contribution from both the reverse shock and the forward shock emission.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation (grants 10073022 and 10225314) and the National 973 Project on Fundamental Researches of China (NKBRSF G19990754).