Additive Self Helicity as a Kink Mode Threshold

TL;DR

This paper introduces additive self helicity as a potential threshold indicator for kink instability in complex magnetic flux tubes, extending the concept of twist helicity beyond thin flux tube approximations.

Contribution

It proposes that additive self helicity can serve as a generalized twist measure for complex equilibria and provides a numerical method for its calculation.

Findings

Additive self helicity increases with twist in flux tubes.

Kink instability coincides with a specific threshold of additive self helicity.

H_A/Φ^2 acts as a generalized twist number for complex flux tubes.

Abstract

In this paper we propose that additive self helicity, introduced by Longcope and Malanushenko (2008), plays a role in the kink instability for complex equilibria, similar to twist helicity for thin flux tubes (Hood and Priest (1979), Berger and Field (1984)). We support this hypothesis by a calculation of additive self helicity of a twisted flux tube from the simulation of Fan and Gibson (2003). As more twist gets introduced, the additive self helicity increases, and the kink instability of the tube coincides with the drop of additive self helicity, after the latter reaches the value of $H_A/\Phi^2\approx 1.5$ (where $\Phi$ is the flux of the tube and $H_A$ is additive self helicity). We compare additive self helicity to twist for a thin sub-portion of the tube to illustrate that $H_A/\Phi^2$ is equal to the twist number, studied by Berger and Field (1984), when the thin flux tube approximation is applicable. We suggest, that the quantity $H_A/\Phi^2$ could be treated as a generalization of a twist number, when thin flux tube approximation is not applicable. A threshold on a generalized twist number might prove extremely useful studying complex equilibria, just as twist number itself has proven useful studying idealized thin flux tubes. We explicitly describe a numerical method for calculating additive self helicity, which includes an algorithm for identifying a domain occupied by a flux bundle and a method of calculating potential magnetic field confined to this domain. We also describe a numerical method to calculate twist of a thin flux tube, using a frame parallelly transported along the axis of the tube.

Figures

Equations

$$H=H_T+H_W.$$

$$L=Tw+Wr.$$

$$T{w}_{(gen)}=H_A/{\Phi }^2.$$

$$H_A\left(\mathbf{B},\mathbf{P}\left(\mathcal{D}\right),\mathcal{D}\right)=\underset{\mathcal{D}}{\int }\left(\mathbf{B}-\mathbf{P}\right)\cdot \left(\mathbf{A}+{\mathbf{A}}_P\right)dV,$$

$${\mathbf{P}\cdot \hat{\mathbf{n}}|}_{\partial \mathcal{D}}=\mathbf{B}\cdot \hat{\mathbf{n}}{|}_{\partial \mathcal{D}},$$

$$\Phi =\underset{z=0,B_z\ge 0}{\int }{B}_{z}dxdy=-\underset{z=0,B_z\le 0}{\int }{B}_{z}dxdy.$$

