text:
# The Magnetohydrodynamic Kelvin-Helmholtz Instability: A Three-Dimensional Study of Nonlinear EvolutionThe Astrophysical Journal, in press

## 1 INTRODUCTION

Strongly sheared boundary flows are ubiquitous in astrophysical environments as diverse as the earth’s magnetopause and supersonic jets. The susceptibility of such boundaries to the Kelvin-Helmholtz (K-H) instability is well-known (e.g., Chandrasekhar 1961). Development of the instability may lead to turbulence, momentum and energy transport, dissipation and mixing of fluids (see, e.g., Maslowe 1985 for a review).

Most astrophysical environments are electrically conducting, so relevant fluids are likely to be magnetized on length and time scales of common interest. Thus, it is important to understand the role of magnetic fields in the K-H instability. The basic linear stability analysis of the magnetohydrodynamic (MHD) K-H instability was carried out long ago (e.g., Chandrasekhar 1961; Miura & Prichett 1982). There is now also a growing literature of the nonlinear evolution of the MHD K-H instability beginning from a variety of possible initial flow configurations, at least in the earlier evolution stages in two dimensions (2D) (e.g., Tajima & Leboeuf 1980; Wang & Robertson 1984; Miura 1984, 1987, 1997; Wu 1986; Dahlburg et al. 1997; Keppens et al. 1999; Keller & Lysak 1999). Fully three dimensional (3D) nonlinear studies are still quite limited and so far have generally not followed flow evolution to anything resembling a final state. They do show that full coupling of magnetic field and flow in the third dimension may quickly introduce obvious dynamical effects, however (e.g., Galinsky & Sonnerup 1994; Keppens & Tóth 1999; Keller, Lysak & Song 1999). The 3D simulations reported below, on the other hand, were continued over many dynamical time scales, so that the ultimate relaxed states for the flows are clear.

Strong magnetic fields, through their tension, are well known to stabilize the K-H instability. However, the considerable potential for much weaker fields to modify the nonlinear instability, and, in particular, to reorganize the subsequent flow, has only recently been emphasized. Malagoli et al. (1996), Frank et al. (1996), Jones et al. (1997) and Jeong et al. (2000) have carried out high resolution 2D and $`2\frac{1}{2}`$D numerical MHD simulations of the full nonlinear evolution of the K-H instability for periodic sections of 2D flows. They have demonstrated clearly that an initially weak magnetic field can fundamentally alter evolution of the K-H instability, either by disrupting the 2D hydrodynamical (HD) flow character or by enhancing dissipation during the nonlinear evolution of this instability.

While specific issues may surely depend on matching details of a simulated configuration to the physical situation imagined, the most basic insights often come from very simple, idealized model flows. In our review here we follow Frank et al. (1996) and Jones et al. (1997) in considering periodic sections of flows. Let us define the 2D computational plane as the $`xy`$ plane, so that for now there is an assumed invariance along the $`z`$ direction. Initially, two uniform, but oppositely directed, velocity fields along the $`x`$direction are separated by a thin, smooth trans-sonic shear layer ($`M_0=U_0/c_s=1`$, with $`U_0`$ the velocity difference across the shear layer). The magnetic field, which is initially aligned within the $`xz`$ plane, is uniform, and so is density. The most important parameter in predicting the outcome of the subsequent MHD K-H instability is the following Alfvénic Mach number of the velocity transition; namely, $`M_{Ax}=U_0/c_{Ax}`$, with $`c_{Ax}=B_x/\sqrt{\rho }`$ the projected Alfvén speed, and $`B_x`$ the magnetic field component aligned with the flow in units giving magnetic pressure, $`p_b=B^2/2`$. We note that Jones et al. (1997) demonstrated when $`M_{Ax}2`$ (for which the instability is not stabilized by magnetic tension) and the initial magnetic field is uniform that the existence of a finite $`B_z`$ is largely irrelevant to evolution of the K-H instability; that is, the orientation of the field in the $`xz`$ plane does not, by itself, matter, except through its isotropic pressure.

When $`M_{Ax}4`$ but with the initial magnetic field aligned with the flow, the field is a little too weak to provide linear stability ($`M_{Ax}2`$). Modest growth of corrugations along the perturbed shear layer generates sufficient magnetic tension to prevent further development of nonlinear evolution. That is, the flow is nonlinearly stabilized. However, when $`M_{Ax}4`$, the magnetic field is too weak to have much, if any, apparent influence during the linear and early 2D nonlinear phases of the K-H instability. Thus, the initial development is largely HD. So, a Kelvin’s Cat’s Eye vortex forms with its axis in the shear layer, but perpendicular to the initial flow. In 2D HD, this structure is stable. Jones et al. (1997), therefore, chose $`M_{Ax}4`$ as a convenient boundary between strong and weak magnetic field behaviors in the 2D MHD K-H instability.

Within the weak field regime, it is also possible to distinguish further two qualitatively different evolutions. Unless an initially uniform field is amplified sufficiently during one rotation of the Cat’s Eye to reduce $`M_{Ax}`$ to values of order unity along the vortex perimeter, magnetic stresses have little immediate dynamical influence on vortex evolution. Then, the magnetic field primarily serves to enhance dissipation of kinetic energy through expulsion of magnetic flux in the $`xy`$ plane (a.k.a. flux annihilation, via tearing mode reconnection). This case of very weak field was called dissipative in Jones et al. (1997). In the discussion below we will label flows initiated in this regime by VWF. The more interesting regime is that where the initial field is too weak to prevent formation of the Cat’s Eye, but strong enough that $`M_{Ax}1`$ at some locations within the Cat’s Eye by the end of a single vortex rotation. Under those circumstances, relaxation of magnetic stresses during reconnection deforms and then disrupts the Cat’s Eye. This was called the weak field regime or the disruptive regime in Jones et al. (1997). We will label these cases below as WF.

In a 2D flow, a magnetic field line is stretched by about an order of magnitude while becoming wrapped around a vortex that it once spanned. That reduces $`M_{Ax}`$ by a similar factor, since the field strength increases proportionally to the length of a flux tube (e.g., Gregori et al. 2000). Thus, it turned out in 2D $`M_{Ax}20`$ is the boundary between the two cases. However, even for $`M_{Ax}>20`$, there can be gradual disruption of the Cat’s Eye at late epoch by an accumulation of small effects from Maxwell stresses. So, this dividing line is not distinct.

In the disruptive, WF case, there is also a dynamical alignment between the magnetic and velocity fields during reconnection along the perimeter of the Cat’s Eye, and local cross helicity ($`|vB|`$) is maximized. In this configuration the 2D flow returns to a laminar form, but is now stable to perturbations smaller than the size of the computational box. Jones et al. (1997) emphasized that 2D vortex disruption was magnetically driven, despite the fact that the initial $`\beta =p_g/p_b>>1`$, where $`p_g`$ is the thermal gas pressure. We note this, since it is very common to ignore dynamical influences from magnetic fields under the condition $`\beta =p_g/p_b>>1`$. That measures only the relative influences of pressure gradients, not the full Maxwell stresses. The Alfvén Mach number, on the other hand, compares more closely Maxwell to Reynolds stresses, so should provide a more direct measure of the immediate dynamical consequences of the magnetic field in nonequilibrium flows. Certainly, that is the case here.

Our objective now is to extend those previous 2D results into fully 3D flows. This step is important, since it is already well established that the Cat’s Eye structure so prominent and stable in plane-symmetric, 2D flows coming from the HD K-H instability is unstable to perturbations along its axis in 3D (Hussain 1984; Bayly 1986; Craik & Criminale 1986). The resultant HD flow becomes turbulent (e.g., Maslowe 1985). So, we should ask how a weak magnetic field will modify that outcome. In addition, since both vorticity and magnetic flux is subject to stretching in 3D but not in 2D, and, since vortex stretching profoundly changes 3D flows when compared to those in 2D, we might expect to find that as soon as the flow evolution deviates from 2D character the characterizations listed above no longer apply. We will find, in fact, that they do still apply, but the domain of initial magnetic field strengths that can significantly influence the flow evolution is extended to weaker fields in 3D. We will also see that the morphologies and statistical properties of magnetic and flow structures expected during 3D nonlinear flow evolution depend on the strength of the initial magnetic field. A preliminary report on some of these calculations is contained in Jones et al. (1999), which also includes some useful animations on a CD ROM. Those same animations are presently posted on the web site: http://www.msi.umn.edu/$``$twj/research/mhdkh3d/nap98.html. The plan of the present paper is as follows: In §2 we will summarize the problem set-up and numerical method. §3 contains detailed discussions of results. A brief summary and conclusion follow in §4.

## 2 THE PROBLEM

The simulations reported here are direct extensions of the Jones et al. (1997) study to fully 3D flows. The equations we solve numerically are those of ideal compressible MHD, where the displacement current and the separation between ions and electrons are neglected as well as the effects of viscosity, electrical resistivity and thermal conductivity. In conservative form, the equations are

$$\frac{\rho }{t}+\stackrel{}{}\left(\rho \stackrel{}{v}\right)=0,$$
(1)
$$\frac{(\rho \stackrel{}{v})}{t}+_j\left(\rho \stackrel{}{v}v_j\stackrel{}{B}B_j\right)+\left(p+\frac{1}{2}B^2\right)=0,$$
(2)
$$\frac{E}{t}+\stackrel{}{}\left[\left(E+p+\frac{1}{2}B^2\right)\stackrel{}{v}\left(\stackrel{}{v}\stackrel{}{B}\right)\stackrel{}{B}\right]=0,$$
(3)
$$\frac{\stackrel{}{B}}{t}+_j\left(\stackrel{}{B}v_j\stackrel{}{v}B_j\right)=0,$$
(4)
along with the constraint $`\stackrel{}{}\stackrel{}{B}=0`$ imposed to account for the absence of magnetic monopoles (e.g., Priest 1984). Gas gas pressure is given by

$$p=\left(\gamma 1\right)\left(E\frac{1}{2}\rho v^2\frac{1}{2}B^2\right)_.$$
(5)
Standard symbols are used for common quantities. The magnetic pressure is $`p_b=B^2/2`$ and the Alfvén speed is $`c_A=B/\sqrt{\rho }`$.

The simulations have been carried out in a cubic computational box of length $`L_x=L_y=L_z=L=1`$. Boundaries are periodic in the directions contained within the shear layer (namely, $`x`$ and $`z`$) and reflecting above and below the shear layer (namely, $`y`$). As before we have simulated flows that are initially uniform except for a hyperbolic tangent velocity shear layer in the $`y`$-coordinate, given as

