FREELY DECAYING TURBULENCE IN FORCE-FREE ELECTRODYNAMICS

Since the Lorentz force density vanishes in FFE, no ${\boldsymbol{E}}\cdot {\boldsymbol{J}}$ work is done on the charge carriers, and the system is formally energy conserving. Nevertheless, time-dependent solutions may develop regions in which the condition E < B is violated, so no frame exists in which the electric field vanishes. This indicates a breakdown of the ideal force-free assumption, and the Ohm's law given by Equation (3) must be modified. Although nonideal force-free Ohm's laws have been proposed in the literature (Gruzinov 2007; Li et al. 2012), it is still commonplace to evolve the ideal system until such a breakdown occurs and when it does to simply reduce the magnitude of ${\boldsymbol{E}}$ to prevent E > B. The physical motivation for this prescription, which we use here, is that energy is being radiated away when charges are accelerated to short out the electric field and restore E < B. The numerical evolution scheme is described in more detail in Section 3.1, and in the Appendix we confirm that it yields numerical convergence of the energy dissipation rate.

FFE shares the same magnetic topological invariants as Newtonian and relativistic MHD, the lowest order of which is the total magnetic helicity given by Equation (1). Its invariance can be seen as stemming from the conservation of magnetic flux through a closed-field loop, which is why it is commonly referred to as characterizing the interlocking of the magnetic field. Although the helicity density ${\boldsymbol{A}}\cdot {\boldsymbol{B}}$ depends on the choice of gauge, its integral ${H}_{{ \mathcal V }}$ over any volume ${ \mathcal V }$ bounded by magnetic surfaces is well defined and also invariant under ideal evolutions, such as Faraday's law, when ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$ (e.g., Brandenburg & Subramanian 2005):

$$\begin{eqnarray}&&{\partial }_{t}{H}_{{ \mathcal V }}=-2{\displaystyle \int }_{{ \mathcal V }}{\boldsymbol{E}}\cdot {\boldsymbol{B}}{d}^{3}x.\end{eqnarray} \tag{ 5 }$$

Equation (5) implies that ${H}_{{ \mathcal V }}$ can still be conserved approximately in the presence of nonideal processes, as long as the volume in which they occur is small.

In principle, a domain that admits a partitioning by magnetic surfaces has an independently conserved helicity associated with each smaller volume. But, in practice, 3D field geometries are too complex for such a partitioning to be possible. This is what led Taylor to conclude that relaxation is only constrained by a single topological invariant: the total helicity. The story is different in two dimensions since the magnetic surfaces are far simpler; their cross sections are nested, closed curves in the x–y plane, and the helicity enclosed by each functions as an independently conserved quantity for as long as that surface retains its identity. But nonideal effects, however small, permit the surfaces to merge with one another, erasing their identities and shuffling up their conserved helicities. Nevertheless, in two dimensions we can still construct a robust topological invariant that is far stricter than the total helicity.

We recall that, in the presence of z-translational symmetry, the in-plane magnetic field is tangent to the isocontours of the magnetic flux function A_z, so each subvolume ${{ \mathcal V }}_{i}(\psi )$ in which A_z > ψ is associated with a conserved helicity ${{ \mathcal H }}_{i}(\psi )={\int }_{{{ \mathcal V }}_{i}(\psi )}{\boldsymbol{A}}\cdot {\boldsymbol{B}}{d}^{3}x$. To whatever extent the helicities of such volumes remain additive with respect to reconnections between their bounding surfaces, the sum ${ \mathcal H }(\psi )={\displaystyle \sum }_{i}{{ \mathcal H }}_{i}(\psi )$ must be a robust topological invariant. This assumption is justified by Equation (5) when the nonideal regions (where ${\boldsymbol{E}}\cdot {\boldsymbol{B}}\ne 0$) are small. Alternatively, consider the helicity (e.g., Berger & Field 1984)

$$\begin{eqnarray}&&{{ \mathcal H }}_{12}=2{{\rm{\Phi }}}_{1}{{\rm{\Phi }}}_{2}\end{eqnarray} \tag{ 6 }$$

of two flux tubes that are once-linked and independently nonhelical. An example is a cylindrical structure whose external, toroidal flux Φ₂ wraps an interior, poloidal flux Φ₁. Two such structures that were joined to share an outer surface would then enclose a poloidal flux 2Φ₁, which is wrapped by the same toroidal flux Φ₂. So, to whatever extent the "joining" is done without destroying magnetic flux, the resulting arrangement has a helicity $2{{ \mathcal H }}_{12}$.

The presence of additional invariants ${ \mathcal H }(\psi )$ is expected to place a more restrictive lower bound on the magnetic energy than would the total helicity H alone. In particular, given that Taylor states of the same helicity but different wavelength will generally span a distinct range of ψ values, there may be no way for one 2D linear equilibrium to evolve into another of a different wavelength while leaving ${ \mathcal H }(\psi )$ unchanged. With this in mind, we anticipate that fully relaxed 2D states will have an energy that exceeds the Taylor minimum of α₀H/2. We will describe the numerical determination of ${ \mathcal H }(\psi )$ in Section 3.3, and in Section 4.3 we confirm that it is conserved by directly measuring it in our 2D simulations.