$$\Theta (\mathbf{r})=\{\begin{array}{rl}1& \text{mbox},if\mathbf{r}\in \mathcal{D}\\ 0& \text{mbox},if\mathbf{r}\notin \mathcal{D}\text{mbox}.\end{array}$$

$$\nabla \cdot {\mathbf{B}}_P={\nabla }^2\chi =0.$$

$$\begin{array}[]{rcl}\left.\frac{\partial\chi}{\partial n}\right|_{\partial{\cal D},z\neq 0}&=&0,\\ \left.\frac{\partial\chi}{\partial z}\right|_{\partial{\cal D},z=0}&=&B_{z}(x,y).\end{array}$$

$$\begin{array}[]{lcr}B_{x}(x_{i},y_{j},z_{k})&=&\cfrac{\chi(x_{i+1/2},y_{j},z_{k})-\chi(x_{i-1/2},y_{j},z_{k})}{x_{i+1/2}-x_{i-1/2}},\\ \end{array},$$

$${\nabla }^2\chi \left(x_i,y_j\right)\approx \frac{1}{h^2}\left(\chi \left(x_{i-1},y_j\right)+\chi \left(x_{i+1},y_j\right)+\chi \left(x_i,y_{j-1}\right)+\chi \left(x_i,y_{j+1}\right)-4\chi \left(x_i,y_j\right)\right)$$

$$\chi \left(x_i,y_j\right)\approx \frac{1}{4}\left(\chi \left(x_{i-1},y_j\right)+\chi \left(x_{i+1},y_j\right)+\chi \left(x_i,y_{j-1}\right)+\chi \left(x_i,y_{j+1}\right)\right)-\frac{h^2}{4}{\nabla }^2\chi \left(x_i,y_j\right).$$

$${\chi }^{[n+1]}\left(x_i,y_j\right)=\frac{1}{4}\left({\chi }^{[n]}\left(x_{i-1},y_j\right)+{\chi }^{[n]}\left(x_{i+1},y_j\right)+{\chi }^{[n]}\left(x_i,y_{j-1}\right)+{\chi }^{[n]}\left(x_i,y_{j+1}\right)\right).$$

$$6{\chi }^{[n+1]}(O)={\chi }^{[n]}(A_1)+{\chi }^{[n]}(A_2)+{\chi }^{[n]}(B_1)+{\chi }^{[n]}(B_2)+{\chi }^{[n]}(C_1)+{\chi }^{[n]}(C_2).$$

$$\begin{array}[]{rcl}\chi^{[n]}(A_{1})&=&\frac{1}{4}\left[\chi^{[n]}(SA_{1})+\chi^{[n]}(TA_{1})+\chi^{[n]}(SA_{1})+\chi^{[n]}(O)\right],\\ \chi^{[n]}(B_{1})&=&\frac{1}{4}\left[\chi^{[n]}(SB_{1})+\chi^{[n]}(TB_{1})+\chi^{[n]}(SB_{1})+\chi^{[n]}(O)\right],\\ \chi^{[n]}(C_{1})&=&\frac{1}{4}\left[\chi^{[n]}(TA_{1})+\chi^{[n]}(TA_{2})+\chi^{[n]}(TB_{1})+\chi^{[n]}(TB_{2})\right]\end{array}$$

$$\begin{array}[]{rcl}6\chi^{[n+1]}(O)&=&2\times\frac{1}{4}\left[\chi^{[n]}(TA_{1})+\chi^{[n]}(TA_{2})+\chi^{[n]}(BA_{1})+\chi^{[n]}(BA_{2})\right]+\\ &+&2\times\frac{1}{4}\left[\chi^{[n]}(TB_{1})+\chi^{[n]}(TB_{2})+\chi^{[n]}(BB_{1})+\chi^{[n]}(BB_{2})\right]+\\ &+&\frac{1}{4}\left[\chi^{[n]}(SA_{1})+\chi^{[n]}(SA_{2})+\chi^{[n]}(SB_{1})+\chi^{[n]}(SB_{2})\right]+\\ &+&4\times\frac{1}{4}\chi^{[n]}(O),\end{array}$$

$$\begin{array}[]{rcl}\chi^{[n+1]}(O)&=&\frac{1}{12}\left[\chi^{[n]}(TA_{1})+\chi^{[n]}(TA_{2})+\chi^{[n]}(BA_{1})+\chi^{[n]}(BA_{2})\right]+\\ &+&\frac{1}{12}\left[\chi^{[n]}(TB_{1})+\chi^{[n]}(TB_{2})+\chi^{[n]}(BB_{1})+\chi^{[n]}(BB_{2})\right]+\\ &+&\frac{1}{24}\left[\chi^{[n]}(SA_{1})+\chi^{[n]}(SA_{2})+\chi^{[n]}(SB_{1})+\chi^{[n]}(SB_{2})\right]+\\ &+&\frac{1}{6}\chi^{[n]}(O).\end{array}$$

$${\mathbf{B}}_0=B_{\psi }\hat{\psi }+B_{\phi }\hat{\phi }=B_t{e}^{-{\varpi }^2/a^2}\left(q\frac{\varpi }{\rho }\hat{\psi }+\frac{a}{\rho }\hat{\phi }\right),$$

$$\begin{array}[]{rcl}A_{x}(x,y,z)&=&\int\limits_{0}^{z}B_{y}(x,y,z^{\prime})dz^{\prime}\\ A_{y}(x,y,z)&=&f(x,y)-\int\limits_{0}^{z}B_{x}(x,y,z^{\prime})dz^{\prime}\\ f(x,y)&=&\int\limits_{0}^{x}B_{z}(x^{\prime},y,0)dx^{\prime}\end{array}$$

$$H_A\left(t,{\varpi }_{max}\right)=\underset{\Theta \left(t,{\varpi }_{max}\right)}{\int }\left(\Theta \mathbf{B}-\mathbf{P}\right)\cdot \left(A+A_P\right)dV$$

$${\hat{\mathbf{u}}}_1=\frac{{\hat{\mathbf{u}}}_0-\left({\hat{\mathbf{u}}}_0\cdot {\text{lˆ}}_1\right){\text{lˆ}}_1}{\left|{\hat{\mathbf{u}}}_0-\left({\hat{\mathbf{u}}}_0\cdot {\text{lˆ}}_1\right){\text{lˆ}}_1\right|},$$

$$H(\mathbf{B},{\mathbf{P}}_{\mathcal{V}},\mathcal{V})\equiv \underset{\mathcal{V}}{\int }{d}^{3}x\mathbf{B}\cdot \mathbf{A}-\underset{\mathcal{V}}{\int }{d}^{3}x{\mathbf{P}}_{\mathcal{V}}\cdot {\mathbf{A}}_P+\underset{z=0}{\int }dxdy{B}_z(x,y,0)\underset{{\mathbf{x}}_0}{\overset{\mathbf{x}}{\int }}d{\mathbf{x}}^{\prime }\left[\mathbf{A}({\mathbf{x}}^{\prime })-{\mathbf{A}}_P({\mathbf{x}}^{\prime })\right]$$

$$H\left({\Theta }_{\mathcal{D}}\mathbf{B},{\mathbf{P}}_{\mathcal{D}},\mathcal{V}\right)+H\left({\mathbf{P}}_{\mathcal{D}},{\mathbf{P}}_{\mathcal{V},\Theta },\mathcal{V}\right)=H\left({\Theta }_{\mathcal{D}}\mathbf{B},{\mathbf{P}}_{\mathcal{V},\Theta },\mathcal{V}\right),$$

Full text

We compare additive self helicity to twist for a thin sub-portion of the tube to illustrate that $H_{A}/\Phi^{2}$ is equal to the twist number, studied by Berger & Field (1984), when the thin flux tube approximation is applicable. We suggest, that the quantity $H_{A}/\Phi^{2}$ could be treated as a generalization of a twist number, when thin flux tube approximation is not applicable. A threshold on a generalized twist number might prove extremely useful studying complex equilibria, just as twist number itself has proven useful studying idealized thin flux tubes. We explicitly describe a numerical method for calculating additive self helicity, which includes an algorithm for identifying a domain occupied by a flux bundle and a method of calculating potential magnetic field confined to this domain. We also describe a numerical method to calculate twist of a thin flux tube, using a frame parallelly transported along the axis of the tube.

1 Introduction

According to a prevalent model coronal mass ejections (CMEs) are triggered by current-driven magnetohydrodynamic (MHD) instability related to the external kink mode (Hood & Priest, 1979; Török et al., 2004; Rachmeler et al., 2009). The external kink mode, in its strictest form, is a helical deformation of an initially symmetric, cylindrical equilibrium, consisting of helically twisted field lines. The equilibrium is unstable to this instability if its field lines twist about the axis by more than a critical angle, typically close to $3\pi$ radians (Hood & Priest, 1979; Baty, 2001). The helical deformation leads to an overall decrease in magnetic energy, since it shortens many field lines even as it lengthens the axis.

Equilibria without symmetry can undergo an analogous form of current-driven instability under which global motion lowers the magnetic energy (Bernstein et al., 1958; Newcomb, 1960). Such an instability implies the existence of another equilibrium with lower magnetic energy. The spontaneous motion tends to deform the unstable field into a state resembling the lower energy equilibrium. Indeed, it is generally expected that there is at least one minimum energy state from which deformation cannot lower the the magnetic energy without breaking magnetic field lines; its energy is the absolute minimum under ideal motion.

Linear stability and instability are determined by the energy change under infinitesimal motions. An equilibrium will change energy only at the second order since first order changes vanish as a requirement for force balance. Ideal stability demands that no deformation decrease the energy at second order, while instability will result if even one energy-decreasing motion is possible. The infinite variety of possible motions make it impractical to establish stability in any but the simplest and most symmetric equilibria.

Based on analogy to axisymmetric systems it is expected that general equilibria, including those relevant to CMEs, are probably unstable when some portion of their field lines are twisted about one another by more than some critical angle. This expectation was mentioned in a study by Fan & Gibson (2003) of the evolution of a toroidal flux rope into a pre-existing coronal arcade. They solved time-dependent equations of MHD in a three-dimensional, rectangular domain. Flux tube emergence was simulated by kinematically introducing an isolated toroidal field through the lower boundary. The toroidal field was introduced beneath a pre-existing arcade slowly enough that the coronal response never approached the local Alfvén speed. Fan and Gibson concluded that the system underwent a current-driven instability after a critical amount of the torus had been introduced. They bolstered this claim by performing an auxiliary run where the kinematic emergence was halted and the system allowed to evolve freely; it settled into an equilibrium.

While twist angle has proven useful in a few cases, it is difficult to demonstrate its utility as a threshold in general, asymmetric equilibria. Indeed, in any but a few very symmetric cases there is no simple, obvious way to define the angle by which the field lines wrap about one another. The local rate of twist is given by the current density, which is after all the source of free energy powering the instability. On the other hand, excessive local current density is not sufficient to drive instability. This fact is illustrated by numerous examples of discontinuous field which are minimum energy states.

It has been suggested that a threshold exists, in general equilibria, for some global quantity such as free magnetic energy or helicity (Zhang et al., 2006; Low, 1994). If this is the case then we expect the instability to lower the value of this global quantity so that it falls below the threshold value in the lower-energy, stable equilibrium. Magnetic helicity is a logical candidate to play this role since it is proportional to total twist angle in cylindrical fields. Relative helicity in particular is a proxy for currents. Helicity is, however, conserved under ideal motion and therefore will not be reduced to a sub-threshold value by an ideal instability.

The total helicity of a thin, isolated flux tube can be written as a sum of two terms called twist and writhe (Berger & Field, 1984; Moffatt & Ricca, 1992),

[TABLE]

The writhe depends on the configuration of the tube’s axis while the twist depends on the wrapping of field lines about one another. A cylindrical tube has a perfectly straight axis and therefore zero writhe helicity. Any ideal motion which helically deforms the entire flux tube will increase the magnitude of the writhe helicity. Since the motion preserves total helicity the change in writhe must be accompanied by an offsetting change in twist helicity. If the writhe has the same sign as the initial twist, then the motion will decrease the twist helicity. In cases where the magnetic energy depends mostly on twist, this motion will decrease the magnetic energy (Linton & Antiochos, 2002). The straight equilibrium is therefore unstable to an external kink mode.

Topologically, the foregoing properties of magnetic field lines could be compared to the properties of thin closed ribbons. One may introduce twist number, writhe number and their combination, called linkage number, is a preserved quantity in the absence of reconnection (Berger & Field, 1984; Moffatt & Ricca, 1992),

[TABLE]

By analogy to the case of a thin isolated flux tube we consider the twist helicity, rather than the total helicity, to be the most likely candidate for a stability threshold. Indeed, within a thin flux tube it is possible to derive a net twist angle among field lines and $H_{T}=\Phi^{2}Tw=\Phi^{2}\Delta\theta/2\pi$, where $\Phi$ is the total magnetic flux through a cross-section of the tube and $\Delta\theta$ is the net twist angle.

Twist and writhe are, however, defined only in cases of thin, isolated magnetic flux tubes, and can no more set the threshold we seek than the net twist angle can.

Recently Longcope & Malanushenko (2008) introduced two generalizations of relative helicity applicable to arbitrary sub-volumes of a magnetic field. They termed both generalized self-helicity, and the two differed only by the reference field used in their computation. The one called additive self-helicity (that we denote $H_{A}$) uses a reference field confined to the same sub-volume as the original field, and can be interpreted as a generalization of the twist helicity to arbitrary magnetic fields. The additive self-helicity of a thin, isolated flux tube is exactly the twist helicity.

Since the additive self-helicity can be computed for arbitrary magnetic fields we propose that it (normalized by the squared flux) is the quantity to which current-driven instability sets an upper limit, which could be considered a generalized twist number:

[TABLE]

The paper is organized as follows. In Section 2, we describe a method for calculating additive self helicity and $Tw_{(gen)}$ numerically. There are two large and nontrivial parts of this calculation, that we describe in 2.1 and 2.2: locating a domain containing a given flux bundle and constructing a potential field in this domain by Jacobi relaxation. In Section 3, we apply the method to a simulation to support our hypothesis, the emerging twisted flux tube from Fan & Gibson (2003). In 3.1 we briefly describe this simulation, and then in 3.2 we show different embedded domains defined by different subportions of the footpoints. In 3.2 we describe, how the twist of Berger & Field (1984) could be calculated for those of the domains for which thin flux tube approximation is applicable. In Section 4 we present the evolution of additive self helicity, unconfined self-helicity, twist (for “thin” domains) and the integrated helicity flux in the simulation. We demonstrate that $Tw_{(gen)}$ increases corresponding to helicity flux, that it drops after it reaches a certain value (about $1.5$) and that this drop coincides with the rapid expansion of the tube due to the kink instability. We also demonstrate that the unconfined self helicity grows only when helicity flux is nonzero and that it stays constant when kink instability happens. We also show that $Tw_{(gen)}$ corresponds to $Tw$ when thin flux tube approximation is applicable.

2 Numerical Solutions

The object of study is a magnetic field ${\bf B}\left({\bf r}\right)$ defined in a domain ${\cal D}$, ${\bf r}\in{\cal D}$, that lies on and above the photosphere, $z\geq 0$. By domain we understand a volume that encloses the field: ${\bf B}\cdot{\bf\hat{n}}=0$ on all boundaries, $\partial{\cal D}$, except at the photosphere, where ${\bf B}\cdot{\bf\hat{n}}=B_{z}(x,y,z=0)$. An example of such a volume is the coronal part of an $\Omega$-shaped loop. The self-helicity is given by

[TABLE]

as defined in Longcope & Malanushenko (2008). Here ${\bf P}$ is the potential magnetic field, whose normal component matches the normal component of ${\bf B}$ on the boundary $\partial{\cal D}$,

[TABLE]

${\bf A}$ and ${\bf A}_{P}$ are the vector potentials of ${\bf B}$ and ${\bf P}$ respectively (as discussed in Finn & Antonsen (1985), helicity, defined this way is gauge-independent).

Once the self-helicity is known, the twist is given by eq. (1) with $\Phi$ being the total signed flux of the footpoints of the configuration:

[TABLE]

In the next two sections we discuss methods of numerically obtaining ${\cal D}$, from given footpoints, and ${\bf P}$.

2.1 Finding the domain.

In order to describe the domain on a grid we introduce the support function:

[TABLE]

This is a function of the given magnetic field ${\bf B}$ and some photospheric area, called the boundary mask. By definition, every field line, initiated at any point on the boundary mask and having the other footpoint somewhere within the mask, is completely inside the domain ${\cal D}$. If the field line traced in both directions from some coronal point ends within the photospheric mask, then this point also belongs to the domain. In numerical computations we replace “point” with a small finite volume, voxel $v_{ijk}$ (3-dimensional pixel). We define a voxel to be inside ${\cal D}$ (equivalent to saying $\Theta\left(v_{ijk}\right)=1$), if there is at least one point inside it that belongs to ${\cal D}$.

The simplest method of constructing the support function would be to trace a field line in both direction from every voxel of the computational grid, set $\Theta=1$ in the voxel if the footpoints both terminate in pixels from the boundary regions, and set $\Theta=0$ otherwise. This, however, is a very time-consuming algorithm, especially for a large arrays of data. Instead we use an algorithm which reduces the computational time by tracing field lines from a subset of voxels. It works by progressively adding voxels to $\Theta$ adjacent to those already known to belong to ${\cal D}$.

We add a voxel centered at $r_{i,j,k}$ to the domain under two different circumstances. 1. A field line initialized somewhere within the volume of the voxel $v_{i,j,k}$, centered at $r_{i,j,k}$, is found to have both footpoints within the boundary mask. 2. A field line initiated in some other voxel, and determined to belong to ${\cal D}$, passes through some portion of the volume $v_{i,j,k}$.

Initially, the domain consists only of footpoint voxels, so the initial step is to trace field lines initiated at the footpoints, assuming, that at least some of these lines will lie in the domain.

We illustrate the method on a simplistic case of a potential magnetic field, confined to a half-space, with $B_{z}=0$ everywhere at the photosphere, except at four pixels, as shown on Fig. 2. We have computed the magnetic field inside a small box of $15\times 15\times 15$ pixels, centered around the photospheric sources. The boundary mask consists of these four voxels at the photosphere with non-zero vertical magnetic field. In this simplistic example the initial guess would be four field lines, initiated at four footpoint voxels, as shown on Fig. 2, left (note that in this particular example a field line, initiated at one voxel, ends at another voxel within the mask and thus is the same as the field line, initiated at that another voxel, so these four initial guesses are really two, not four field lines). The voxels of the initial guess are shown with crosses.

In an algorithm, this would be the first step:

• Step 1:

Make the initial guess: trace field lines from the footpoints.

As soon as an initial step is made, the next step is to assume, that the immediate neighbourhood of voxels known to be in ${\cal D}$ are likely to be also in the domain. Thus, in the next (iterative) search the following steps are performed:

• Step 2:

Locate voxels on the boundary of the current domain.

• Step 3:

For every voxel on the boundary: trace a field line and check whether it is in the domain.

• If yes:

Add the voxel to the domain. Add all voxels along the line to the domain. Exclude them from the boundary (there is no need to check them again).

• If no:

Mark the voxel as “questionable”. (If there is a field line, which passes through the voxel and does not belong to the domain, then at least part of the voxel is outside of the domain. Since its immediate neighbourhood is in the domain, then it is possible that part of it is also in the domain.)

• Loop:

Repeat steps 2-3 until all the voxels in the boundary are “questionable” and no new voxels are added.

When the iterative search does not find any new voxels, we make the final check of the boundary voxels. The idea is to trace field lines from all corners of such “questionable” voxels to see, which corners (and thus which part of a voxel) belongs to the domain. We consider this to be optional check, which may improve the precision of the definition of the domain by at most one layer of voxels.

This last search may also give information about the normal to the domain surface. If it is known that some corners of a voxel are in the domain and some are not, it is possible to approximate the boundary as a plane separating those two groups of corners.

• Step 4, optional:

For each voxel, marked previously as “questionable”, check the corners (by tracing field lines) to see which of them are in the domain and which aren’t. Keep this information.

2.2 Constructing the Confined Potential Field ${\bf P}$

Once the domain has been determined, the next step is to construct the potential magnetic field confined to it. We use a common relaxation method on a staggered grid in orded to account for the complex boundaries of ${\cal D}$.

We introduce a scalar potential ${\bf B}_{p}=\nabla\chi$ and look for the solution of the Laplace’s equation for $\chi$

[TABLE]

By the definition of ${\cal D}$, field lines never cross $\partial{\cal D}$, except at the lower boundary, $z=0$. Thus, boundary conditions for ${\bf B}_{p}$ could be written as: $\left.{\bf B}_{p}\cdot{\bf\hat{n}}\right|_{\partial{\cal D},z\neq 0}=0$ and $\left.{\bf B}_{p}\cdot{\bf\hat{n}}\right|_{\partial{\cal D},z=0}={\bf B}_{z}(x,y)$. This is equivalent to Neumann boundary conditions for $\chi$:

[TABLE]

The Algorithm for the Relaxation Method

We use the Jacobi iterative method (see, for example LeVeque, 1955) to solve for the potential field. Here we briefly summarize the algorithm and further explain in details. The $n+1$-th iteration is

$\forall{\bf r}\in{\cal D}$: calculate a new iteration $\chi^{[n+1]}$ as a solution of the equation $\chi^{[n+1]}-\chi^{[n]}=Kh^{2}\nabla^{2}\chi^{[n]}$, where $h$ is the grid spacing. The Laplacian $\nabla^{2}\chi^{[n]}\left({\bf r}\right)$, found ausing standard finite difference methods, is equivalent to an average over some stencil of neighbouring points minus the central value; $K$ is a constant that depends on the exact shape of the stencil. 2. 2.

$\forall{\bf r}_{b}\in\partial{\cal D}$: set $\chi^{[n+1]}\left({\bf r}_{b}\right)$ so as to satisfy boundary conditions (BCs). 3. 3.

Repeat steps 1–2, until the difference between $\chi^{[n]}\left({\bf r}\right)$ and $\chi^{[n+1]}\left({\bf r}\right)$ is sufficiently small in some sense (namely, until $||\chi^{[n+1]}-\chi^{[n]}||<\epsilon$, where $\epsilon$ is pre-defined small number).

**Staggered Mesh

**

The functions $B_{x}(x,y,z)$, $B_{y}(x,y,z)$ and $B_{z}(x,y,z)$ are defined on the same mesh points $(x_{i},y_{j},z_{k})$. If we are interested in finding $\chi(x,y,z)$, so that $B_{x}=\frac{\partial\chi}{\partial x}$, $B_{y}=\frac{\partial\chi}{\partial y}$ and $B_{z}=\frac{\partial\chi}{\partial z}$, it advantageous to define $\chi$ in between the original mesh points and calculate the derivatives using finite difference as following:

[TABLE]

and so on for $B_{y}$ and $B_{z}$. $\chi$, then, would only be defined in the middle of the faces of cubic voxels, i.e., at points $(i\pm 1/2,j,k)$, $(i,j\pm 1/2)$ and $(i,j,k\pm 1/2)$.

Such a mesh, called a “cartesian staggered mesh”, is known to have better numerical properties, such as immunity from decoupling of variables and having a smaller numeric dispersion (Perot, 2000, see, for example).

The finite difference approximation of a Laplacian at one point can be interpreted as a weighted average over a stencil of several points minus the value at that point. For example, in the 2D case the second order approximation to $\nabla^{2}\chi(x,y)$ on a uniform Cartesian grid at the point $\left(x_{i},y_{j}\right)$ could be computed over a 5-point stencil:

[TABLE]

(here $h$ is the spacing of the grid). It could be rewritten as

[TABLE]

The Jacobi method uses this equation to iteratively update the value at the point, constantly assuming $\nabla^{2}\chi=0$. In the case of the 5-points stencil the updated value would be

[TABLE]

In our case of a 3D staggered mesh, choosing a stencil becomes more complicated. We propose a 13-point scheme, shown on the right of Fig. 3 (black dots). To motivate this stencil, we derive it from the “unstaggered” one (Fig. 3, left, gray dots). In an “unstaggered” finite differencing scheme the $[n+1]$-th iteration in Jacobi method would be expressed as

[TABLE]

But for the staggered mesh $\chi$ is undefined at these nodes. This can be resolved by setting $\chi$ at each “gray” point to be equal to the average of its 4 closest neighbours,

[TABLE]

and so on. Then we may substitute this in the original expression and get:

[TABLE]

which is eqivalent to

[TABLE]

With these weights th “farthest” nodes $S[AB]_{[12]}$ have half the influence on the laplacian, of the “closer” nodes. Note also, that the sum of the weights is one.

Boundary Conditions

Boundary conditions (given by Eqn. 5) in the staggered mesh is particularly easy if one assumes that the boundary surface passes inside of boundary voxels, rather than on their sides. Suppose, for example, that the boundary plane normal to ${\bf\hat{z}}$ passes through the center of the voxel $v_{ijk}$. Then the BC for this voxel would be that $B_{z}\left(i,j,k\right)=0$, or simply $\chi\left(i,j,k+\frac{1}{2}\right)=\chi\left(i,j,k-\frac{1}{2}\right)$.

To motivate such choice of the boundary, we note that boundary voxels, by definition, are the voxels part of which is inside of ${\cal D}$ while part is outside. Such a conclusion is made about voxels, some of whose corners are inside of ${\cal D}$, and some of the corners are outside of ${\cal D}$ (this information about the domain is obtained in the step 4 of the algorithm, described in section 2.1). We approximate the boundary inside of each boundary voxel as a plane, that passes through the center of the voxel and that separates its “exterior” part from its “interior” part. Such approximation will err by no more that $1/\sqrt{2}$ voxel’s length off the real location of the boundary. We also find it easier to work in terms of faces rather than corners, since this is where $\chi$ is defined. (We say, that a face is “exterior” to the domain if more than two of its corners are not in the domain, i.e., for a voxel, we say, that if only one corner or only one edge are “exterior”, we do not consider it a subject to BC’s).

There are several ways to orient such a boundary plane inside a voxel, based on the behaviour of the boundary in the immediate surrounding of the voxel.

The voxel has only one face outside of the domain. Then we consider the boundary parallel to that face of the voxel (see Fig. 4, left). If, say, the boundary is parallel to the face between faces $A$ and $A_{1}$ (see Fig. 4, bottom left), then the normal field to the boundary is ${\bf B}\cdot\bf\widehat{AA}_{1}$ (hereafter $\bf\widehat{AA}_{1}$ denotes a unit vector along the line from $A$ to $A_{1}$, which might be $\pm{\bf\hat{x}}$, $\pm{\bf\hat{y}}$ or $\pm{\bf\hat{z}}$), and BC would be formulated as

\begin{array}[]{lcccccc}\chi_{A}&=&1\cdot\chi_{A^{\prime}}&+&0\cdot\chi_{B^{\prime}}&+&0\cdot\chi_{C^{\prime}}\end{array}. 2. 2.

The voxel has two adjacent faces outside of the domain. Then we approximate the boundary as a plane, that cuts off these two faces, as shown on Fig. 4, middle. If faces $A$ and $B$ are outside and faces $A_{1}$ and $B_{1}$ are inside of the domain, then we consider the normal field to be ${\bf B}\cdot\frac{1}{\sqrt{2}}\left(\bf{\widehat{AA}_{1}}+\bf{\widehat{BB}_{1}}\right)$ and set BC’s as

\begin{array}[]{lcccccc}\chi_{A}&=&0\cdot\chi_{A^{\prime}}&+&1\cdot\chi_{B^{\prime}}&+&0\cdot\chi_{C^{\prime}},\\ \chi_{B}&=&1\cdot\chi_{A^{\prime}}&+&0\cdot\chi_{B^{\prime}}&+&0\cdot\chi_{C^{\prime}}.\\ \end{array} 3. 3.

Similarly, if three mutually adjacent faces of the voxel are outside of the domain (and three others are inside), as shown on Fig. 4, right, then, analogously, we assume that the normal field is ${\bf B}\cdot\frac{1}{\sqrt{3}}\left(\bf{\widehat{AA}_{1}}+\bf{\widehat{BB}_{1}}+\bf{\widehat{CC}_{1}}\right)$ and BC’s could be set in the following way:

\begin{array}[]{lcccccc}\chi_{A}&=&0\cdot\chi_{A^{\prime}}&+&\frac{1}{2}\cdot\chi_{B^{\prime}}&+&\frac{1}{2}\cdot\chi_{C^{\prime}},\\ \chi_{B}&=&\frac{1}{2}\cdot\chi_{A^{\prime}}&+&0\cdot\chi_{B^{\prime}}&+&\frac{1}{2}\cdot\chi_{C^{\prime}},\\ \chi_{C}&=&\frac{1}{2}\cdot\chi_{A^{\prime}}&+&\frac{1}{2}\cdot\chi_{B^{\prime}}&+&0\cdot\chi_{C^{\prime}}.\\ \end{array}

(Note that in this case there are really three variables and one equation to satisfy; thus, there are different solutions to $\chi$. But each of those solutions would be valid, as long as it satisfies ${\bf B}\cdot{\bf\hat{n}}=0$.) 4. 4.

“Everything else”: the voxel has three or more non-adjacent faces that are outside of the domain, but still is on the boundary. It is considered an extraneous voxel and is removed from the boundary.

3 The Experiment

The method described above was tested on a simple quadrupole example, and the values of self-helicity it gives are in a good agreement with theoretical predictions (Longcope & Malanushenko, 2008). That work, however, does not consider any sort of stable equilibrium and does not study any kinking instability thresholds, similar to those developed in Hood & Priest, 1981.

The objective of the current work is to test whether the parameter $H_{A}/\Phi^{2}$ behaves like a total twist in the sense that it has a critical value above which a system is unstable to a global disruption. To do so, we use the numerical simulation of kink instability in an emerging flux tube from Fan & Gibson, 2003 .

3.1 Simulation Data

The initial configuration is a linear arcade above the photosphere, into which a thick, non-force-free torus was emerged. Inside the torus the field lines wind around its minor axis and the field magnitude drops with distance from the minor axis. The exact shape of the magnetic field, in the coordinates shown on Fig. 5, is the following:

[TABLE]

where $a=0.1L$ is the minor radius, $R=0.375L$ is the major radius, $q=-1$, $B_{t}=9B_{0}$, $L$ is the length scale of the domain (further in our calculations $L=1$), $B_{0}$ is the characteristic strength of the photospheric arcade the torus is emerging into, and the time is given in the units of Alfven time, $\tau_{A}=L/v_{A}$. The field strength drops as $e^{-\varpi^{2}/a^{2}}$ with $\varpi$ being the distance from the minor axis. At $\varpi=3a$ magnetic field was artificially set to 0.

The torus is “emerged” from underneath the photosphere with a constant speed. There is a mass flow across the photosphere in the area, and the emerging tube is driven into the domain by an electric field at the boundary. This is made in the following way: for each time step (starting at t=0 and until the axis of the torus has emerged, t=54) the vertical photospheric field is set to that from the appropriate slice of the torus’s field. Dynamical equations are then solved in order for the field above $z=0$ to relax, so that at every time step the resulting configuration is a force-free equilibrium. The unsigned photospheric flux as a function of time is shown on Fig. 6.

A visual representation of characteristics times is shown on Fig. 7. Different rows correspond to different times: $t=15$ – the tube is about to emerge; $t=24$ – the minor axis of the torus has emerged; $t=32$ – the bottom of the torus has emerged; $t=45$ – the tube undergoes acceleration; $t=54$ – the major axis of the torus has emerged, the torus has stopped emerging, the tube starts getting a significant writhe; $t=58$ – the tube escapes the domain; the simulation is over. Note that the torus starts to kink at $t\geq 45$ and keeps kinking until it escapes the computational domain at $t=58$.

3.2 Computing $H_{A}$ For Given Volume And The Potential Field.

We define different domains, ${\cal D}$, with the same field, by making a different choices of boundary mask. We were interested in how different portions of the torus, namely, the “core” and the outer layers behave during the instability.

Our masks are defined to be within the photospheric intersection of the emerging torus, $\varpi\leq\varpi_{max}$. By choosing different values of $\varpi_{max}$ we construct domains, containing different portions of the emerging flux tube The footpoints of domains with different $\varpi_{max}$ are shown on Fig. 8. The shape is distorted with respect to the original cross-section of a torus due to reconnection with the arcade, current sheet formation and due to near horizontality of some field lines.

We found domains for masks with $\varpi_{max}\in\left[0.5,1.0,2.0\right]R$ at different times during the emergence. We computed $\Theta\left(\varpi_{max},t\right)$ and then constructed a potential field confined to it. The results are shown in Fig. 8, Fig. 9 and Fig. 10.

For each $t$ and $\varpi_{max}$ we calculated vector potentials of the actual field, $\Theta\left(t,\varpi_{max}\right){\bf B}\left({\bf r},t\right)$, and the reference field ${\bf P}\left({\bf r},t,\varpi_{max}\right)$. To do this we used a gauge in which one of the components of the vector potential (in our case, $A_{z}$) is identically zero. The other two could be found with a straight-forward computation:

[TABLE]

In terms of these elements the addirive self helicity the additive self-helicity:

[TABLE]

is computed.

3.3 Measuring Twist in Thin Flux Tube Approximation

To make contact with previous work we compare the additive self helicity to the twist helicity in our flux bundeles. It can be shown analytically that in the limit of a vanishingly thin flux tube these quantities are identical. Here we must compute twist helicity for flux bundles of non-vanishing width. We do this in terms of a geometrical twist related to twist helicity.

One cannot really speak of twist, or of an axis, in the domains defined above. First, the thickness and the curvature radius of the flux bundles are comparable to their lengths. Secondly, the magnetic field and the twist vary rapidly over the cross section of the bundle.

The domains constructed from the smaller masks, $\varpi_{max}=0.5$ and $\varpi_{max}=1.0$, may, however, be suitable for approximation as thin tubes. Even in these cases the approximation may suffer near the top part at later times: at $t=50$ the radius of curvature becomes comparable to the width, and later, during kinking the radius of the tube becomes comparable to the length (see Fig. 8 and Fig. 7).

We define an axis for the flux bundle by first tracing many field lines within it. Then we divide each field line into $N$ equal segments ($N$ is the same for all lines) of length $L_{i}/N$, where $L_{i}$ is the length of the $i^{\rm th}$ line. If the bundle were an ideal cylinder, the midpoints of the $n^{\rm th}$ segment from every line would lie on a single plane; provided the bundle is thin the these midpoints will lie close to a plane. We define the $n^{\rm th}$ point on an axis by the centroid of these approximately co-planar points. The set of $N$ centroids forms the axis of our tube.

We then define the tangent vector $\textit{{\^{l}}}_{i}$ along this axis, and a plane normal to this vector and thus normal to the flux tube (at least in the thin flux tube approximation). If the tube has some twist in it, then the point where one field line intersects the plane will spin about the axis as the plane moves along the tube. Such spinning must be defined relative to a reference vectore on the plane which “does not spin”. The net angle by whcih the intersection point spins, relative to the non-spinning vecotr, is the total twist angle of the tube. In a thin tube all field lines will spin by the small angle; in our general case we compute an everage angle.

We produce a non-spinning reference vector using an orthonormal triad, arbitrarily defined at one end of the tube, and carried along the axis by parallel transport. For a curve with tangent unit vector l̂, the parallel transport of a vector ${\bf u}$ means $\textit{{\^{l}}}\cdot\left(\partial{\bf u}/\partial l\right)=0$. To impliment this numerically an arbitrary unit vector ${\bf\hat{u}}_{0}$ is chosen at one end of the axis perpendicular to the tangent, ${\bf\hat{u}}_{0}\cdot\textit{{\^{l}}}_{0}=0$. The third member of the triad is ${\bf\hat{v}}_{0}={\bf\hat{u}}_{0}\times\textit{{\^{l}}}_{0}$. At the next point, ${\bf\hat{u}}_{1}$ is chosen by projecting ${\bf\hat{u}}_{0}$ onto a plane normal to $\textit{{\^{l}}}_{1}$ and normalizing it

[TABLE]

(see Fig. 11). Then ${\bf\hat{v}}_{1}={\bf\hat{u}}_{1}\times\textit{{\^{l}}}_{1}$, and the procedure is repeated for every segment along the axis.

4 Results

Based on the analogy between $Tw$ and $H_{A}/\Phi^{2}$, it would be natural to introduce quantity analogous to $L$ and $Wr$ in a similar way. We propose that $L$ in the general (non-“thin”) case might be analogous to the unconfined self-helicity, introduced in Longcope & Malanushenko, 2007, and $Wr$ is similar to the helicity of the confined potential field relative to the unconfined potential field. From equation (3) of Longcope & Malanushenko, 2007

[TABLE]

(where ${\bf x}={\bf r}(x,y,0)$ and ${\bf P}_{\cal V}$ is a potential field confined to ${\cal V}$ that matches boundary conditions ${\bf B}\cdot{\bf\hat{n}}|_{\partial{\cal V}}={\bf P}_{\cal V}\cdot{\bf\hat{n}}|_{\partial{\cal V}}$) by plugging it into $H(\Theta_{\cal D}{\bf B},{\bf P}_{\cal D},{\cal V})$ and $H({\bf B},{\bf P}_{{\cal V},\Theta},{\cal V})$ and adding them together it immediately follows, that

[TABLE]

where ${\cal D}\subset{\cal V}$ and $\Theta_{\cal D}$ is a support function of ${\cal D}$. By ${\bf P}_{\cal D}$ we mean the potential field confined to ${\cal D}$ (and identically zero outside of ${\cal D}$) that matches boundary conditions ${\bf B}\cdot{\bf\hat{n}}|_{\partial{\cal D}}={\bf P}_{\cal D}\cdot{\bf\hat{n}}|_{\partial{\cal D}}$, and by ${\bf P}_{{\cal V},\Theta}$ we mean the potential field, confined to ${\cal V}$, that matches boundary conditions $\Theta_{\cal D}{\bf B}\cdot{\bf\hat{n}}|_{\partial{\cal V}}={\bf P}_{{\cal V},\Theta}\cdot{\bf\hat{n}}|_{\partial{\cal V}}$. As long as ${\cal D}$ is fully contained in ${\cal V}$, which is constant in time, the quantity $H_{\rm unc,{\cal V}}/\Phi^{2}\equiv H\left(\Theta_{\cal D}{\bf B},{\bf P}_{{\cal V},\Theta},{\cal V}\right)/\Phi^{2}$ will behave like $L$ and $H\left({\bf P}_{\cal D},{\bf P}_{{\cal V},\Theta},{\cal V}\right)/\Phi^{2}$ would then behave like $Wr$.

Fig. 14 compares the generalized twist number, $H_{A}/\Phi^{2}$, with helicity, unconfined to the flux bundle’s volume, but confined to the computational domain of the simulation: $H_{\rm unc,box}$. In this case ${\cal V}$ is the computational domain, a rectangular box. The behaviour of all quantities matches expecation: $H_{\rm unc,box}/\Phi^{2}$ increases as the torus emerges, and stays nearly constant after the emergence is complete (the slight decrease is due to the reconnection with the arcade field). The generalized twist number, $H_{A}/\Phi^{2}$ also increases with the emergence, but decreases between $t=50$ and $t=58$ – the time when the torus kinks (see Fig. 7). For different $\varpi_{max}$ the decrease seems to start at a slightly different time.

Fig. 14 demonstrates as well, that the general behaviour of $H_{\rm unc,{\cal V}}/\Phi^{2}$ is qualitatively similar whether the volume ${\cal V}$ over which unconfined helicity is computed is the computational domain or the half space. To compute the unconfined helicity in the half space, $H_{\rm unc,Z_{+}}$, we integrate the helicity flux in the way described in (DeVore, 2000) and used in (Fan & Gibson, 2004). The helicity flux is computed relative to the potential field in half space, and thus, the helicity flux, obtained in this way, might be considered a “confined to a half space”.

Longcope and Malanushenko (2008) show that $H_{\rm unc,box}=H_{\rm unc,Z_{+}}$ when the volumes, ${\cal V}$ and $Z_{+}$ and the vertical field, $B_{z}(z=0)$, all share a reflectional symmetry. This situation occurs in the simulation only for $t\geq 54$ when the torus is fully emerged and its major axis is at the photosphere. At these times the vertical component of the field is the toroidal component of the torus, which is symmetric about $y=0$. Due to reconnection with the arcade, however, the footpoints of ${\cal D}$ may not share this symmetry, in which case the photospheric field $\Theta_{\cal D}B_{z}$ is not precisely symmetric. If the two helicities were ever to coincide, it would be at $t=54$, so we choose constant of integration by setting $H_{\rm unc,box}=H_{\rm unc,Z_{+}}$ at that time. The time histories of both unconfined helicities are plotted in Fig. 14. The discrepancy between the two before $t=54$ arises from the non-vanishing helicity of ${\bf P}_{{\cal V},\Theta}$ relative ${\bf P}_{Z_{+},\Theta}$ owing to a photospheric field, $B_{z}$, lacking reflectional symmetry. In spite of the discrepancy, we draw from each curve the same basic conclusion, that the kink deformation of ${\cal D}$ does not change $H_{unc,{\cal V}}$.

Fig. 15 compares the generalized twist number to the traditional twist number described above. The twist number was computed only for the thinner subvolumes of the torus, $\varpi_{max}=0.5R$ nd $\varpi_{max}=R$. Fig. 15 shows agreement quite well for $\varpi_{max}=R$ and less well for $\varpi_{max}=0.5R$. The reason might be the following: the smaller the subvolume, the the fewer points does it have, so that, first, there are fewer field lines to be traced to measure twist, and second, the potential field, obtained by relaxation is numerically less precise. Nevertheless, the magnitudes and the general behaviors do agree.

Fig. 15 also shows the twist number measured for the potential field in a subvolume ${\bf P}$, is zero to measurement error. Note, that a significant portion of the torus is emerged, its length is not large enough (relative to the thickness) for the thin tube approximation to be valid. As the twist of the potential field should theoretically be zero (as well as generalized twist), this plot also gives an idea of the magnitude of the error of twist measurements; at most times the error is less than 15% of the value.

5 Discussion

We have demonstrated that, at least in one MHD simulation, the quantity, $Tw_{(gen)}$, defined in terms of the additive self helicity shows a threshold beyond which the system became dynamically unstable. The simulation we considered, originally studied by Fan & Gibson (2003), is a three-dimensional, numerical solution of the time- dependent, non-linear evolution of an emerging flux system. The original study established that the system became unstable to a current-driven (kink) mode at some point during its evolution. In this work we have shown that the quantity $Tw_{(gen)}$ increases until the instability ($Tw_{(gen)}\simeq 1.5$) at which time it drops. This drop occurs as a natural consequence of the instability itself.

The quantity we propose as having a threshold, $Tw_{(gen)}$, is computed using a version of the self helicity previous defined by Longcope & Malanushenko (2008). The present work has provided a detailed method for computing this quantity for any complex bundle of field lines within a magnetic field known on a computational grid. We also demonstrate that for the very special cases when that bundle can be approximated as a thin flux tube, $Tw_{(gen)}$ is approximately equal to the traditional twist number, $Tw$. In the case of thin flux tubes which are also dynamically isolated, free magnetic energy is proportional to $(Tw)^{2}$. Their free energy may be spontaneously reduced if and when it becomes possible to reduce the magnitude of $Tw$ at the expense of the writhe number, $Wr$, of the tube’s axis.

All this supports the hypothesis that $Tw_{(gen)}$ could be treated as a generalization of $Tw$. Such a generalization might be extremely useful in predicting the stability of magnetic equilibria sufficiently complex that they cannot be approximated as thin flux tubes. The case we studied, of a thick, twisted torus of field lines (Fan & Gibson, 2003), appears to become unstable when $Tw_{(gen)}$ exceeds a threshold value between $1.4$ and $1.7$. This value happens to be similar to the threshold on $Tw$ for uniformly twisted, force-free flux tubes, $Tw\approx 1.6$, as $\Delta\theta\approx 3.3\pi$ (Hood & Priest, 1979).

Previous investigations have shown that the threshold on $Tw$ depends on details of the equilibrium such as internal current distribution (Hood & Priest, 1979). It is reasonable to expect the same kind of dependance for any threshold on $Tw_{(gen)}$, so we cannot claim that $Tw_{(gen)}<1.7$ for all stable magnetic field configurations. To investigate such a claim is probably intractable, but useful insights may be obtained by applying the above analysis to magnetic equilibria whose stability to the current-driven instability is already known. The paucity of closed-form, three-dimensional equilibria in the literature, and far fewer stability analyses of them, suggests this may be a substantial undertaking.

• Baty (2001)

Baty, H. 2001, A&A, 367, 321

• Berger & Field (1984)

Berger, M. A., & Field, G. B. 1984, JFM, 147, 133

• Bernstein et al. (1958)

Bernstein, I. B., Frieman, E. A., Kruskal, M. D., & Kulsrud, R. M. 1958, Proc. Roy. Soc. Lond., A244, 17

• DeVore (2000)

DeVore, C. R. 2000, ApJ, 539, 944

• Fan & Gibson (2003)

Fan, Y., & Gibson, S. E. 2003, ApJ, 589, L105

• Fan & Gibson (2004)

Fan, Y., & Gibson, S. E. 2004, ApJ, 609, 1123

• Finn & Antonsen (1985)

Finn, J., & Antonsen, T. M., Jr. 1985, Comments Plasma Phys. Controlled Fusion, 9, 111

• Hood & Priest (1979)

Hood, A. W., & Priest, E. R. 1979, Solar Phys., 64, 303

• Hood & Priest (1981)

Hood, A. W., & Priest, E. R. 1981, Geophys. Astrophys. Fluid Dynamics, 17, 297

• LeVeque (1955)

LeVeque, R. J. 1955, Finite Difference Methods for Ordinary and Partial Differential Equations, Steady-State and Time-Dependent Problems (SIAM)

• Linton & Antiochos (2002)

Linton, M. G., & Antiochos, S. K. 2002, ApJ, 581, 703

• Longcope & Malanushenko (2008)

Longcope, D. W., & Malanushenko, A. 2008, ApJ, 674, 1130

• Low (1994)

Low, B. C. 1994, Physics of Plasmas, 1, 1684

• Moffatt & Ricca (1992)

Moffatt, H. K., & Ricca, R. L. 1992, Proc. Roy Soc. Lond. A, 439, 411

• Newcomb (1960)

Newcomb, W. A. 1960, Ann. Phys., 10, 232

• Perot (2000)

Perot, B. 2000, Journal of Computational Physics, 159, 58

• Rachmeler et al. (2009)

Rachmeler, L. A., DeForest, C. E., & Kankelborg, C. C. 2009, ApJ, 693, 1431

• Török et al. (2004)

Török, T., Kliem, B., & Titov, V. S. 2004, A&A, 413, L27

• Zhang et al. (2006)

Zhang, M., Flyer, N., & Low, B. C. 2006, ApJ, 644, 575