$$\stackrel{}{v}_0=\frac{U_0}{2}\mathrm{tanh}\left(\frac{yL_y/2}{a}\right)\widehat{x}$$
(6)
with $`a=L/25`$. The equilibrium flow is directed in the $`x`$ direction for $`y>0.5`$, and $`+x`$ direction for $`y<0.5`$. The velocity difference across the shear layer is unity, $`U_0=1`$. The sonic Mach number of the transition is unity, $`M_s=1`$, and the adiabatic index, $`\gamma =5/3`$. With this configuration K-H unstable modes will have zero phase velocities in the computational reference frame. Whereas an initial $`B_z`$ has no appreciable influence on $`2\frac{1}{2}`$D flows, because field lines could not be stretched in that dimension, we expect field line stretching in this dimension to be important in 3D. Thus, the initial magnetic field is oblique to the flow direction, with $`\theta =30\mathrm{°}`$, but parallel to the shear plane with strengths corresponding to $`M_{Ax}=2.5,5,14.3,50,143,500,`$ and $`1.43\times 10^3`$. For comparison, $`\beta =(2/\gamma )(M_{Ax}\mathrm{cos}\theta \mathrm{°}/M_s)^2`$ here. See Table 1 for further details. A random perturbation of small amplitude has been added to the velocity to initiate the instability.

All cases have been simulated with grids having $`64^3`$ and $`128^3`$ zones (labeled ‘$`l`$’ standing for ‘low’ and ‘m’ standing for ‘medium’, respectively in Table 1) to explore basic properties including resolution issues. For three representative cases with $`M_{Ax}=14.3,50,1.43\times 10^3`$, the calculations have been repeated again with $`256^3`$ zones (labeled ‘h’ standing for ‘high’ in Table 1).

Each simulation has been run up to time $`t=2050`$ (see Figs. 3 and 8 for the end time). For comparison, the sound crossing time for the box is unity. With the initial perturbation applied, the Cat’s Eye forms by about $`t6`$ in those cases where it develops. The nominal subsequent turnover time for the Cat’s Eye is also $`t=t_e6`$. We note for reference that our computing time units here are 2.51 longer than those in Frank et al. (1996) and Jones et al. (1997), since $`L=2.51`$ was set to there to match earlier papers, but they are same as those in Jeong et al. (2000). To aid comparison we mention that in our present units the growth times associated with the modes having wavelength $`\lambda =L`$, the box size, would typically be $`t_g0.60.7`$. Thus, our simulations extend typically $`3570`$ K-H linear growth times for such modes. Here, we have applied random velocity perturbations to the initial equilibrium in the present simulations. Hence, modes with shorter wavelengths develop first if unstable. However, they generally merge on time scales comparable to the growth time of the mode with wavelength which is the sum of the wavelengths of merged modes. So the mentioned growth time is still a reasonable estimate of the time required for instabilities to become a significant influence on the flow.

The ideal MHD equations have been solved using a multi-dimensional MHD code based on the explicit, finite difference “Total Variation Diminishing” or “TVD” scheme. That method is an MHD extension of the second-order finite-difference, upwinded, conservative gasdynamics scheme of Harten (1983), as described by Ryu & Jones (1995). The multi-dimensional version of the code, along with a description of various one and 2D flow tests is contained in Ryu, Jones, & Frank (1995). This version of the code contains a fast Fourier transform-based “flux cleaning” routine that maintains the $`\stackrel{}{}\stackrel{}{B}=0`$ condition at each time step within machine accuracy.

## 3 RESULTS

The seven cases listed in Table 1 include examples that in 2D exhibit “dissipative”, “disruptive” and “nonlinear stabilizing” behaviors, using the terminology defined in §1. Cases that are linearly stable were not considered, since we expect no new behaviors. Cases 1-4, with $`50M_{Ax}1.43\times 10^3`$, would in 2D have been dissipative in character, Cases 5-6, with $`5M_{Ax}14.3`$, would have been disruptive, while Case 7, with $`M_{Ax}=2.5`$, would have been nonlinearly stabilizing.

Initially all the cases evolve in ways consistent with the 2D description in §1 for the same field strength. For Case 7 that is basically the end of the story, since the flow is nonlinearly stabilized and remains laminar through the entire simulation. Just as in the 2D simulation described by Frank et al. (1996) the final result is a broadened shear layer that has been subjected to very minimal kinetic energy dissipation.

In all the remaining cases a Cat’s Eye develops and is subsequently destroyed. The Cat’s Eye is a 2D structure, so up to the point of its formation all the flows are still quasi-2D. Recall that in 2D the Cat’s Eye remained stable in the very weak field, VWF case, since Maxwell stresses were not built up enough to disrupt it. But in 3D this structure is HD unstable (see $`\mathrm{\S }`$3.1). However, by and large, we still find the flow characters of dissipative (for the WF case) and disruptive (for the VWF case) carry over into 3D flows, as described in detail in §3.1 and 3.2 below. One interesting deviation is that the 3D magnetic field in Case 4, with $`M_{Ax}=50`$, is significantly more disruptive than it would be in 2D. Case 4 would have been defined as dissipative in 2D, and it exhibits similar properties in our low and medium resolution 3D simulations. On the other hand, in higher resolution 2D simulations it showed some evidence for long-term disruptive magnetic field influence through accumulated small flow distortions. Those tendencies are much more consequential in 3D, so, we will describe the in detail the 3D behaviors observed in that case in §3.3.

### 3.1 Very Weak Field (VWF) Cases: Turbulence

Since they come close to HD behaviors, and thus offer a useful benchmark, we begin with discussion of Cases 1, 2 and 3 ($`M_{Ax}=1.43\times 10^3,500,143`$), which would all fall under the VWF or dissipative descriptions in 2D. Fig. 1 shows at two times the spatial distributions of magnetic field strength ($`|B|`$) and vorticity magnitude ($`|\omega |`$) for Case 1$`\mathrm{}`$ and Case 1m, while Fig. 2 shows at three times the same information for the analogous high resolution, Case 1h. One can see that the global behaviors are qualitatively consistent in all three simulations. Quantitatively, the comparisons are quite similar to what we described previously in Frank et al. (1996) and in Jones et al. (1997). That is, as expected, simulations with higher resolution capture finer structures.

In 3D HD flow the Cat’s Eye violently breaks up within approximately one eddy turnover, and the flow becomes highly disordered, with very little evidence of the initial shear field; i.e., decaying isotropic turbulence develops (e.g., Maslowe 1985). For the VWF MHD cases we studied (Cases 1, 2, and 3) this behavior is also seen. We can identify the root causes for the disruption of the Cat’s Eye as follows: First, Hussain (1984) pointed out the importance of the growth of coherent vortex tubes that span the Cat’s Eye. These features, called rib vortices by Hussain, are clearly present in the early snapshots of Figs. 1 and 2. Hussain pointed out that the ribs are anchored in saddle points within the flow at the ends of the Cat’s Eye, so they are subjected to rapid and intense vortex stretching. This leads to non-axial stresses on the flows. A second effect is fluid elements caught in the Cat’s Eye vortex move along non-circular, or roughly elliptical, paths, so that they feel time varying shear forces. Such fluid elements are known to be subject to the elliptical instability, when motion perpendicular to the elliptical path is allowed in 3D (Bayly 1986; Craik & Criminale 1986). Together these effects unstably distort the initially 2D character of the flow, so that the Cat’s Eye breaks up violently in less than a single turnover. Fig. 2 shows that by $`t=6`$ the interior of the Cat’s Eye in Case 1h is filled with a knotted tangle of thin vortex tubes. They are the remnants of smaller size vortices, which are developed and gone through the processes of disruption and merger earlier. Fig. 2 also illustrates for Case 1h that between $`t=6`$ and $`t=8`$ the Cat’s Eye has already been severely distorted. By $`t=20`$ the entire flow pattern has broken down into an apparently isotropic distribution of vortex tubes.

At this point it is useful to look closely at the simultaneous evolution of the magnetic field. Fig. 2 shows us that the regions of strong magnetic field generally match the regions of strong vorticity (although there are vortex tubes which are the remnants of the initial vorticity in the problem or the early activity inside the Cat’s Eye as described more in below, and so do not match with strong magnetic field). Note that magnetic fields in this case are essentially passive, and the flows are nearly ideal. It is well-known that under these circumstances the magnetic field and vorticity obey the same evolution equation (e.g., Shu 1992). So the coincidence of strong magnetic field regions with strong vorticity region is, to a certain degree, expected. We emphasize, however, that the coincidence does not necessarily imply that the two vector fields are aligned. In fact, that is distinctly not the case along the rib vortices as they first form. Rather, magnetic field lines passing diagonally across the shear layer are initially stretched around the forming Cat’s Eye and become embedded within the flow fields of the rib vortices. At first the fields in the rib vortices are simply stretched over the ribs, but, over the course of the Cat’s Eye formation, those field lines become twisted around the ribs, like twisted-pair electrical wires. This effect enhances significantly the amount of stretching those lines undergo compared to their 2D analogs. Thus, the magnetic fields embedded in the rib vortices are substantially stronger than other field lines merely stretched around the Cat’s Eye perimeter. That explains the coincidence of strong magnetic field regions with strong vorticity region at $`t=20`$ in Fig. 2. It also enhances the role of weak magnetic fields in 3D compared to 2D, as we shall see.

As magnetic field lines become twisted around rib vortices they soon develop topologies unstable to reconnection (e.g., Lysak & Song 1990), however. The reconnection preserves helicity (e.g., Ruzmaikin & Akhmetiev 1994), although other topological field measures, such as twist, writhe and kink that contribute to helicity may change (e.g., Berger & Field 1984; Bazdenkov & Sato 1998). The product by the time the Cat’s Eye begins to be disrupted is a set of twisted magnetic flux tubes around the perimeter of the Cat’s Eye that do align themselves with the original rib vortices; i.e.,in the $`xy`$ plane.

Simultaneously, the complex motions in the Cat’s Eye interior, which were mentioned above, lead to extensive magnetic reconnection events that produce a 3D version of magnetic flux expulsion mentioned in §1 (Weiss 1966; see Jones et al. 1997 for discussion of that process). The product, when the Cat’s Eye begins to break up, is a region of very weak and tangled magnetic field inside the Cat’s Eye that on average trends in the $`z`$ direction. The total magnetic flux through the full computational box is constant, of course, but this series of events has separated the magnetic flux embedded in the Cat’s Eye into relatively strong flux tubes perpendicular to the axis of the Cat’s Eye and wrapped around it, plus largely disordered magnetic flux inside the Cat’s Eye with a mean field aligned with the axis of the Cat’s Eye. This dichotomy is retained in the magnetic field structures at the end of end of Case 4h, in fact, as we will address that in §3.3.

The evolution of energy partitioning is illustrated for the VWF cases along with Case 4 in Fig. 3. Keep in mind that because we use a periodic box in the $`x`$ and $`z`$ directions and hard walls in the $`y`$ direction, the system is effectively closed and total energy is conserved. It is, therefore, a necessarily decaying dynamical system, since there is finite numerical dissipation. There is an abrupt, but almost imperceptible ($`<1\%`$) decrease in kinetic energy within the flow around $`t6`$, caused directly by formation of the Cat’s Eye. Beginning with Cat’s Eye disruption, however, there is a steady, steep decay of this quantity. By the end of this simulation the kinetic energy has dropped by about two orders of magnitude or more. This is in sharp contrast to the analogous 2D version of this K-H derived flow, where after formation of the Cat’s Eye the kinetic energy is virtually constant on these time scales. Turbulent decay in 3D is the reason for the difference, of course. Our result is, indeed, consistent with studies of 3D decaying MHD and HD turbulence (e.g., Mac Low et al. 1998; Stone et al. 1998; Porter et al. 1994), which also showed rapid dissipation of kinetic energy.