There is also the question of what should happen to configurations in which H = 0 but ${ \mathcal H }(\psi )\ne 0$ since such a configuration could not attain arbitrarily low energy while respecting each of the invariants ${ \mathcal H }(\psi )$. Although behavior like this may be entirely possible, the particular initial conditions used in this study (described in Section 3.2 and summarized in Table 1) generally have values of $| { \mathcal H }(\psi )|$ that are much smaller than 2U_B/α₁, so we have not been able to observe it yet. Indeed, we will see in Section 4.2 that our 2D states with zero net helicity still decay toward very small energy.

Table 1.  Summary of Runs

Grid Resolution		α0	See Figure
5123	1	16	1(b)
5122	1	16	2
5123	1	$\sqrt{257}\approx 16$	2
5122	0	$\sqrt{257}\approx 16$	2
30722	1	256	1(a), 4, 5, 8
40962, 61442, 81922, 122882	1	256	12
163842	1	256	7, 12
2562, 5122, 10242, 20482, 40962	1	8, 16, 32, 64, 128	6
40962	1	32	3
122882	0	$\sqrt{65641}\approx 256$	2, 11, 10(a)
122882	1	$\sqrt{65641}\approx 256$	2, 11, 10(b)
10243	0	$\sqrt{2305}\approx 48$	11, 10(c)
10243	1	$\sqrt{2305}\approx 48$	11, 10(d)
163842	0	$\sqrt{65641}\approx 256$	11
163842	1	$\sqrt{65641}\approx 256$	11
7683	0	$\sqrt{2305}\approx 48$	11
7683	1	$\sqrt{2305}\approx 48$	11

Note. Shown is a summary of the runs along with the figures that reference them. Helical runs have  = 1, and nonhelical runs have  = 0. Those whose initial data is the 2D ABC field have α₀ values that are exact integers, whereas randomized initial data have α₀ values that are not exact integers.

Download table as:  ASCIITypeset image

We evolve the FFE system using a fourth-order finite difference scheme. We use standard fourth-order difference operators to evaluate all the field gradients and standard fourth-order Runge–Kutta time stepping to advance the solution in time.

The FFE system requires three vector constraints to be maintained: no monopoles, ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$, perfect conductivity, ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$, and the existence of a frame in which ${\boldsymbol{E}}$ vanishes, E < B. The first two constraints are formally preserved by FFE but can be violated numerically at the level of truncation error. Our scheme maintains the solenoidal constraint using the hyperbolic divergence cleaning scheme proposed by Dedner et al. (2002) and later used in FFE simulations (e.g., Palenzuela et al. 2010). This amounts to supplementing Faraday's law with a magnetic monopole current ${{\boldsymbol{J}}}_{B}=-{\rm{\nabla }}{\rm{\Psi }}$, where the scalar field Ψ evolves according to the damped wave equation ${\partial }_{t}{\rm{\Psi }}=-{\rm{\nabla }}\cdot {\boldsymbol{B}}-{\tau }^{-1}{\rm{\Psi }}$, with τ being a nonphysical timescale for quenching the magnetic monopole. Here, ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$ is maintained exactly by disregarding the part of the truncation error that would give rise to a component of ${\boldsymbol{E}}$ in the direction of ${\boldsymbol{B}}$. Numerical noise introduced by finite difference operations can lead to the unphysical growth of modes whose wavelengths are comparable to the numerical grid spacing. Our scheme suppresses these unphysical modes using Kreiss–Oliger dissipation, a form of low-pass filtering. Each of the procedures just mentioned supplements the FFE equations with terms at or below the level of the truncation error, so they do not modify the formal convergence order of our numerical scheme.

This numerical scheme was used in East et al. (2015), and convergence results, as well as comparisons to relativistic MHD simulations and analytical methods, can be found in that reference. It has been implemented as part of the Mara (Zrake & MacFadyen 2011) suite of relativistic turbulence codes, which has many run-time postprocessing capabilities that allow us to perform spectral and statistical analysis of the solution at a high cadence while minimizing strain on the host architecture's file system.

We start our simulations with a monochromatic magnetic field, where all of the power is at a single wavenumber magnitude α₀, and with a vanishing electric field. The general expression for our initial conditions is