Three points are noticed from the kinetic energy plot of Fig. 3. First, in the medium resolution simulations, the decay rate increases with increasing initial magnetic field, from Case 1 to Case 4. This is the the result of enhanced dissipation through reconnection in MHD turbulence. Hence, even an obviously very weak magnetic field with $`M_{Ax}50`$ (or $`\beta 2.25\times 10^3`$ in our setup) does play an important role as an agency of enhancing dissipation. This character of increasing dissipation was also observed in 2D, although the dissipation there occurred through reconnection around the stable Cat’s Eye instead of reconnection driven by turbulent motion (Jones et al. 1997). Second, in Case 1 the decay is faster in the high resolution simulation than in the medium resolution calculation. This is because higher resolution allows a greater number of smaller scale structures to form and, so, reconnection events are more frequent. Again, a similar behavior was observed in the 2D VWF cases (Jones et al. 1997). Finally, the evolution curve of kinetic energy in Case 4 is quite different in the two simulations with different resolution. This happens because in the high resolution simulation, the magnetic field is amplified enough locally to play a more important dynamical role. The details are described in §3.3.

The magnetic energy plot of Fig. 3 shows the following behavior. Initially the magnetic energy increases at the expense of the kinetic energy, but stops increasing before energy equipartition is reached. After that, the magnetic energy starts decreasing, but the rate of decline is smaller than that of the kinetic energy. During this period, the flow character is close to that of HD turbulence. But eventually, progressing from smaller scales to larger scales (see the discussion on energy power spectrum below) the magnetic energy catches up the kinetic energy, and the character of MHD turbulence is fully established. Then, both energies decay with the same rate. The turbulence developed by the K-H instability in a closed system is a decaying, quasi-isotropic turbulence, and the symmetry of this flow does not support dynamo action. Hence, the magnetic energy must decay on some time scale along with the kinetic energy, and both energies should convert into the thermal energy. Over a very long time dependent on the effective magnetic Reynolds number at the dissipation scale (see below for more discussion), but much longer than our simulation, the magnetic field in this closed system should return to something resembling the initial configuration.

In simulations of ideal MHD flows, resistivity, $`\eta `$, is provided by numerical truncation and diffusion at the grid-cell level. So it does not have a constant value, but depends on the size of structures considered. In a numerical code based on a second-order scheme, such as the TVD scheme, the effective numerical resistivity is inversely proportional to the square of the scale, $`\mathrm{}`$, $`\eta \mathrm{}^2`$ (Ryu et al. 1995). As a result, the effective magnetic Reynolds number is proportional to the square of the scale, $`R_m\mathrm{}^2`$. We can use the evolution of magnetic energy in our simulations, in fact, to estimate heuristically the effective magnetic Reynolds numbers as follows. The magnetic energy decay rate for non-ideal decaying incompressible MHD turbulence (e.g., Biskamp 1993) is just

$$\frac{dE_m}{dt}=\eta j^2d^3x=\eta (\times B)^2d^3x2\eta \frac{E_m}{L_3^2}_,$$
(7)
where $`L_3`$ is the thickness of current sheet. From this we can write roughly that the magnetic field decay time is

$$t_{dm}\frac{L_3^2}{2\eta }\frac{R_mL_3}{2v}_,$$
(8)
where

$$\eta \frac{vL_{diss}}{R_m}$$
(9)
and $`L_{diss}L_3`$ are used since $`L_{diss}`$ represents the scale on which energy dissipation by reconnection occurs, that is the typical thickness of current sheets. Here, $`v`$ represents the typical flow velocity across the current sheet. For the Case 1h simulation using $`256^3`$ grid zones, for instance, we estimate from Fig.3 that $`t_{dm}20`$ and $`\sqrt{v^2}\sqrt{2E_k/\rho }5\times 10^2`$ at $`30t50`$. So for an effective magnetic Reynolds number corresponding to the typical scale of current sheet thickness, $`L_310^2L`$ (see Fig. 6 for an estimate of $`L_3`$), we obtain $`R_m200`$. Note that the smallest values for the scale $`L_3`$ correspond to $`23`$ grid zones, so they are numerically limited. The inertial range of turbulence in a simulation should require $`R_m10^3`$. Hence, applying the inverse square effective dissipation behavior of our second-order scheme, an inertial range is possible on scales greater than roughly $`8`$ zones. That is, in our simulations, turbulence can be approximately represented on scales larger than $`8`$ zones.

Additional insights about the evolution of fluid and magnetic field properties can be gleaned from the 3D power spectra of the kinetic end magnetic energies, $`E_k(k)`$ and $`E_m(k)`$, respectively, defined as follows. The Fourier amplitude of the kinetic energy is calculated as

$$A_{k,j}(\stackrel{}{k})=\frac{1}{L_xL_yL_z}H_xH_yH_z\sqrt{\rho }v_j\mathrm{exp}\left[i(k_xx+k_yy+k_zz)\right]d^3x,$$
(10)
and similarly the Fourier amplitude of the magnetic energy is calculated as

$$A_{m,j}(\stackrel{}{k})=\frac{1}{L_xL_yL_z}H_xH_yH_zB_j\mathrm{exp}\left[i(k_xx+k_yy+k_zz)\right]d^3x,$$
(11)
where $`j\{x,y,z\}`$. Here, $`H_x`$, $`H_y`$ and $`H_z`$ are the Hanning window functions (Press et al. 1986), which are given as

$$H_x=\frac{1}{2}\left[1\mathrm{cos}\left(2\pi \frac{x}{L_x}\right)\right],$$
(12)
and similarly for $`H_y`$ and $`H_z`$. Windowing in $`y`$ is used because the flow is not periodic in that direction. Then, it seems desirable to window in $`x`$ and $`z`$, too, in order avoid artificial anisotropies in Fourier space. Assuming isotropy on the scales of interest (see Fig. 5 and discussion in below), the energy power spectra are given as

$$E_k(k)=\frac{L_xL_yL_z}{2(2\pi )^3W}\left(|A_{k,x}(k)|^2+|A_{k,y}(k)|^2+|A_{k,z}(k)|^2\right)k^2,$$
(13)
and

$$E_m(k)=\frac{L_xL_yL_z}{2(2\pi )^3W}\left(|A_{m,x}(k)|^2+|A_{m,y}(k)|^2+|A_{m,z}(k)|^2\right)k^2,$$
(14)
where

$$W=\frac{1}{L_xL_yL_z}H_x^2H_y^2H_z^2d^3x.$$
(15)
Note that with the above definition

$$E_k(k)𝑑k=\frac{1}{L_xL_yL_zW}H_x^2H_y^2H_z^2\frac{1}{2}\rho v^2d^3x,$$
(16)
and

$$E_m(k)𝑑k=\frac{1}{L_xL_yL_zW}H_x^2H_y^2H_z^2\frac{1}{2}B^2d^3x,$$
(17)
that is, $`E_k(k)`$ and $`E_m(k)`$ are the kinetic and magnetic energies per unit $`k`$, respectively.

Fig. 4 shows the above two energy power spectra, along with their sum, $`E_{k+m}(k)E_k(k)+E_m(k)`$, for Case 1h at $`t=6,8,12,20,32`$ and $`50`$. Also, for comparison, we include lines with $`Ek^{5/3}`$ and $`Ek^3`$, representing the canonical forms for inertial range 3D and 2D isotropic HD turbulence, respectively (e.g., Lesieur 1997). The vertical dotted lines indicate the scales of $`L/4`$, one fourth of the box size, and 8 zones ($`\mathrm{log}k=0.602`$ and $`\mathrm{log}k=1.505`$). Structures with $`\mathrm{}L`$ (more specifically, $`\mathrm{}L/4`$ according to our tests) have been strongly affected by the finite box size, while structures with $`\mathrm{}8`$ zones have been severely dissipated by numerical diffusion (see above). So we may regard the region only between the two vertical dotted lines as an approximately inertial range. We have seen in Fig. 2 that at the earlier times, $`t=6`$ and $`8`$, there is still considerable large-scale, non-isotropic organization to the flows (e.g., the Cat’s Eye). But later, at $`t=20`$, the flow looks to the human eye as though it is isotropic turbulence. Fig. 4 supports that impression. In particular, if we examine $`E_{k+m}(k)`$ in the inertial range, we see that the power-law slope starts with a value close to $`3`$, as expected from the 2D flow character of the Cat’s Eye. Then, over time, the $`E_{k+m}(k)`$ becomes flatter, but the slope is still steeper than $`5/3`$, until the flow develops into something very close to decaying isotropic turbulence by $`t20`$. After that, the amplitude decays with time, but the form remains relatively unchanged to the end of the simulation.

Another point to emphasize is that at the early epochs $`E_k(k)`$ dominates $`E_m(k)`$ on all scales. This fact is consistent with our earlier conclusion that in the VWF cases, the flow is initially essentially HD in character. But, as complex flow structures develop, magnetic field is amplified by flux stretching. By $`t20`$, $`E_m(k)`$ has caught up $`E_k(k)`$ on small scales. By the end of the simulation, $`E_m(k)E_k(k)`$ over the most of inertial range except on the largest scales. Hence, by this time the flow of Case 1h shows the character of MHD turbulence.

The evolution to a quasi-isotropic flow character in the VWF case can be seen by looking at Fig. 5. This shows for Case 1h at a sequence of times $`v_x(y)_{x,z}`$, which is the average of $`v_x`$ over the $`xz`$ plane. $`v_x(y)_{x,z}𝑑y`$ is always very close to 0, from the symmetry of the initial conditions, although, since the initial perturbations were random, there is no exact symmetry in $`y`$ required. The shear in $`v_x`$, represented by $`dv_x(y)_{x,z}/dy`$, keeps decreasing rapidly. By $`t=40`$, not only $`dv_x(y)_{x,z}/dy0`$, but also $`v_x(y)_{x,z}0`$ for all $`y`$. This indicates there is no residual shear left, and the flow has become isotropic.

There are a number of quantitative ways to characterize the structural evolution of the flows in MHD simulations. The following quantities are particularly simple and useful: The mean magnetic curvature radius, $`L_1`$,

$$L_1\sqrt{\frac{B^4}{\left[\left(\stackrel{}{B}\stackrel{}{}\right)\stackrel{}{B}\right]^2}}_,$$
(18)
the flow Taylor microscale, $`L_2`$,

$$L_2\sqrt{\frac{v^2}{\left(\stackrel{}{}\times \stackrel{}{v}\right)^2}}_,$$
(19)
the magnetic Taylor microscale, $`L_3`$,

$$L_3\sqrt{\frac{B^2}{\left(\stackrel{}{}\times \stackrel{}{B}\right)^2}}_,$$
(20)
and the magnetic intermittency $`I`$,

$$I\frac{B^4}{B^2^2}_,$$
(21)
(e.g., Lesieur 1997 for $`L_2`$; Ethan T. Vishniac, private communication 1999 for others). The first of these, $`L_1`$, measures how sharply the magnetic field lines are bent. $`L_2`$ and $`L_3`$ measure the transverse dimensions, or thicknesses, of vortex tubes and current sheets, respectively. $`I`$ measures spatial contrast in the magnetic field strength distribution; i.e., $`I>>1`$ signifies the presence of magnetic voids and relatively intense flux tubes.

Fig. 6 shows the evolution of the above quantities for the three high resolution simulations Case 1h, 4h and 5h, which are VWF, VWF/WF (transitional) and WF cases, respectively (see the next two subsections for discussions on Cases 4h and 5h). Initially the magnetic field is uniform, so that $`L_1`$ and $`L_3`$ are infinite, while $`I=1`$. The initial shear layer gives $`L_2\sqrt{3La/4}0.17`$. We note that for the boundary conditions used the mean vector magnetic field, $`\stackrel{}{B}`$, is exactly constant; i.e., there is no dynamo action. So, any net increase in magnetic energy must also lead to the increase in the magnetic intermittency, $`I>1`$. The rest of this paragraph focuses on Case 1h, the VWF case. One can see from $`L_1`$, $`L_3`$ and $`I`$ that very quickly, on a time scale $`t12`$, the magnetic field is drawn into thin structures. Especially, reduction of $`L_1`$ signals formation of highly bent or twisted field regions. At the same time $`L_2`$ is reduced, because vortices of smaller scales form. The modest increase in $`L_2`$ just before $`t4`$ is due to the merger of smaller scale vortices, whose remnants are seen at $`t=6`$ in Fig. 2. By $`t68`$, when the Cat’s Eye is formed and begins to break up, magnetic intermittency, $`I`$, is already very large. Field curvature, $`L_1`$, stays low, but shows a peak around $`t8`$. This is because the field is wrapped into the Cat’s Eye. However, $`L_1`$ increases due to partial field relaxation during reconnection just prior to the Cat’s Eye break-up. During Cat’s Eye break-up the field becomes twisted and tangled, so that $`L_1`$ is again reduced. After the Cat’s Eye breaks up and the memory of the initial shear is gone in this case ($`t20`$), there begins a gradual relaxation in all of $`L_{1,2,3}`$; that is, the size of vortex tubes increases and the field becomes less strongly curved, while the thickness of current sheets increases and the magnetic intermittency decreases. The last of these remains relatively steady for Case 1h, near $`I2`$, for $`t30`$, while the others slowly increase to the end of that simulation. This behavior reflects the fact at late time that magnetic field is gradually relaxed by straightening itself, but flux tubes remain stable structures. These properties match the finding of the slow decay of magnetic energy for Case 1h in Fig. 3.

### 3.2 Weak Field (WF) Cases: Magnetic Reorganization of the Cat’s Eye

This subsection discusses Cases 5 and 6 ($`M_{Ax}=14.3,5`$), which in 2D were categorized as the WF or disruptive cases. The character of 3D flow and magnetic field evolution in these cases is perhaps best illustrated in the morphologies of Fig. 7. It shows the spatial distributions of magnetic field strength ($`B`$) and vorticity magnitude ($`\omega `$) in the high resolution simulation Case 5h at three epochs. Initially, the Cat’s Eye forms in this case, so the morphologies at $`t=6`$, although more sheet-like, carry some resemblance with those of Case 1h in Fig. 2. But, in Case 1h the analogous images at $`t=20`$ showed a completely disordered arrangement of magnetic and vortex tubes. Here those features are clearly laid out in patterns aligned to the original flow. Furthermore, a closer examination shows both the magnetic flux tubes and the vortex tubes to have a remarkable, sheet-like morphology with their minimum extent in the $`y`$ direction. Thus, they follow and resemble the original shear layer itself. There is also a good correspondence between strong vorticity regions and strong magnetic field regions, as one might expect in the presence of self-organization.

An interesting point for this simulation is that 2D cuts through fixed $`z`$’s resemble very much the 2D simulations of WF cases (see the images in Frank et al. 1996 and Jones et al. 1997). This is the result of the sheet-like morphology, and the indication that the flow and magnetic field evolution, although all three dimensions are available to the flows, the behavior is essentially 2D in character. So the WF cases in 3D evolve towards some degree of self-organized shear just as in 2D; that is, the Maxwell stresses developed during formation of the Cat’s Eye noticeably reorganize the flow and lead to significant alignments between magnetic and velocity fields. As for the VWF cases (Cases 1-3), the magnetic field itself becomes organized during Cat’s Eye development through the action of rib vortices into relatively strong-field flux tubes parallel to the original velocity field, separated from relatively weak fields trending along the Cat’s Eye axis. Subsequently, the reorganized velocity field aligns with the stronger magnetic field and retains a clear memory of the original velocity shear. Again, that contrasts with the vector fields in the weakest field VWF or HD cases in 3D, which become essentially isotropic in nature, excepting that the mean vector magnetic field must remain unchanged, due to the symmetry.

The above point is obvious in Fig. 8, which shows at a sequence of times $`v_x(y)_{x,z}`$ for the WF cases as well as the nonlinearly stable case (Case 7). The bottom panel in Fig. 8 shows at $`t=30`$ the shear strength, $`dv_x(y)_{x,z}/dy`$, in the original midplane of the shear layer. Two points are made from the figure. First, not only in the nonlinearly stable case but also in the WF cases, there left is still a well defined shear layer which is also laminar. This is the result of reorganization. Second, the residual shear strength at this time clearly scales with the initial magnetic field strength (or more importantly with $`B_{x0}`$). In HD, linear shear is stable against linear perturbations but unstable to 3D finite-amplitude perturbations (Bayly et al. 1988). But the magnetic field has the stabilizing effects, just as in the MHD K-H instability case. Stronger field can stabilize flows with larger linear shear. The linear correlation of the residual shear with the initial field strength is the direct consequence.

Among other things the self-organization and associated laminarity in the WF cases substantially slow the rate of kinetic energy dissipation, since it reduces the energy transfer to small, dissipation scales. That point is clearly made by comparing Fig. 3 and Fig. 9, which illustrate the evolution of energy partitioning for the VWF, VWF/WF (transitional) and WF cases. We make three points from Fig. 9. First, there is less kinetic energy dissipation in the stronger field Cases 6 than in Case 5, as expected from the above discussion on the residual shear. Second, in both WF Cases 5-6, about half of the initial kinetic energy is still present at $`t=20`$, and the decay rate has reduced significantly from what it was during the time of Cat’s Eye disruption. So, this flow pattern should continue for a moderately long time, but not as long as we found in 2D, since small scale structures in the third dimension can still form and enhance the dissipation. Third, a comparison of Cases 5h and 5m shows a good match between them. This indicates that small scale structures do not play a major role, although they do exist. At the same time, by this measure, we can state safely that the simulations are reasonably well resolved in the WF cases.

The overall behaviors of the mean magnetic curvature radius, $`L_1`$, the flow Taylor microscale, $`L_2`$, the magnetic Taylor microscale, $`L_3`$, and the magnetic intermittency $`I`$ in Case 5h are similar as those in Case 1h, as seen in Fig. 6. Three differences are noticed. First, the $`L`$’s remain small in Case 1h since there is little dynamical self-organization, But in Case 5h, self-organization relaxes the magnetic field as well as vortices. As a result, $`L`$’s increase after the Cat’s Eye starts to break apart. Second, the small peak in $`L_2`$ around $`t4`$ is missing in the WF case. This is because initially the formation of smaller scale vortices inside the Cat’s Eye is not allowed due to the magnetic field, although weak. This agrees with the visual impression that structures are absent inside the Cat’s Eye at $`t=6`$ in Fig. 7. Finally, $`I`$ approaches unity after $`t20`$, indicating the magnetic field has returned, more or less, to the initial uniform configuration.

### 3.3 Case 4: A Transitional VWF/WF Case with Eventual Reorganization

Case 4 with $`M_{Ax}=50`$ begins with a magnetic field too weak in 2D to have any immediate direct dynamical role, although through an accumulation of small magnetic field-induced perturbations even the 2D version of this case eventually begins to be distorted. Thus, in 2D we would have classified this as a dissipative case with a “footnote”. Here we will use the label VWF/WF. On the face of it, the 3D Case 4 looks during formation of the Cat’s Eye like the VWF cases discussed earlier, resembling the morphologies in Figs. 1 and 2 at $`t=6`$. There is even a briefly chaotic flow pattern right after the Cat’s Eye breaks up. However, in the high resolution simulation, Case 4h, slowly, over time, the flow begins to reorganize, so that by the end of the simulation residual shear becomes dominant, while the magnetic field has organized into one predominant flux tube parallel to the flow. This behavior is clearly visible in Fig. 10. So, effectively, this case behaves like the WF, disruptive cases discussed in the immediately preceding subsection. This case provides an evidence that the range of dynamically influential magnetic fields is greater in 3D than in 2D.

The causes of that difference can be seen clearly by a closer examination of magnetic field evolution during formation of the Cat’s Eye. The key is evident in Fig. 11, which at $`t=8`$ shows regions where the Alfvén Mach number is less than unity in Case 4h. This is just as the Cat’s Eye begins to fall apart. The regions with $`M_A<1`$ are all coincident with rib vortices that initially were HD in character. Now, however, they are magnetically dominated. Within the rib vortices at this time the flows are sub-Alfvénic, with mostly $`0.1M_A<1`$. An image showing the regions with small $`\beta =p_g/p_b`$ would be almost identical in appearance to Fig. 11, as well. The smallest values of $`\beta 1`$, so it is really the tension force rather than the pressure force that is revealing the magnetic field’s role.

In 2D we would have expected $`M_A`$ to drop by about one order of magnitude from its initial value, since the magnetic field lines around the vortex perimeter are stretched by about that much due to formation and rotation of the Cat’s Eye. That is insufficient to produce the properties seen in Fig. 11 and consistent with the observation of Jones et al. (1997) that a 2D $`M_{Ax}50`$ flow would not lead to magnetic dominance. In 3D, however, formation of the rib vortices provides a new mechanism to enhanced field amplification, as mentioned earlier. In particular, field lines become wrapped around the rib vortices, so that they become twisted as well as stretched around the Cat’s Eye. Fig. 11 shows this effect by tracing field lines within one rib vortex. Those field lines are clearly twisted around the structure, so that this feature is a legitimate flux tube with sufficient magnetic tension to begin a self-organization of the flow field. A close examination of the magnetic field distribution at this time reveals strengths about ten times amplified over the initial field as expected; namely, i.e., $`|B|0.2`$ along the entire length of each of the flux tubes visible in Fig. 11. But, each tube also contains a core down much of its length that has $`|B|0.4`$, an additional enhancement we attribute to twisting. The fraction of the flow under magnetic control at this time is still small, so the influence is not immediately obvious. It is, however, crucial to the eventual character of the flow. Evidently, if field amplification in any significant region is able to reduce the Alfvén Mach number to less than unity before the Cat’s Eye is HD disrupted, some memory of the original shear will be retained and, through self-organization, the magnetic and flow fields will align, and the flow may be smoothed. The evolution of these features for Cases 4h and 5h is clearly seen in animations of vorticity and magnetic pressure published in the CD ROM along with Jones et al. (1999), and currently posted at the web site given at the end of §1.