$$\begin{eqnarray}{\boldsymbol{B}}({\boldsymbol{x}}) & = & \displaystyle \sum _{| {\boldsymbol{k}}| ={\alpha }_{0}}({\alpha }_{0}{{\boldsymbol{\Psi }}}_{{\boldsymbol{k}}}+\epsilon i{\boldsymbol{k}}\times {{\boldsymbol{\Psi }}}_{{\boldsymbol{k}}}){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\\ {\boldsymbol{k}}\cdot {\hat{{\boldsymbol{\Psi }}}}_{{\boldsymbol{k}}} & = & 0\\ {\hat{{\boldsymbol{\Psi }}}}_{{\boldsymbol{k}}} & = & {\hat{{\boldsymbol{\Psi }}}}_{-{\boldsymbol{k}}}^{*}\end{eqnarray} \tag{ 7 }$$

where the parameter is chosen to be either one or zero, corresponding to helical or nonhelical configurations, respectively. Helical initial configurations where α₀ > 1 are unstable equilibria (see Section 4.1), whereas the nonhelical configurations are out of equilibrium. Some of our initial conditions are randomized, having ${{\boldsymbol{\Psi }}}_{{\boldsymbol{k}}}={\hat{{\boldsymbol{e}}}}_{{\boldsymbol{k}}}{e}^{i{\phi }_{{\boldsymbol{k}}}}$ where ${\hat{{\boldsymbol{e}}}}_{{\boldsymbol{k}}}$ is a random unit vector in the plane orthogonal to ${\boldsymbol{k}}$, and ${\phi }_{{\boldsymbol{k}}}$ is a random phase. We also make use of a special case of Equation (7) known as the "ABC" solution (Arnold 1965; Dombre et al. 1986). In general this is given by

$$\begin{eqnarray}&&{{\boldsymbol{B}}}^{{\rm{ABC}}}({\boldsymbol{x}})=\left(\begin{array}{c}{B}_{3}\mathrm{cos}{\alpha }_{0}z-{B}_{2}\mathrm{sin}{\alpha }_{0}y\\ {B}_{1}\mathrm{cos}{\alpha }_{0}x-{B}_{3}\mathrm{sin}{\alpha }_{0}z\\ {B}_{2}\mathrm{cos}{\alpha }_{0}y-{B}_{1}\mathrm{sin}{\alpha }_{0}x\\ \end{array}\right),\end{eqnarray} \tag{ 8 }$$

which is highly ordered and fully helical, meaning that ${\boldsymbol{B}}={\alpha }_{0}{\boldsymbol{A}}$ (in the Coulomb gauge). In this study we will make frequent use of the case with B₁ = B₂ = 1 and B₃ = 0, which we refer to as the 2D ABC configuration.

Our results are based on simulations having a range of initial frequencies α₀ and numerical resolutions, which we will refer to by the number of grid points in each linear dimension N. In general, the quality of our results improves when we are able to simulate larger values of α₀ with more separation between the initial length scale and the domain length scale. However, features (of size ${\alpha }_{0}^{-1}$) in our initial condition need to be resolved by a certain number of grid points in order to obtain robust solutions. In the Appendix we show that 32 cells per ${\alpha }_{0}^{-1}$ are sufficient to keep the error in the global helicity conservation smaller than 1%. In 2D, we will present simulations with α₀ as large as 256 with resolutions up to 16384². In 3D, we will present simulations with α₀ as large as 48 and resolution 1024³.

We define the power spectral density of the electric, magnetic, and helicity fields, respectively, as

$$\begin{eqnarray}{P}_{E}({k}_{i}) & = & \displaystyle \frac{1}{{\rm{\Delta }}{k}_{i}}\displaystyle \sum _{{k}_{i}\lt | {\boldsymbol{q}}| \lt {k}_{i}+{\rm{\Delta }}{k}_{i}}{{\boldsymbol{E}}}_{{\boldsymbol{q}}}\cdot {{\boldsymbol{E}}}_{{\boldsymbol{q}}}^{*}/2,\\ {P}_{B}({k}_{i}) & = & \displaystyle \frac{1}{{\rm{\Delta }}{k}_{i}}\displaystyle \sum _{{k}_{i}\lt | {\boldsymbol{q}}| \lt {k}_{i}+{\rm{\Delta }}{k}_{i}}{{\boldsymbol{B}}}_{{\boldsymbol{q}}}\cdot {{\boldsymbol{B}}}_{{\boldsymbol{q}}}^{*}/2,\\ {P}_{H}({k}_{i}) & = & \displaystyle \frac{1}{{\rm{\Delta }}{k}_{i}}\displaystyle \sum _{{k}_{i}\lt | {\boldsymbol{q}}| \lt {k}_{i}+{\rm{\Delta }}{k}_{i}}{{\boldsymbol{A}}}_{{\boldsymbol{q}}}\cdot {{\boldsymbol{B}}}_{{\boldsymbol{q}}}^{*}\end{eqnarray} \tag{ 9 }$$

where ${{\boldsymbol{E}}}_{{\boldsymbol{q}}}$, ${{\boldsymbol{B}}}_{{\boldsymbol{q}}}$, and ${{\boldsymbol{A}}}_{{\boldsymbol{q}}}$ are, respectively, the electric field, the magnetic field, and the vector potential Fourier harmonics of wavenumber ${\boldsymbol{q}}$. We normalize the Fourier harmonics so that the volume-integrated electric and magnetic field energies U_E and U_B and the magnetic helicity H are given by

$$\begin{eqnarray*}{U}_{E} & = & \displaystyle \sum _{i}{P}_{E}({k}_{i}){\rm{\Delta }}{k}_{i},\\ {U}_{B} & = & \displaystyle \sum _{i}{P}_{B}({k}_{i}){\rm{\Delta }}{k}_{i},\\ H & = & \displaystyle \sum _{i}{P}_{H}({k}_{i}){\rm{\Delta }}{k}_{i}.\end{eqnarray*}$$

We also define the characteristic frequency of each field k_E, k_B, and k_H as

$$\begin{eqnarray}&&{k}_{X}=\displaystyle \frac{{\displaystyle \sum }_{i}{P}_{X}({k}_{i}){k}_{i}{\rm{\Delta }}{k}_{i}}{{\displaystyle \sum }_{i}{P}_{X}({k}_{i}){\rm{\Delta }}{k}_{i}}\end{eqnarray} \tag{ 10 }$$

where X is one of E, B, or H. The most probable wavenumber, where P_X(k) is maximal, is denoted by ${\tilde{k}}_{X}$.

In two dimensions, we track the "helicity mass" function discussed in Section 2.2:

$$\begin{eqnarray}&&{ \mathcal H }(\psi )=\int {\rm{\Theta }}({A}_{z}({\boldsymbol{x}})-\psi ){\boldsymbol{A}}\cdot {\boldsymbol{B}}{d}^{3}x,\end{eqnarray} \tag{ 11 }$$

where Θ is the Heaviside step function. In practice, this diagnostic is more easily computed as the "helicity density" function $d{ \mathcal H }/d\psi$, which we calculate by binning the lattice points according to their value of A_z and assigning the weight ${\boldsymbol{A}}\cdot {\boldsymbol{B}}$. We also create the volume distribution $d{ \mathcal V }/d\psi$ by binning points according to A_z with uniform weights and the helicity distribution over volume $d{ \mathcal H }/d{ \mathcal V }=\frac{d{ \mathcal H }}{d\psi }/\frac{d{ \mathcal V }}{d\psi }$.

Our helical initial conditions are linear force-free equilibria. Clearly they are stable when α₀ = 1 since such states are global energy minima for a given magnetic helicity. The question of the ideal stability of the periodic shorter-wavelength (α₀ > 1) Taylor states has a conflicted history (e.g., Moffatt 1986; Galloway & Frisch 1987; Er-Riani et al. 2014), which has recently been resolved in East et al. (2015). In that study, numerical simulations of both FFE and relativistic MHD revealed that α₀ > 1 Taylor states are linearly unstable to ideal perturbations. The instability is marked by exponential growth of the electric field on roughly the Alfvén wave crossing time of the initial structure size ${\alpha }_{0}^{-1}$ and saturates when the medium attains the Alfvén speed, which for the magnetically dominated case is c, implying the existence of regions where E ≈ B. This instability affects generic states, and the only counterexamples that were found were one-dimensional ABC states (e.g., B₁ = 1, B₂ = 0, B₃ = 0) that are stable for all values of α₀. Such states are pure plane waves having circular polarization and are force-free by virtue of having uniform magnetic pressure. All of our helical initial conditions are short wavelength and either 2D or 3D, and turbulent relaxation begins after the saturation of the ideal instability.

Here we discuss the magnetic energy associated with the end-state magnetic configurations. Since the Taylor states have ${\boldsymbol{B}}={\alpha }_{0}{\boldsymbol{A}}$, their energy is given simply by α₀H/2. In other words, their energy is α₀ times larger than the theoretical lower limit imposed by assuming the state reaches α = 1 at constant H. Whether or not dynamical relaxation processes settle with the field in this global energy minimum remains an open question, but here we provide some evidence to support the following conjecture: "Force-free magnetic relaxation starting from periodic Taylor states ends in the lowest-energy configuration allowed by helicity conservation if and only if the domain is 3D."

We have carried out a suite of calculations belonging to one of four categories, being either 2D or 3D and either helical or nonhelical. Those that are nonhelical are, by construction, out of equilibrium at t = 0 and could decay until they reach zero energy because helicity conservation does not place any lower bound on their magnetic energy. The helical ones are initially at an unstable equilibrium, enter a period of turbulent relaxation, and settle in a force-free equilibrium of lower energy. We considered initial conditions that are both of the randomized type given by Equation (7) and the ABC type given by Equation (8).

Our results are summarized in Figure 2, which shows the time series U_B for each of six different runs. As expected, the nonhelical configurations decay toward zero energy in both 2D and 3D settings. In three dimensions, the helical configurations all terminate in the lowest-energy state allowed by helicity conservation.² Both a randomized setup with ${\alpha }_{0}^{2}=257$ and the 2D ABC setup with ${\alpha }_{0}^{2}=256$ showed the same general behavior. The latter setup was intentionally chosen to be identical with the 2D setup apart from the inclusion of a low-level (one part in 10⁸) white-noise perturbation introduced to break the z-translational symmetry.

Zoom In Zoom Out Reset image size

Both of the helical 2D runs terminate their relaxation with an energy that is decreased to only 30% of its initial value. For example, a randomized 2D initial condition with α₀ ≈ 256 settles in a state that has roughly 77 times more magnetic energy than the Taylor minimum-energy state. This is not unique because actually all of our helical 2D runs where α₀ ≫ 1 (including α₀ = 16) settle in a state whose energy is decreased to roughly 30%. In the Appendix we show that, for the 2D ABC α₀ = 256 setup, the final energy is numerically converging to a value very near 30%. This means that the terminal states in two dimensions are not linear equilibria. In other words, they do approximately solve the force-free condition of Equation (2), but the value of α may vary from one field line to another.

The fact that the 2D configurations do not relax into linear equilibria stems from invariance of the helicity distribution ${ \mathcal H }(\psi )$ that we introduced in Section 2.2. Figure 3 confirms that it does indeed remain constant over time to a very good approximation, even while the volume distribution over the magnetic potential $d{ \mathcal V }/d\psi$ changes significantly. The feature evident in the bottom panel of Figure 3, where in the relaxed state most of the volume is occupied by the extreme values of the magnetic potential, can be connected to the formation of current layers, which we discuss in the next section.

Zoom In Zoom Out Reset image size

As an illustration that invariance of the helicity distribution might prevent 2D relaxation from attaining longer-wavelength linear equilibria, consider that the magnetic potential function A_z of any 2D linear equilibrium satisfies

$$\begin{eqnarray*}&&{A}_{z}={\left(\displaystyle \frac{H}{2{\alpha }_{0}{L}^{3}}\right)}^{1/2}\displaystyle \frac{{B}_{z}}{{\bar{B}}_{z}}\end{eqnarray*}$$

where ${\bar{B}}_{z}$ is the root mean square value of B_z, and H is the total helicity. Note that ${ \mathcal H }(\psi )$ can be characterized by its domain, namely the global maximum of $| {A}_{z}|$, which we denote by ψ_max. For two different linear equilibria with frequencies α₀ and ${\alpha }_{0}^{\prime }$ to have the same helicity distribution, it would be necessary for them to at least share the same values of H and ψ_max, requiring that

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\alpha }_{0}}{{\alpha }_{0}^{\prime }}\right)}^{1/2}=\displaystyle \frac{{B}_{z,{\rm{max}}}/{\bar{B}}_{z}}{{B}_{z,{\rm{max}}}^{\prime }/{\bar{B}}_{z}^{\prime }},\end{eqnarray} \tag{ 12 }$$