The ways in which the above physics impacts on energy evolution are shown for this case in Fig. 3. We see that in the high resolution simulation energy dissipation is intermediate between the quasi-HD Case 1h, which became turbulent, and Case 5h (in Fig. 9), which quickly developed into a smooth flow. We have noted previously that in medium and low resolution simulations Case 4 behaves as a VWF flow, since then numerical dissipation prohibits enough amplification of magnetic field to allow it to dominate dynamics. In addition, we can see that in Case 4h the rate of kinetic energy decay drops significantly after $`t25`$. By that time the flow has begun to organize strongly, and initially numerous magnetic flux tubes, twisted by vortical motion, have merged into a single, relatively intense structure. Note in this respect from Fig. 11 that the magnetic energy, $`E_m`$, in Fig. 3 is relatively constant from that time on, as well.

Fig. 12 shows images of the magnetic flux structures of Case 4h at $`t=40`$. The dominant flux tube is obvious. It contains most of the original flux that passed through the $`x=0`$ and $`x=L_x`$ faces of the computational box. Originally all the field lines ran obliquely in the $`xz`$ plane, but now most of the magnetic energy is concentrated in this one structure, aligned in the $`x`$ direction alone. On the other hand, for the periodic boundary conditions applied here the magnetic flux through each of the individual faces of the computational box is preserved. So, there must have been a topological change in the magnetic field along the way. That fact can also be seen in Fig. 12, where we see most of the flux through the $`z=0`$ and $`z=L_z`$ faces is now provided by field lines that meander around this dominant flux tube, and essentially orthogonal to it. That is, the magnetic flux has separated as a result of the instability and self-organization into two distinct domains.

Examination of the various structure measures in Fig. 6 augments the sense that Case 4h represents a transition between the quasi-HD Case 1 and the immediately reorganized Case 5. The magnetic field curvature measure, $`L_1`$, and the magnetic Taylor microscale, $`L_3`$, evolve in very similar ways for Cases 4 and 5, reflecting the fact that the magnetic field in each case is strong enough to smooth the flow and coalesce into one dynamically important flux tube. On the other hand, there is a much closer match between Cases 1 and 4 with regard to the flow Taylor microscale, $`L_2`$, reflecting the development of chaotic flow in each of these cases. The magnetic intermittency, $`I`$, stays moderately large, approaching $`23`$ by the end of the simulation, confirming the formation of one big flux tube.

## 4 SUMMARY AND DISCUSSION

Through high resolution MHD simulations using up to $`256^3`$ grid zones, we have studied the 3D nonlinear evolution of the compressible MHD K-H instability. As in our earlier work, we have considered periodic sections of flows that contain a thin shear layer, but are otherwise uniform. The initially uniform magnetic field is parallel to the shear plane, but $`30\mathrm{°}`$ oblique to the flow itself. Its strength spans the range corresponding to $`M_x=2.51.43\times 10^3`$. The sonic Mach number of the flow transition was initially unity.

The most important consequence of this work is confirmation in 3D flows of the conclusion from our earlier 2D work (Frank et al. 1995; Jones et al. 1997; Jeong et al. 2000) that even apparently weak magnetic fields corresponding to $`M_{Ax}4`$ can be important to nonlinear evolution of the K-H instability. The role of weak magnetic fields has been manifested in the following two ways. First, in the case of very weak field (VWF or dissipative case) with $`M_{Ax}50`$, dissipation is enhanced through magnetic reconnection. In this case, the instability remains essentially HD in character to the end. That is, the Cat’s Eye is destroyed by HD secondary instabilities and the flows are developed into mostly isotropic turbulence. But the decay rate of the turbulence increases. Second, in the case of weak field (WF or disruptive case) with $`4M_{Ax}50`$, the Cat’s Eye is destroyed by the magnetic stress of the field which has been amplified by stretching and twisting on the perimeter of the vortex, once the Alfvén Mach number of flows around the Cat’s Eye drops to unity or less. The flows, in this case, are eventually self-organized into relatively smooth ones with linear shear, which are stable against further instabilities.

There are two noticeable differences in the results of the current 3D work from that of our previous 2D work. First, in the VWF cases, while the Cat’s Eye remains stable in 2D, it is destroyed in 3D by inherently 3D instabilities (Hussain 1984; Bayly 1986; Craik & Criminale 1986). This is a HD process, which was fully studied before. The second difference, a 3D MHD process, is additional amplification of magnetic field by twisting inside rib vortices developed hydrodynamically around the Cat’s Eye. That leads to an increased role for magnetic stresses in destroying the Cat’s Eye. In both 2D and 3D WF cases magnetic field caught in the initially quasi-HD roll-up of the Cat’s Eye vortex tubes within the initial shear layer is amplified by field-line stretching. In 2D the field strength on the perimeter of the vortex is increased by about an order of magnitude, representing the increased length of field lines dragged around the forming Cat’s Eye, before they become subject to magnetic reconnection. In 3D that effect is further enhanced by the development of rib vortices spanning the Cat’s Eye that twist field lines into flux tubes, which then span the Cat’s Eye and apply a tension force to the plasma. In our 3D simulation, twisting of the field increases the maximum field strength by more than a factor of two over the 2D effect. For our rather idealized uniform density and vector field configurations, the boundary field strength for the WF case decreases to the value corresponding to $`M_{Ax}50`$ in 3D from $`M_{Ax}20`$ in 2D.

Two additional interesting points can also be made. First, in the WF cases, the magnetic energy reaches its maximum just as magnetic stresses begin to destroy the Cat’s Eye. However, it returns very close to its initial value, becoming almost uniform magnetic field again, on the time scale same as the flow becomes organized. The latter is barely longer than the time needed to form the Cat’s Eye and disrupt it. Hence, the field effectively plays the role of a catalyst. Second, in the VWF cases, where the flows become turbulent and, so, are not re-organized, the amplified magnetic field during the development of turbulence is retained long afterwards. The magnetic energy decays slowly, until it catches up to the kinetic energy. Then the flows approach decaying MHD turbulence, so that the magnetic energy along with the turbulent kinetic energy decays at an enhanced rate.

In the cases of strong field with $`M_{Ax}4`$, the development of the MHD K-H instability is essentially 2D in character, even though variation along the third dimension is allowed. When $`M_{Ax}2`$, the MHD flow is linearly stable and the instability is not initiated. When $`2M_{Ax}4`$, the shear layer is initially corrugated, but the enhanced magnetic tension due to the corrugated magnetic field stabilizes the instability before the Cat’s Eye forms; that is, the flow is nonlinearly stable.

The model configurations studied have no global helicity, and, thus, are not capable of dynamo action. Indeed the mean vector field is a constant throughout the simulations. So, the enhancement of magnetic energy comes from twisting followed by stretching of magnetic field lines, and/or from maintenance of significant magnetic intermittency (non-uniformity), not through generation of a large-scale field. We have seen that in the transitional VWF/WF case, however, the enhancement of magnetic energy is maintained through concentration of magnetic flux in flux tubes. Beyond the dynamical impact of such flux concentrations, that tendency could also be significant astrophysically for a different reason. In particular, measures of magnetic field strength, such as Faraday rotation, Zeeman splitting and the intensity of synchrotron emission become biased to those localized structures, so it becomes important to establish the intermittency of the field to understand associated observations.

We note that our simulations have been done in an idealized box with periodic boundaries along the $`x`$ and $`z`$-directions and reflecting boundaries along the $`y`$-direction. As a result, they have the following practical limitations. First, some astrophysical systems subject to the K-H instability, such as jets, contain continuous supplies of kinetic energy, while our simulations conserve the total energy. This limitation may be overcome by considering the convective K-H instability, which employs inflow/outflow boundaries along the $`x`$-direction, as Wu (1986) did in his 2D simulations. However, since one must then use a much larger computational domain to contain the evolving structures, the simulations become substantially more expensive, and it would be currently possible to follow only up to the early stage of the nonlinear development of the instability. A second constraint is due to the periodic nature of the $`z`$-boundaries, along the shear plane but perpendicular to the flow direction. From this symmetry the axis of the Cat’s Eye is constrained to be perpendicular to the initial flow direction. If the Cat’s Eye were allowed to rotate relative to the background flow, it would interact with the flow by crumpling or corrugating. This might lead somewhat different initial nonlinear behaviors in the WF and VWF/WF cases. Relaxation of that symmetry must await later work. Finally, the reflecting boundaries along the $`y`$-direction imposes another limitation of our simulations. However, in 2D, we saw that when activities are limited around the shear boundary as in the WF case, the effects of the reflecting boundaries are minimal (see, Frank et al. 1996; Malagoli et al. 1996). In addition, in the VWF case, flow develops into turbulence, so we expect the boundary effects would not be very important to the local properties of the flow.

Thus, we encourage restraint in direct application of our results for interpretation of observational features. Our intent is rather to provide more general physical insights into boundary layer dynamical processes within astrophysical objects subject to the K-H instability, such as jets associated with young stellar objects, accreting binaries, or larger scale flows from active galaxies (e.g., Ferrari et al. 1980), strongly sheared flows in the solar corona (e.g., Kopp 1992) and the earth’s magnetopause separating the magnetosphere from the solar wind (e.g., Miura 1984). In addition, our findings have several broader implications in astrophysics. The most obvious is that relatively weak magnetic fields may be able to reduce the development of turbulence from the K-H instability and diminish the tendency for mixing and related kinds of transport across slip surfaces. This work, thus, augments earlier suggestions that weak magnetic fields may inhibit turbulent diffusion (e.g., Vainshtein & Rosner 1991).

The work by DR was supported in part by KOSEF through grant and 981-0203-011-2. The work by TWJ was supported in part by the NSF through grants INT95-11654 and AST96-19438, by NASA grant NAG5-5055 and by the University of Minnesota Supercomputing Institute. The work by AF was supported by NSF Grant AST-0978765 and NASA Grant NAG5-8428. We thank the anonymous referee for clarifying comments.


format: text

text:
# 1 Introduction

## 1 Introduction

The subject of this paper is a problem that arises when a large body of data from diverse experiments is analyzed according to a theoretical model that has many adjustable parameters. Consider a generic data fitting problem based on experimental measurements $`\{D_I,I=1,\mathrm{},N\}`$ with errors $`\{\sigma _I\}`$. The data are to be compared to predictions $`\{T_I\}`$ from a theoretical model with unknown parameters $`\{a_i,i=1,\mathrm{},n\}`$. A common technique for comparing data with theory is to compute the $`\chi ^2`$ function defined by

$$\chi ^2=\underset{I=1}{\overset{N}{}}\left(\frac{D_IT_I}{\sigma _I}\right)^2,$$
(1)
or a generalization of that formula if correlations between the errors are known in terms of a set of correlation matrices. The physics objectives are (i) to find the best estimate of the parameters $`\{a_i\}`$ and their uncertainties, and (ii) to predict the values and uncertainties of physical quantities $`\{X^{(\alpha )},\alpha =1,2,\mathrm{}\}`$ that are functions of the $`\{a_i\}`$.

If the errors are randomly distributed, and the correlations well determined, then standard statistical methods of $`\chi ^2`$ minimization  apply, and established fitting tools like the CERN Library program MINUIT  can be employed. However, real problems are often more complex. This is particularly true in a “global analysis,” where the large number of data points $`\{D_I\}`$ do not come from a uniform set of measurements, but instead consist of a collection of results from many experiments on a variety of physical processes, with diverse characteristics and errors. The difficulties are compounded if there are unquantified theoretical uncertainties, if the number of theoretical parameters $`n`$ is large, or if the best parametrization cannot be uniquely defined a priori. All of these difficulties arise in the global analysis of hadronic parton distribution functions (PDFs) , which originally motivated this investigation. Several groups have addressed the question of estimating errors for the PDF determinations . But the problem is clearly more general than that application.

Of the many issues that confront a global analysis, we address in this paper two, for which we have been able to significantly improve on the traditional treatment. The improvements allow a more reliable determination of the uncertainties of $`\{a_i\}`$ and $`\{X^{(\alpha )}\}`$ in complex systems for which conventional methods may fail. To define these problems, we assume the system can be described by a global fitting function $`\chi _{\mathrm{global}}^2`$, or $`\chi ^2`$ for short, that characterizes the goodness-of-fit for a given set of theory parameters $`\{a_i\}`$. This $`\chi ^2`$ distills all available information on the theory and on the global experimental data sets, including their errors and correlations. One finds the minimum value $`\chi _0^2`$ of $`\chi ^2`$, and the best estimate of the theory parameters are the values $`\left\{a_i^0\right\}`$ that produce that minimum. The dependence of $`\chi ^2`$ on $`\{a_i\}`$ near the minimum provides information on the uncertainties in the $`\{a_i\}`$. These are usually characterized by the error matrix and its inverse, the Hessian matrix $`H_{ij}`$, where one assumes that $`\chi ^2`$ can be approximated by a quadratic expansion in $`\{a_i\}`$ around $`\left\{a_i^0\right\}`$. Once the Hessian is known, one can estimate not only the uncertainties of $`\{a_i\}`$, but also the uncertainty in the theoretical prediction for any physical quantity $`X`$, provided the dependence of $`X`$ on $`\{a_i\}`$ can be approximated by a linear expansion around $`\left\{a_i^0\right\}`$, and is thus characterized by its gradient at $`\left\{a_i^0\right\}`$ (cf. Sec. 2).

The first problem we address is a technical one that is important in practice: If the uncertainties are very disparate for different directions in the $`n`$-dimensional parameter space $`\{a_i\}`$, i.e., if the eigenvalues of $`H_{ij}`$ span many orders of magnitude, how can one calculate the matrix $`H_{ij}`$ with sufficient accuracy that reliable predictions are obtained for all directions? To solve this problem, we have developed an iterative procedure that adapts the step sizes used in the numerical calculation of the Hessian to the uncertainties in each eigenvector direction. We demonstrate the effectiveness of this procedure in our specific application, where the standard tool fails to yield reliable results.

The second problem we address concerns the reliability of estimating the uncertainty $`\mathrm{\Delta }X`$ in the prediction for some physical variable $`X`$ that is a function of the $`\{a_i\}`$: How can one estimate $`\mathrm{\Delta }X`$ in a way that takes into account the variation of $`\chi ^2`$ over the entire parameter space $`\{a_i\}`$, without assuming the quadratic approximation to $`\chi ^2`$ and the linear approximation to $`X`$ that are a part of the error matrix approach? We solve this problem by using Lagrange’s method of the undetermined multiplier to make constrained fits that derive the dependence of $`\chi ^2`$ on $`X`$. Because this method is more robust, it can be used by itself or to check the reliability of the Hessian method.

Section 2 summarizes the error matrix formalism and establishes our notation. Section 3 describes the iterative method for calculating the Hessian, and demonstrates its superiority in a concrete example. Section 4 introduces the Lagrange multiplier method and compares its results with the Hessian approach to the same application. Section 5 concludes.

## 2 Error Matrix and Hessian

First we review the well-known connection between the error matrix and the Hessian matrix of second derivatives. We emphasize the eigenvector representations of those matrices, which are used extensively later in the paper.

The basic assumption of the error matrix approach is that $`\chi ^2`$ can be approximated by a quadratic expansion in the fit parameters $`\{a_i\}`$ near the global minimum. This assumption will be true if the variation of the theory values $`T_I`$ with $`\{a_i\}`$ is approximately linear near the minimum. Defining $`y_i=a_ia_i^0`$ as the displacement of parameter $`a_i`$ from its value $`a_i^0`$ at the minimum, we have

$`\chi ^2`$ $`=`$ $`\chi _0^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \underset{i,j}{}}H_{ij}y_iy_j,`$ (2)
$`H_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{^2\chi ^2}{y_iy_j}}\right)_0,`$ (3)
where the derivatives are evaluated at the minimum point $`y_i=0`$ and $`H_{ij}`$ are the elements of the Hessian matrix.111We include a factor $`1/2`$ in the definition of $`H`$, as is the custom in high energy physics. There are no linear terms in $`y_i`$ in (2), because the first derivatives of $`\chi ^2`$ are zero at the minimum.

Being a symmetric matrix, $`H_{ij}`$ has a complete set of $`n`$ orthonormal eigenvectors $`V_i^{(k)}v_{ik}`$ with eigenvalues $`ϵ_k`$:

$`{\displaystyle \underset{j}{}}H_{ij}v_{jk}`$ $`=`$ $`ϵ_kv_{ik}`$ (4)
$`{\displaystyle \underset{i}{}}v_{ij}v_{ik}`$ $`=`$ $`\delta _{jk}.`$ (5)
These eigenvectors provide a natural basis to express arbitrary variations around the minimum: we replace $`\{y_i\}`$ by a new set of parameters $`\{z_i\}`$ defined by

$$y_i=\underset{j}{}v_{ij}\sqrt{\frac{1}{ϵ_j}}z_j.$$
(6)
These parameters have the simple property that

$$\mathrm{\Delta }\chi ^2=\chi ^2\chi _0^{\mathrm{\hspace{0.17em}2}}=\underset{i}{}z_i^{\mathrm{\hspace{0.17em}2}}.$$
(7)
In other words, the surfaces of constant $`\chi ^2`$ are spheres in $`\{z_i\}`$ space, with $`\mathrm{\Delta }\chi ^2`$ the squared distance from the minimum.

The orthonormality of $`v_{ij}`$ can be used to invert the transformation (6):

$`z_i=\sqrt{ϵ_i}{\displaystyle \underset{j}{}}y_jv_{ji}.`$ (8)
The Hessian and its inverse, which is the error matrix, are easily expressed in terms of the eigenvalues and eigenvector components:

$`H_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}ϵ_kv_{ik}v_{jk}`$ (9)
$`(H^1)_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{1}{ϵ_k}}v_{ik}v_{jk}.`$ (10)
Now consider any physical quantity $`X`$ that can be calculated according to the theory as a function of the parameters $`\{a_i\}`$. The best estimate of $`X`$ is the value at the minimum $`X_0=X(a_i^0)`$. In the neighborhood of the minimum, assuming the first term of the Taylor-series expansion of $`X`$ gives an adequate approximation, the deviation of $`X`$ from its best estimate is given by

$$\mathrm{\Delta }X=XX_0\underset{i}{}\frac{X}{y_i}y_i=\underset{i}{}X_iz_i$$
(11)
where

$$X_i\frac{X}{z_i}$$
(12)
are the components of the $`z`$-gradient evaluated at the global minimum, i.e., at the origin in $`z`$-space.

Since $`\chi ^2`$ increases uniformly in all directions in $`z`$-space, the gradient vector $`X_i`$ gives the direction in which the physical observable $`X`$ varies fastest with increasing $`\chi ^2`$. The maximum deviation in $`X`$ for a given increase in $`\chi ^2`$ is therefore obtained by the dot product of the gradient vector $`X_i`$ and a displacement vector $`Z_i`$ in the same direction with length $`\sqrt{\mathrm{\Delta }\chi ^2}`$, i.e., $`Z_i=X_i\sqrt{\mathrm{\Delta }\chi ^2/\mathrm{\Sigma }_jX_j^2}`$. For the square of the deviation, we therefore obtain the simpler formula

$$(\mathrm{\Delta }X)^2=(XZ)^2=\mathrm{\Delta }\chi ^2\underset{i}{}X_i^{\mathrm{\hspace{0.17em}2}}.$$
(13)
The traditional formula for the error estimate $`(\mathrm{\Delta }X)^2`$ in terms of the original coordinates $`\{y_i\}`$ can be derived by substituting $`X_i=\frac{X}{z_i}={\displaystyle \underset{j}{}}\frac{X}{y_j}\frac{y_j}{z_i}`$ in (13) and using (6) and (10). The result is

$$(\mathrm{\Delta }X)^2=\mathrm{\Delta }\chi ^2\underset{i,j}{}\frac{X}{y_i}(H^1)_{ij}\frac{X}{y_j}.$$
(14)
This standard result can of course also be derived directly by minimizing $`\chi ^2`$ in (2) with respect to $`\{a_i\}`$, subject to a constraint on $`X`$.

Equations (13) and (14) are equivalent if the assumptions of a linear approximation for $`X`$ and a quadratic approximation for $`\chi ^2`$ are exact. But in practice, the numerical accuracy of the two can differ considerably if these conditions are not well met over the relevant region of parameter space. To calculate the error estimate $`\mathrm{\Delta }X`$, we prefer to use Eq. (13) using derivatives $`X_i`$ calculated by finite differences of $`X`$ at the points $`z_i=\pm \frac{1}{2}\sqrt{\mathrm{\Delta }\chi ^2}`$ (with $`z_j=0`$ for $`ji`$). This is generally more accurate, because it estimates the necessary derivatives using an appropriate step size, and thus reduces the effect of higher order terms and numerical noise.

In a complex problem such as a global analysis, the region of applicability of the approximations is generally unknown beforehand. A situation of particular concern is when the various eigenvalues $`\{ϵ_i\}`$ have very different orders of magnitude—signaling that the function $`\chi ^2`$ varies slowly in some directions of $`a_i`$ space, and rapidly in others. The iterative method described in the next section is designed to deal effectively with this situation.