where B_z,max is the global maximum of $| {B}_{z}|$. Note that, in general, ${\bar{B}}_{z}\leqslant {B}_{z,{\rm{max}}}$. For the particular case of the 2D ABC state in which B₁ = B₂ = 1, Equation (12) leads to the requirement that ${\alpha }_{0}^{\prime }\geqslant {\alpha }_{0}/2$, and therefore its relaxation into a linear state with two times smaller energy is impossible. In other words, there is no way for the 2D ABC state represented in Figure 3, whose α₀ = 32, to evolve into another linear equilibrium whose frequency is smaller than 16 while preserving ${ \mathcal H }(\psi )$. We suspect this argument could be generalized further, but for now we leave it as a conjecture that the wavelength of a linear 2D equilibrium may be uniquely specified by its helicity distribution, which if true would render it impossible for one linear 2D equilibrium to relax into another of lower energy.

During the turbulent relaxation of helical 2D configurations, the solution consists of oppositely signed flux domains separated from one another by a network of rotational force-free current layers. The flux domains (black and white regions of Figure 1) have nearly uniform B_z and are thus relatively current-free. Across the current layers, the magnetic field direction rotates through approximately 180°, while its magnitude (and thus magnetic pressure) remains fixed. One such current layer is shown in Figure 4, where a one-dimensional profile has been taken along the x axis, passing through the layer where it is aligned with the y axis.