## 3 Iterative Procedure

In practical applications, the Hessian matrix $`H_{ij}`$ is calculated using finite differences to estimate the second derivatives in (3). A balance must be maintained in choosing the step sizes for this, since higher-order terms will contribute if the intervals are too large, while numerical noise will dominate if the intervals are too small. This noise problem may arise more often than is generally realized, since the theory values $`\{T_I\}`$ that enter the $`\chi ^2`$ calculation may not be the ideally smooth functions of the fit parameters that one would associate with analytic formulas. For in complex theoretical models, the $`\{T_I\}`$ may be computed from multiple integrals that have small discontinuities as functions of $`\{a_i\}`$ induced by adaptive integration methods. These numerical errors forbid the use of a very small step size in the finite difference calculations of derivatives. Furthermore, as noted above, the eigenvalues of $`H_{ij}`$ may span a wide range, so excellent accuracy is needed especially to get the smaller ones right.

### The Procedure

We want to evaluate $`H_{ij}`$ by sampling the values of $`\chi ^2`$ in a region of parameter space where Eq. (2) is a good approximation. In principle, the parameters $`\{z_i\}`$ are the natural choice for exploring this space; but of course they are not known in advance. We therefore adopt the following iterative procedure:

1. Define a new set of coordinates $`\{\xi _i\}`$ by

$`y_i={\displaystyle \underset{j}{}}u_{ij}t_j\xi _j`$ (15)
where $`u_{ij}`$ is an orthogonal matrix and $`\{t_i\}`$ are scale factors. In the first iteration, these are chosen as $`u_{ij}=\delta _{ij}`$ and $`t_i=1`$, so that $`\xi _i=y_i`$. This makes the first round of iteration similar to the usual procedure of taking derivatives with respect to $`a_i`$. The iterative method is designed such that with successive iterations, $`u_{ij}`$, $`t_i`$, and $`\xi _i`$ converge to $`v_{ij}`$, $`\sqrt{1/ϵ_i}`$, and $`z_i`$ respectively.
2. Calculate the effective second derivative matrix $`\mathrm{\Phi }_{ij}`$ defined by

$`\chi ^2`$ $`=`$ $`\chi _0^{\mathrm{\hspace{0.17em}2}}+{\displaystyle \underset{ij}{}}\mathrm{\Phi }_{ij}\xi _i\xi _j`$ (16)
$`\mathrm{\Phi }_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2\chi ^2}{\xi _i\xi _j}}`$ (17)
using finite differences of the $`\xi _i`$. The step size in $`\xi _i`$ is chosen to make the increase in $`\chi ^2`$ due to the diagonal element $`\mathrm{\Phi }_{ii}\xi _i^2`$ equal to a certain value $`\delta \chi ^2`$. The choice of $`\delta \chi ^2`$ is determined by the particular physics application at hand. Naively, one might expect $`\delta \chi ^21`$ to be the right choice. That would indeed be appropriate for a $`\chi ^2`$ function obeying ideal statistical requirements. But when the input to $`\chi _{\mathrm{global}}^2`$ is imperfect, a reasonable choice of $`\delta \chi ^2`$ must be based on a physics judgement of the appropriate range of that particular $`\chi ^2`$ function. We therefore leave the choice of $`\delta \chi ^2`$ open in this general discussion.222Cf. discussion in the following subsection on a sample problem. In any case, if the final results are to be trustworthy, they must not be sensitive to that choice.

We calculate each off-diagonal second derivative by evaluating $`\chi ^2`$ at the four corners of the rectangle $`(+\delta _i,+\delta _j)`$, $`(\delta _i,\delta _j)`$, $`(+\delta _i,\delta _j)`$, $`(\delta _i,+\delta _j)`$, where $`\delta _i`$ is the step size. This is a modification of the technique used in MINUIT . For the sake of efficiency, the MINUIT subroutine HESSE estimates off-diagonal elements using only one of those corners, together with values at $`(\delta _i,0)`$ and $`(0,\delta _j)`$ that are already known from calculating the diagonal elements of the Hessian. Our method is slower by a factor of 4, but is more accurate because it fully or partly cancels some of the contributions from higher derivatives. The first derivatives $`\chi ^2/\xi _i`$ are also calculated at this stage of the iteration and used to refine the estimate of the location of the minimum.
3. Compute the Hessian according to $`\mathrm{\Phi }_{ij}`$,

$$H_{ij}=\underset{m,n}{}\frac{\mathrm{\Phi }_{mn}u_{im}u_{jn}}{t_mt_n}.$$
(18)
4. Find the normalized eigenvectors of the Hessian, as defined by Eqs. (4) and (5).
5. Replace $`u_{ij}`$ by $`v_{ij}`$, $`t_j`$ by $`\sqrt{1/ϵ_j}`$, and go back to step 1. The steps are repeated typically 10–20 times, until the changes become small and $`\mathrm{\Phi }_{ij}`$ converges to $`\delta _{ij}`$.

This iterative procedure improves the estimate of the Hessian matrix, and hence of the error matrix, because in the later iterations it calculates the Hessian based on points that sample the region where $`\mathrm{\Delta }\chi ^2`$ has the magnitude of physical interest.

### Results from a Sample Application

As an example, we apply the iterative procedure to the application that motivated this study—the global analysis of PDFs —and compare the results with those obtained from MINUIT. The experimental input for this problem consists of $`N=1295`$ data points from 15 different experimental data sets involving four distinct physical processes. All the potential complexities mentioned earlier are present in this system. The theory is the quark parton model, based on next-to-leading order perturbative Quantum Chromodynamics (QCD). The model contains $`n=16`$ parameters $`a_i`$ that characterize the quark and gluon distributions in the proton at some low momentum scale $`Q_0`$. From a calculational point of view, the theoretical model consists of the numerical integration of an integro-differential equation and multiple convolution integrals that are evaluated mostly by adaptive algorithms. The fitting function $`\chi _{\mathrm{global}}^2`$ in this case combines the published statistical and systematic errors of the data points in quadrature. The only correlated errors incorporated are the published overall normalization uncertainties of the individual experiments. The fitting program is the same as that used to generate the CTEQ parton distributions . The global $`\chi ^2`$ minimum for this system defines the CTEQ5M1 set of PDFs, for which $`\chi _0^21200`$ . We find that the eigenvalues $`\{ϵ_i\}`$ of the Hessian for this system range over 5–6 orders of magnitude (distributed approximately exponentially).

The value of $`\mathrm{\Delta }\chi ^2`$ that corresponds to a given confidence level is well defined for an ideal experiment: e.g., $`\mathrm{\Delta }\chi ^2<1`$ defines the $`68\%`$ confidence region. But in a real-world global analysis, the experimental and theoretical values in Eq. (1) include systematic errors, and the uncertainties $`\sigma _I`$ include subjective estimates of those errors, so the relation between $`\mathrm{\Delta }\chi ^2`$ and confidence level requires further analysis. From independent detailed studies of the uncertainties , we estimate that an appropriate choice of $`\delta \chi ^2`$ for the iterative calculation is around $`10`$ in our application, and only the region $`\mathrm{\Delta }\chi ^2>100`$ can be ruled out for the final fits.

The error matrix approach relies on a quadratic approximation to $`\chi ^2`$ in the neighborhood of the minimum. To test whether that approximation is valid, we plot $`\chi ^2`$ as a function of distance along a particular direction in $`\{a_i\}`$ space, as shown in Fig. 1. The direction chosen is a typical one—specifically it is the direction of the eigenvector with median eigenvalue. The dotted curve in Fig. 1 is the exact $`\chi ^2`$ and the solid curve is the quadratic approximation (2). The approximation is seen to provide a rather good description of the function. Even at points where $`\chi ^2`$ has increased by $`50`$, the quadratic approximation reproduces the increase to $`20`$% accuracy.

To correctly measure the curvature of the quadratic approximation, it is important to fit points that are displaced by an appropriate distance. This can be seen from Fig. 2, which shows the difference between the two curves in Fig. 1 in the central region. The difference displays a small cubic contribution to $`\chi ^2`$. It also reveals contributions that vary erratically with a magnitude on the order of $`0.03`$. These fluctuations come from the noise associated with switching of intervals in the adaptive integration routines. Because the fluctuations are small, they do not affect our results in principle. But they do require care in estimating the derivatives. In particular, they would make finite-difference estimates based on small intervals extremely unreliable. The iterative method avoids this problem by choosing a suitable scale for each eigenvector direction when evaluating the Hessian.

Figures 1 and 2 show the behavior of $`\chi ^2`$ along a single typical direction in the 16 dimensional parameter space. Fig. 3 shows a complementary test of the iterative method for all possible directions. We have chosen 1000 directions at random in $`\{z_i\}`$ space. We displace the parameters away from the minimum in each of these directions by a distance that makes $`\mathrm{\Delta }\chi ^2=5`$. We then compute the value of $`\mathrm{\Delta }\chi ^2`$ predicted by the quadratic approximation (2), using the Hessian calculated by the iterative method and, for comparison, by the routine HESSE within the MINUIT package. The results are displayed in Fig. 3 as histograms, with $`\mathrm{\Delta }\chi ^2`$ on the horizontal axis and the number of counts on the vertical axis. If $`\chi ^2`$ were quadratic in $`\{a_i\}`$, then a perfect computational method would yield a delta function at $`\mathrm{\Delta }\chi ^2=5`$. Fig. 3 shows that:

* For the solid histogram—the result of the iterative procedure—the quadratic approximation is close to the exact result in all directions, and hence Eq. (2) is a pretty good representation of $`\chi ^2`$. Quantitatively, the middle $`68`$% of the distribution is contained in the region $`5.4\pm 0.6`$.
* For the dotted histogram—based on the general purpose program MINUIT—the distribution is also spread around the expected value of $`5`$, but it is very broadly distributed. This estimate of the Hessian is therefore unsatisfactory, because we might be interested in a quantity whose gradient direction is one for which the Hessian computed by MINUIT is widely off the mark. A major source of this problem is the numerical noise visible in Fig. 2: MINUIT uses a small step size to calculate the derivatives, and gets misled by the small-scale discontinuities in $`\chi ^2`$. For some directions, $`\mathrm{\Delta }\chi ^2`$ even becomes negative because the errors in one or more of the small eigenvalues are big enough to allow their calculated values to become negative. (Within MINUIT, this circumstance elicits a warning message, and a constant is added to all the eigenvalues, which in the context of Fig. 3 corresponds to shifting the dotted distribution to the right.)

Figure 4 shows the results of a similar study, in which the 1000 random directions are chosen only from the subspace of $`\{z_i\}`$ that is spanned by the 10 directions with the largest eigenvalues $`ϵ_i`$. The larger eigenvalues correspond to directions in which $`\chi ^2`$ rises most rapidly, or in other words, directions in which the parameters are more strongly constrained by data. Because the distance moved away from the minimum in $`\{a_i\}`$ space is smaller in this case, the quadratic approximation is generally better, so it is not surprising that the histograms are more sharply peaked than in Fig. 3. But the advantage of the iterative method remains apparent.