Zoom In Zoom Out Reset image size

It is evident that the current layers have a characteristic frequency, which we denote by α_c and determine empirically as follows. Since the solution is near a force-free equilibrium, the current is ${\boldsymbol{J}}\approx \alpha {\boldsymbol{B}}$ for some spatially dependent frequency

$$\begin{eqnarray*}&&\alpha ={\boldsymbol{B}}\cdot {\rm{\nabla }}\times {\boldsymbol{B}}/{B}^{2}.\end{eqnarray*}$$

We anticipate that the probability density function P(α) will have two local maxima: one at α = 0, corresponding to the potential flux domain interiors, and the other at the frequency α_c, marking the frequency of the current layers. Figure 5 confirms this to be the case. It shows P(α) at logarithmically spaced times throughout the relaxation and reveals the location of the second peak once the solution is sufficiently close to a force-free equilibrium. For this particular run, with α₀ = 256 and N = 3072, the value of α_c ≈ 128.

Zoom In Zoom Out Reset image size

It turns out to be a coincidence that in this case α_c ≈ α₀/2. We have performed a family of calculations, varying α₀ between 8 and 128 and varying N between 128 and 8192. For each run, we recorded the value of α_c, time-averaged over roughly 100 snapshots between time t = 10 and t = 16, which was late enough that the second peak in P(α) had emerged in each run. There was no secular evolution of α_c. Figure 6 reveals that current layers become increasingly narrow with higher resolution but also with increasing initial frequency α₀. The scaling is consistent with the expression

$$\begin{eqnarray}&&{\alpha }_{c}={k}_{1}^{3/4}{\alpha }_{0}^{1/4}\end{eqnarray} \tag{ 13 }$$

where k₁ = N/30 is the turbulence cutoff frequency, which has been determined by modeling the magnetic energy spectrum (see Figure 7) and is insensitive to initial conditions, depending only on the numerical scheme and grid resolution. We emphasize that Equation (13) is strictly empirical, in that it matches the data shown in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The scaling given by Equation (13) indicates that, in the infinite resolution limit, the current layers will have zero characteristic length. We are not able to say whether the scaling could be associated with a physical property of 2D FFE at a finite Reynolds number or if it might depend on the details of the numerical scheme. This question could be resolved by imposing the turbulence cutoff frequency k₁ explicitly and then varying the numerical resolution. In other words, solutions to a resistive FFE system at a given conductivity parameter will be necessary.

The flux domain interiors contain another type of coherent structure. They are long-lived bubbles of toroidal magnetic field that are confined by ambient magnetic pressure. The circular object to the left of the current layer in Figure 4 is an example. They may be referred to as flux tubes or magnetic vortices, but we will refer to them as magnetic bubbles because we believe they are objects similar to those studied by Gruzinov (2010). The bubbles are helical, force-free magnetic field structures, having ${\boldsymbol{J}}$ very well aligned with ${\boldsymbol{B}}$. But they are not linear equilibria: the value of α varies with distance r from the axis, and all the bubbles we examined had similar radial profiles α(r). The current flow is parallel to the background magnetic field near the core, but an equal and opposite return current flows through the sheath, so they could also be thought of as coaxial electric current channels oriented along the symmetry axis.

Force-free configurations with axial and z symmetry satisfy the ordinary differential equations

$$\begin{eqnarray}{B}_{\phi }\alpha +{B}_{z}^{\prime } & = & 0\\ {B}_{\phi }{r}^{-1}-{B}_{z}\alpha +{B}_{\phi }^{\prime } & = & 0\end{eqnarray} \tag{ 14 }$$

where α(r) needs to be specified to make the solution unique, as well as the boundary conditions ${B}_{z}(r)\to {B}_{\infty }$ and ${B}_{\phi }(r)\to 0$ as $r\to \infty$. Requiring the total current $I=\int {B}_{\phi }{rd}\phi$ to be zero only forces B_ϕ(r) to decrease faster than 1/r. Alternatively, the magnetic pressure u_B(r) may be specified instead of α(r). One simple assumption, consistent with measurements of the bubbles, is that the magnetic pressure enhancement, relative to the ambient pressure ${B}_{\infty }^{2}/2$, is Gaussian. This gives a pressure profile

$$\begin{eqnarray}&&{u}_{B}(\tilde{r})=\displaystyle \frac{1}{2}{B}_{\infty }^{2}(1+{{fe}}^{-{\tilde{r}}^{2}})\end{eqnarray} \tag{ 15 }$$

in the dimensionless radius $\tilde{r}={\alpha }_{b}r$ for a bubble of radius ${\alpha }_{b}^{-1}$, where f is the magnetic pressure enhancement relative to the background and could be any positive number. The corresponding solution of Equation (14) is

$$\begin{eqnarray*}{B}_{\phi }(\tilde{r}) & = & {B}_{\infty }{f}^{1/2}\tilde{r}{e}^{-{\tilde{r}}^{2}/2}\\ {B}_{z}(\tilde{r}) & = & {B}_{\infty }\sqrt{1+{{fe}}^{-{\tilde{r}}^{2}}(1-{\tilde{r}}^{2})}\end{eqnarray*}$$

which has the profile

$$\begin{eqnarray}&&\alpha (\tilde{r})={\alpha }_{b}\displaystyle \frac{(2-{\tilde{r}}^{2}){e}^{-{\tilde{r}}^{2}/2}}{\sqrt{(1-{\tilde{r}}^{2}){e}^{-{\tilde{r}}^{2}}+{f}^{-1}}}.\end{eqnarray} \tag{ 16 }$$

In Figure 8 we show the radial profile of magnetic pressure, azimuthally averaged around the bubble's axis. This is the same object depicted to the left of the current layer in Figure 4. We also show the expression given in Equation (15) with its best-fit model parameters. For this particular object, the best-fit fractional magnetic energy enhancement was f = 0.41, and the best-fit inverse radius was α_b = 162. The lower panel of Figure 8 shows the radial profile of α and also the model-predicted value given by Equation (16) with the best-fit parameters.

Zoom In Zoom Out Reset image size

So far we have not detected these structures in 3D simulations, but that does not mean they never happen. It is crucial to address whether they are even stable in three dimensions, and if they are, then under what conditions they may be attractors.

The inverse cascade can be characterized by the rate with which the magnetic coherence scale ${k}_{B}^{-1}$ (Equation (10)) migrates toward longer wavelengths. We observe inverse cascading in FFE for every setting we have considered, whether the turbulence is 2D or 3D, helical or nonhelical. Figure 9 shows the evolution of the magnetic energy spectrum P_B(k) over time, for representative helical and nonhelical runs in 2D. We note how, in both cases, the remaining magnetic energy resides at an increasing scale over time. We also note that the spectral energy density at long wavelengths is an increasing function of time, indicating that selective decay of the short-wavelength modes alone cannot explain growth of the spatial coherency. This further generalizes the observations made of nonhelical 3D MHD turbulence in both Newtonian (Brandenburg et al. 2015) and relativistic (Zrake 2014) settings.

Zoom In Zoom Out Reset image size

Figure 10 shows the evolution of the characteristic magnetic frequency k_B for one model from each of the four categories. For the helical runs we also show the evolution of the helicity frequency k_H. Except for the general feature that, over time, both k_B and k_H move to smaller frequency, there is no single decay law that characterizes all these settings. Some runs exhibit a power-law time dependence of U_B, k_B, or k_H, but not all. Power-law dependence is considered evidence of self-similarity in the relaxation process, meaning the solution evolves only by rescaling itself in space and time (e.g., Landau & Lifshitz 1987, p. 227). Self-similar evolution may occur when the characteristic scales are all sufficiently smaller than the domain size so that processes around the coherence scale are not yet contaminated by the requirement of periodicity at the domain scale. For this reason, large values of α₀, and thus high domain resolution, may be required for self-similarity to emerge. For 3D, freely decaying, nonhelical relativistic MHD turbulence, Zrake (2014) found a power-law dependence of ${k}_{B}\propto {t}^{-2/5}$ for α₀ ≈ 48 initial conditions. For the same conditions in FFE, k_B evolves with a power-law index of −0.48, as shown in the second column of Figure 10. However, the dependence of k_B on time is not convincingly power law, indicating that a scale separation of α₀ ≈ 48 may not be sufficient to yield a decisive measurement of the decay index. In the future, we plan to follow up on this with higher-resolution simulations.