### Comment

Information from the iteratively-improved Hessian provides a useful tool for refining the choice of functional forms used to parametrize a continuous degree of freedom in the theoretical model. Explicitly, the relevant formulas are as follows.

The length squared of the displacement vector in the space of fit parameters is

$`{\displaystyle \underset{i}{}}(a_ia_i^0)^2={\displaystyle \underset{i}{}}y_i^{\mathrm{\hspace{0.17em}2}}={\displaystyle \underset{i}{}}{\displaystyle \frac{z_i^{\mathrm{\hspace{0.17em}2}}}{ϵ_i}}`$ (19)
while $`\mathrm{\Delta }\chi ^2=_iz_i^{\mathrm{\hspace{0.17em}2}}`$ by (7). Hence the directions in which the parameters are well determined (the steep directions) correspond to eigenvectors of the Hessian with large eigenvalues, while the shallow directions in which they are weakly determined correspond to small eigenvalues.

The extreme values for any particular $`a_i`$ are

$`a_i=a_i^0\pm \mathrm{\Delta }a_i`$ (20)
where

$`(\mathrm{\Delta }a_i)^2=\mathrm{\Delta }\chi ^2{\displaystyle \underset{j}{}}{\displaystyle \frac{v_{ij}^{\mathrm{\hspace{0.17em}2}}}{ϵ_j}}.`$ (21)
Equation (21) can be used to see if each parameter is appropriately well constrained. Furthermore, the individual terms in the sum show the contributions to $`\mathrm{\Delta }a_i`$ from the various eigenvectors, so if a parametrization leads to a poorly defined minimum because it allows too much freedom—which is indicated by a failure of the iteration to converge for the smallest eigenvalues of the Hessian—it is easy to see which of the parameters are most responsible for the too-shallow directions.

## 4 Lagrange Multiplier Method

The Hessian, via its inverse which is the error matrix, provides a general way to propagate the uncertainties of experimental and theoretical input to the fit parameters $`\{a_i\}`$, and thence on to a measurable quantity $`X(\{a_i\})`$ by Eqs. (13) or (14). But these equations are based on assuming that $`\chi ^2`$ and $`X`$ can be treated as quadratic and linear functions of $`\{a_i\}`$ respectively. In this section we describe a different approach, based on the mathematical method of the Lagrange undetermined multiplier, which avoids those assumptions.

### The Procedure

Let $`X_0`$ be the value of $`X`$ at the $`\chi ^2`$ minimum, which is the best estimate of $`X`$. For a fixed value of $`\lambda `$, called the Lagrange multiplier, one performs a new minimization with respect to the fit parameters $`\{a_i\}`$, this time on the quantity

$$F=\chi ^2+\lambda (XX_0),$$
(22)
to obtain a pair of values $`(\chi ^2(\lambda ),X(\lambda ))`$. (The constant term $`\lambda X_0`$ here is not necessary, because it does not affect the minimization; but it makes the minimum value of $`F`$ easier to interpret.) At this new minimum, $`\chi ^2(\lambda )`$ is the lowest possible $`\chi ^2`$ for the corresponding value $`X(\lambda )`$ of the physical variable $`X`$. Thus one achieves a constrained fit in which $`\chi ^2`$ is minimized for a particular value of $`X`$.

By repeating the minimization for many values of $`\lambda `$, one maps out the parametrically-defined curve $`(\chi ^2(\lambda ),X(\lambda ))`$. Since $`\lambda `$ is just the parameter for this curve, its value is of no particular physical significance. The relevant range for $`\lambda `$ can be found by trial-and-error; or it can be estimated using the Hessian approximation, which predicts that $`\lambda 2\mathrm{\Delta }\chi ^2/\mathrm{\Delta }X`$. In that approximation, $`F`$ goes down by the same amount that $`\chi ^2`$ goes up.

One way to understand the Lagrange Multiplier method is to imagine that the quantity $`X`$ is simply one of the fitting parameters, say $`a_1`$. The variation of $`\chi ^2`$ with $`a_1`$ could be mapped out by minimizing $`\chi ^2`$ with respect to $`\{a_2,\mathrm{},a_n\}`$ for a sequence of values of $`a_1`$. (That operation is indeed so useful that MINUIT provides a procedure MINOS to carry it out.) In the more general case that $`X`$ is a function of all $`\{a_i\}`$, one wants to similarly minimize $`\chi ^2`$ for fixed values of $`X`$. That is exactly what the Lagrange Multiplier method does, since including the undetermined multiplier term in (22) renders the $`\{a_i\}`$ independent in spite of the constraint on $`X`$.

A more phenomenological way to understand the Lagrange Multiplier method is to imagine that $`X`$ has just been measured, with result $`X_{\mathrm{new}}\pm \sigma _{\mathrm{new}}`$. To decide whether this hypothetical new measurement is consistent with the old body of data, one would add a term $`[(X_{\mathrm{new}}X)/\sigma _{\mathrm{new}}]^2`$ to $`\chi _{\mathrm{global}}^2`$ of Eq. (1) and redo the minimization. The added contribution to $`\chi ^2`$ consists of a constant, a linear term in $`X`$, and a quadratic term in $`X`$. This is equivalent to Eq. (22), because a constraint on $`X^2`$ is equivalent to a constraint on $`X`$ itself.

The essential feature of the Lagrange Multiplier method is that, for a given $`\mathrm{\Delta }\chi ^2`$, it finds the largest range of $`X`$ allowed by the global data set and the theoretical model, independent of any approximations. The full parameter space $`\{a_i\}`$ is explored in the minimization procedure, not just the immediate neighborhood of the original $`\chi ^2`$ minimum as in the Hessian method, and no approximations based on a small deviation from the original minimum are needed.

The only drawback to the Lagrange Multiplier method is that it can be slow computationally, since it requires a separate series of minimizations for each observable $`X`$ that is of interest.

### Example

We now look at an example of the Lagrange Multiplier method from our application, the uncertainty of parton distribution functions. For the physical quantity $`X`$, we consider the cross section $`\sigma _W`$ for $`W^\pm `$ production in $`p\overline{p}`$ collisions at the energy $`\sqrt{s}=1.8`$ TeV of the Tevatron collider at Fermilab. We want to estimate $`\sigma _W`$ and the uncertainty on that estimate, based on the global analysis of parton distributions.

The points in Fig. 5 show $`\chi _{\mathrm{global}}^2`$ as a function of $`\sigma _WB`$ in nanobarns, where $`B=0.106`$ is the branching ratio assumed for $`We\nu `$. These points are obtained by the Lagrange Multiplier method using $`\lambda =0`$, $`\pm 1000`$, $`\pm 2000`$, $`\pm 3000`$, $`\pm 4000`$. They are thus discrete cases of $`\chi _{\mathrm{global}}^2`$ versus $`\sigma _W`$, without approximations.

The smooth curve in Fig. 5 is the parabola given in Eq. (14), using the Hessian computed by the iterative method and treating $`\sigma _W`$ in the linear approximation. The comparison between this curve and the discrete points from the Lagrange Multiplier calculation tests the quality of the quadratic and linear approximations and the reliability of the iterative calculation of $`H_{ij}`$. For this application, we conclude that the improved Hessian method works very well, since the difference between the points and the curve is small, and indicates only a small cubic contribution. If the two results did not agree, the correct result would be the one given by the Lagrange Multiplier method.

To estimate $`\mathrm{\Delta }X`$, the uncertainty of $`X`$ consistent with the global analysis of existing data, one needs to specify what range of $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$ is allowed. As discussed earlier, the acceptable limit of $`\chi _{\mathrm{global}}^2`$ depends on the nature of the original definition of this fitting function, in the context of the specific system under study. For the case of $`\sigma _W`$, Ref.  estimates $`\mathrm{\Delta }\chi _{\mathrm{global}}^2100`$, which translates into $`\mathrm{\Delta }\sigma _W/\sigma _W\pm 3`$% according to Fig. 5.

## 5 Conclusion

We have addressed some computational problems that arise in a global phenomenological analysis, in which a complex theory with many parameters confronts a large number of data points from diverse experiments.

The traditional error-matrix analysis is based on a quadratic approximation to the function $`\chi ^2`$ that measures the quality of the fit, in the neighborhood of the minimum that defines the best fit. The iterative method proposed in Sec. 3 improves the calculation of the Hessian matrix which expresses that quadratic approximation, for a complex system in which general-purpose programs may fall short. The inverse of this improved version of the Hessian matrix is an improved version of the error matrix. It can be used to estimate the uncertainty of predictions using standard error matrix formulas.

Our iterative procedure for calculating the Hessian is implemented as an extension to the widely-used CERN fortran library routine MINUIT . The code is available from http://www.pa.msu.edu/$``$pumplin/iterate/, or it can be requested by e-mail from
pumplin@pa.msu.edu. Included with this code is a test example which demonstrates that the iterative method is superior to standard MINUIT even for a $`\chi ^2`$ function that has no numerical noise of the type encountered in Fig. 2.

The Lagrange Multiplier method proposed in Sec. 4 calculates the uncertainty on a given physical observable directly, without going through the error matrix. It thus avoids the assumption that the theoretical quantities can be approximated by linear functions of the search parameters, which is intrinsic to the Hessian approach.

For simplicity, we have discussed only the problem of obtaining error estimates on a single quantity $`X`$. It is straightforward to generalize our methods to find the region allowed simultaneously for two or more variables by a given $`\mathrm{\Delta }\chi ^2`$. For example, in the case of two variables $`X^{(1)}`$ and $`X^{(2)}`$, the allowed region according to the Hessian method is the interior of an ellipse. The Lagrange multiplier method can be generalized for this case by adding two terms, $`\lambda _1X^{(1)}+\lambda _2X^{(2)}`$, to $`\chi ^2`$.

Although the Lagrange Multiplier procedure is conceptually simple and straightforward to implement, it is slow computationally because it requires many full minimizations to map out $`\chi ^2`$ as a function of $`X`$, and this must be done separately for each quantity $`X`$ whose error limits are of interest. In contrast, once the Hessian has been determined from the global analysis, it can be applied to any physical observable. One needs only to compute the gradient $`X/a_i`$ of the observable $`X`$ and substitute into Eq. (14); or better, to compute the gradient $`X_i=X/z_i`$ and substitute into Eq. (13). For computational efficiency, the iteratively calculated Hessian is therefore the method of choice, provided its linear approximations are sufficiently accurate. Whether or not that is the case can be determined by comparing Hessian and Lagrange Multiplier results. We use both methods in a detailed study of the uncertainties in the CTEQ5 parton distribution functions  that is based on the work presented here.

## Acknowledgments

We thank R. Brock, D. Casey, J. Huston, and F. Olness for discussions on the uncertainties of parton distribution functions. We thank M. Botje for comments on an earlier version of the manuscript. This work was supported in part by NSF grant PHY-9802564.