Zoom In Zoom Out Reset image size

In freely decaying, 2D, nonhelical force-free turbulence with very high resolution (12288²) and ${\alpha }_{0}\approx 256$, a clear power-law behavior of k_B is observed with an index of −0.38. The helical 2D case is different. As shown in Figure 10 (column 3), the energy drops to its terminal value of ≈0.3 U_B (t = 0) long before ${k}_{B}^{-1}$ approaches the domain scale. At times later than t = 0.15, the merging of flux domains (evident in Figure 1) moves the magnetic energy to progressively larger scales while suffering slower and slower dissipative losses. During this epoch, both the helicity and magnetic frequencies decay with a power-law index of roughly −0.55.

In general, freely decaying magnetic turbulence must exhibit an inverse cascade when the magnetic helicity is near maximal (H ≈ U_B/k_B) if ${\partial }_{t}H=0$ is to be satisfied (Frisch et al. 1975; Christensson et al. 2001; Cho 2011). However, the same conclusion cannot be drawn from magnetic helicity conservation alone when H ≪ U_B/k_B, as is the case for our nonhelical runs. Observation of the nonhelical inverse cascade was thus seen as a surprise (Zrake 2014; Brandenburg et al. 2015). It was suggested in Zrake (2014) that the nonhelical inverse cascade may reflect the tendency for aligned current structures to attract one another. Brandenburg et al. (2015) found that the net transfer of energy from small to large scales was about twice larger in helical than in nonhelical, freely decaying, Newtonian MHD turbulence.

In all of our initial conditions, the magnetic energy is concentrated around a single frequency, P_B(k) ∝ δ(k − α₀). As turbulence develops, the energy redistributes itself over all available scales, with the bulk of the energy around k_B (which decreases over time, as discussed in Section 4.6) and a power-law tail extending to the spectral cutoff frequency k₁, which lies consistently near N/30 when the grid resolution is N. We have determined the index s of the power-law tail by fitting the logarithm of the spectral distributions P_E(k) and P_B(k) to the model function

$$\begin{eqnarray}&&f(x)={{Ae}}^{-{(x-{x}_{0})}^{2}/2{\sigma }^{2}}+{xs}\end{eqnarray} \tag{ 17 }$$

where $x=\mathrm{log}k$, and A, x₀, σ, and s are model parameters. We have found it expedient not to try to model the high-frequency cutoffs; we just use frequency bins between the peak frequency ${\tilde{k}}_{B}$ and the cutoff k₁ to obtain our fit. A representative spectral fit is shown in Figure 7.

Figure 11 shows P_E(k) and P_B(k) for one model from each of the categories 2D/3D and helical/nonhelical along with the best-fit model given by Equation (17). We use randomized initial conditions for each model, and the resolution is 16384² in 2D and 1024³ in 3D (coinciding dash-dotted curves are taken from simulations whose resolution is 12288² and 768³ to indicate numerical consistency). The spectra represent the solution at an intermediate time, when k_B ≈ 8, and we find the power-law indices to be quite stable³ during a window of time beginning on the snapshot chosen for each model. The power-law index of the electric field spectral energy distribution P_E(k) is between 0.96 and 1.18 for the various models, with both of the helical models having an index of 1.18. The power-law index of the magnetic field spectral energy distribution P_B(k) is different in two and three dimensions. In both helical and nonhelical 3D settings, it is notably close to the Kolmogorov value of 5/3. Note that Zrake (2014) and Brandenburg et al. (2015) both measured an index near 2 for freely decaying nonhelical MHD turbulence. The helical 2D model has an index of 1.93, and it is tempting to conclude that the true value is 2, but actually the best-fit parameter for all resolutions up to 16384² is significantly smaller than 2, always being between 1.92 and 1.96. The nonhelical 2D model is found to have an index of 1.41.

Zoom In Zoom Out Reset image size

We have examined turbulent relaxation in FFE with the motivation of elucidating the physical mechanism behind the extremely fast time variability that is characteristic of astrophysical gamma-ray emitters, including the Crab Nebula, many blazars, and nearly all gamma-ray bursts. The extreme energetics and temporal intermittency of gamma radiation from these sources require a mechanism in which plasma promptly converts the majority of its magnetic energy into high-energy particles and radiation. Furthermore, the emitting regions are thought to be strongly magnetized and are known to lie a great distance from the primary mover (pulsar, progenitor star, or black hole). These facts are highly suggestive that electromagnetic outflows may contain persistent magnetic structures with copious free energy supplies, whose spontaneous disruption could be linked to the observed flaring events. A scenario like this was referred to in Blandford et al. (2014) as magnetoluminescence. For it to be plausible, it is necessary (1) that metastable, force-free (or hydromagnetic) equilibria can exist far from any supporting boundaries, (2) that such objects can form under realistic astrophysical conditions, and (3) that upon their disruption, magnetic energy is promptly and completely dissipated.

In this paper, we have begun to address points (1) and (3) and found results that are at least partially encouraging. All of our periodic 3D simulations exhibit prompt relaxation into the Taylor minimum energy state, supporting the idea that magnetic energy can be dissipated completely in a light-travel time. But the same result suggests that persistent, metastable structures are not a generic outcome of turbulence in FFE. This is not to say that such behavior is impossible, as we have only considered a small class of initial conditions. In fact, Smiet et al. (2015) very recently identified magnetic arrangements that, in 3D periodic domains, in full MHD, relax to nonlinear hydromagnetic equilibria. So (1) is possible, at least in the hydromagnetic case. The volatility of such objects and the generality of conditions under which they may arise remain important questions for the future.

In this study we have begun to address the question of whether relativistic, force-free magnetic relaxation in periodic domains generically ends in a Taylor state, and we found evidence to support the view that it does, provided the domain is 3D. But thus far we have only considered a restricted class of initial conditions, namely, isotropic, monochromatic fields that are either linear force-free equilibria (helical) or completely nonhelical.

Several studies in full MHD have now identified settings that relax to more general hydromagnetic equilibria for which ${\boldsymbol{J}}\times {\boldsymbol{B}}={\rm{\nabla }}p$, where the Lorentz force density balances the gradient of gas pressure p. Examples include stratified 3D environments (Braithwaite 2006, 2008) where the field is helical and 2D periodic settings where the field is incompressible and nonhelical (Gruzinov 2009). Amari & Luciani (2000) and Braithwaite (2015) have both provided examples of 3D hydromagnetic relaxation that first develop current layers and then proceed to a smooth configuration via resistive processes. Both studies used boundary conditions where at least one of the directions was not periodic. Very recently, Smiet et al. (2015) has found instances of 3D hydromagnetic relaxation ending in smooth, nonlinear equilibria with nonuniform pressure, even when the boundary is periodic or open. Such boundary conditions are highly relevant for the astrophysical processes mentioned in Section 5.1, and an important question is whether their results possess any force-free analogs. In particular, if stable and nonlinear force-free equilibria do exist in three dimensions away from boundaries, then how likely are they to arise under realistic astrophysical conditions?

Many studies of relativistic plasma and MHD processes are carried out assuming either translational or rotational symmetry, including simulations of FFE (McKinney 2006b; Tchekhovskoy et al. 2008), relativistic MHD (Barkov & Komissarov 2008; Komissarov et al. 2009; Komissarov & Barkov 2009; Mizuno et al. 2011), and particle-in-cell (PIC) simulations (Spitkovsky 2008; Keshet et al. 2009). As we have seen, 2D magnetic relaxations are far more topologically constrained and persist in configurations of much higher energy than equivalent 3D relaxations. Furthermore, axisymmetric calculations are expected to exhibit artificialities similar to the slab-symmetric ones studied here because they share the same topological simplifications. In particular, the axisymmetric magnetic surfaces, now toroidal shells that are labeled by their value of the azimuthal vector potential component A_ϕ, each enclose a conserved magnetic helicity.

Our results indicate that as a field's complexity increases, so does the discrepancy between the energy of its most relaxed state in 2D and 3D. So, axisymmetry may be appropriate when the field is near a stable force-free equilibrium, when little or no energy resides in higher-order radial or angular modes. But when these modes are populated, the imposed symmetry could make them artificially persistent.

Numerical studies of pulsar magnetospheres and magnetar flares are quite challenging, and much of the progress in this field has been obtained by restricting them to axisymmetry. Simulations using both FFE (Parfrey et al. 2012; Yu 2012; Parfrey et al. 2013) and PIC (Cerutti et al. 2015) show the development of a complex structure in the meridional plane. Similarly, a short-wavelength toroidal magnetic structure has been found in axisymmetric FFE simulations of black hole accretion. For example, the results of Parfrey et al. (2014) suggest that even disordered (as opposed to large-scale) magnetic fields advected inward by an accretion disk could facilitate angular momentum extraction from a black hole. The differences in 2D versus 3D found here suggest that extending these studies to be fully 3D, as that becomes computationally feasible, will be an exciting frontier and will likely reveal qualitatively new dynamics.

Another setting where 2D and 3D calculations are expected to differ is shock-generated turbulence, which may account for the relatively high magnetizations inferred from nonthermal emission spectra of astrophysical shock fronts—both in nonrelativistic settings such as supernova remnants and in relativistic settings such as gamma-ray burst afterglows. Amplification of the magnetic field by turbulence in and around the shock and its subsequent decay in the postshock flow have been studied extensively in both two (Mizuno et al. 2011) and three (Inoue et al. 2011) dimensions with relativistic MHD and also in 2D and 3D PIC simulations (Sironi & Spitkovsky 2009; Sironi et al. 2013). In this study we have observed a power-law decay of the magnetic energy in all settings, but with a steeper index in 3D than in 2D. This should be kept in mind because even minor differences in the decay law can have an impact on the efficiency of first-order Fermi processes (Lemoine et al. 2006; Niemiec & Ostrowski 2006; Niemiec et al. 2006; Pelletier et al. 2009), as well as on interpretations of gamma-ray burst afterglows (Gruzinov & Waxman 1999; Rossi & Rees 2003; Lemoine 2014).