The formation and dissipation of current sheets and shocks due to compressive waves in a stratified atmosphere containing a magnetic null

TL;DR

This study uses numerical models to analyze how magnetohydrodynamic waves propagate, form current sheets, and generate shocks in a stratified atmosphere with a magnetic null, revealing insights into solar atmospheric heating and spicule formation.

Contribution

It demonstrates the formation of current sheets and shocks due to wave refraction and mode conversion near a magnetic null, with implications for solar atmospheric heating mechanisms.

Findings

Null current sheets undergo multiple reconnection episodes.

Shock and Ohmic heating may account for 1-10% of coronal bright point radiative losses.

Wave dynamics are similar to mechanisms proposed for type-II spicule generation.

Abstract

We study the propagation and dissipation of magnetohydrodynamic waves in a set of numerical models that each include a solar--like stratified atmosphere and a magnetic field with a null point. All simulations have the same magnetic field configuration but different transition region heights. Compressive wave packets introduced in the photospheric portion of the simulations refract towards the null and collapse it into a current sheet, which then undergoes reconnection. The collapsed null forms a current sheet due to a strong magnetic pressure gradient caused by the inability of magnetic perturbations to cross the null. Although the null current sheet undergoes multiple reconnection episodes due to repeated reflections off the lower boundary, we find no evidence of oscillatory reconnection arising from the dynamics of the null itself. Wave mode conversion around the null generates a series of slow mode shocks localized near each separatrix. The shock strength is asymmetric across each separatrix, and subsequent shock damping therefore creates a tangential discontinuity across each separatrix, with long--lived current densities. A parameter study of the injected wave energy to reach the null confirms our previous WKB estimates. Finally, using current estimates of the photospheric acoustic power, we estimate that the shock and Ohmic heating we describe may account for $\approx1-10\%$ of the radiative losses from coronal bright points with similar topologies, and are similarly insufficient to account for losses from larger structures such as ephemeral regions. At the same time, the dynamics are comparable to proposed mechanisms for generating type--II spicules.

Figures

Tables

Table 1. Table 1: Simulation Grid Results: Fraction of injected energy to reach the null’s wave conversion layer, defined in terms of the equipartition radius r E subscript 𝑟 𝐸 r_{E} and scaleheight H E subscript 𝐻 𝐸 H_{E} , as a function of the injection location x 0 subscript 𝑥 0 x_{0} and the transition region height y t ​ r subscript 𝑦 𝑡 𝑟 y_{tr} .

	$x_0\left(\mathrm{Mm}\right):$	$-6.0$	$-5.0$	$-4.0$	$-3.0$	$-2.0$	$-1.5$	$-1.0$	$-0.5$	$0.0$	$0.5$	$1.0$	$1.5$	$2.0$	$3.0$
$y_{tr}\left(\mathrm{Mm}\right)$	$r_E,H_e\left(\mathrm{Mm}\right)$
3.00	0.77, 0.48							0.17
3.15	0.45, 0.24							0.09
3.30	0.28, 0.14	0.0005	0.0009	0.0006	0.008	0.02	0.03	0.04	0.03	0.0001	0.03	0.04	0.03	0.02	0.008
3.45	0.19, 0.09							0.02

Equations

$$v_{y}(x,y,t)=v_{d}\exp\Bigl{[}-\frac{(x-x_{d})^{2}}{2w_{x}^{2}}-\frac{(y-c_{s}t)^{2}}{2w_{y}^{2}}\Bigr{]}$$

$${\rho }_1={\rho }_0{v}_y/c_s\text{and}{ϵ}_1={\rho }_1(\gamma -1)\frac{{ϵ}_0}{{\rho }_0}.$$

$$\partial_{t}\Biggl{[}\frac{1}{2}\rho_{0}v_{1}^{2}+\frac{{P}_{1}^{2}}{2\rho_{0}c_{s}^{2}}+\frac{\lvert\mathbf{b}_{1}^{2}\rvert}{2}\Biggr{]}\\ +\nabla\cdot\Bigl{[}{P}_{1}\mathbf{v}_{1}+(\mathbf{B}_{0}\times\mathbf{v}_{1})\times\mathbf{b}_{1}\Bigr{]}=0.$$

$$\frac{d\mathbf{x}}{dt}={\mathbf{v}}_g\text{ and }\frac{d\mathbf{k}}{dt}=-|k|\nabla v_{\varphi }.$$

$$\mathcal{D}(\omega,\mathbf{k},\mathbf{x})=0=\Biggl{(}\omega_{A}^{2}-k^{2}v_{A}^{2}\cos^{2}\theta\Biggr{)}\times\\ \Biggl{(}\omega_{\pm}^{2}-k^{2}\Bigl{[}\frac{1}{2}(v_{A}^{2}+c_{s}^{2})\pm\frac{1}{2}\sqrt{v_{A}^{4}+c_{s}^{4}-2v_{A}^{2}c_{s}^{2}\cos 2\theta}\Bigr{]}\Biggr{)},$$

$$\rho \frac{d\mathbf{v}}{dt}=-\nabla P+\mathbf{j}\times \mathbf{B}+\rho \mathbf{g},$$

$$\mathbf{j}\times \mathbf{B}$$

$$=-{\nabla }_{\perp }\frac{B^2}{2}+B^2(\hat{\mathbf{b}}\cdot \nabla \hat{\mathbf{b}})$$

$$=-{\nabla }_{\perp }\frac{B^2}{2}+B^2\mathbf{t}\equiv {\mathbf{f}}_{mp}+{\mathbf{f}}_{mt}$$

$$\delta \rho (s;t)=\frac{\rho (s;t)-{\rho }_0(s)}{{\rho }_0(s)},$$

$$S(s)=\frac{1}{2}\Bigl{(}\tanh\frac{s_{c}-s}{s_{w}}+1\Bigr{)}\Bigl{(}\exp(a+{bs})+c\Bigr{)}.$$

$$\displaystyle W(t)=\int_{\mathcal{A}}\Bigl{(}E_{K}(t)+E_{A}(t)+E_{M}(t)\Bigr{)}d\mathcal{A},$$

$$\displaystyle W_{\text{flux}}(t)=\int_{0}^{t}\int_{\partial\mathcal{A}}\Bigl{(}\mathbf{F}_{A}(t^{\prime})+\mathbf{F}_{B}(t^{\prime})\Bigr{)}\cdot\hat{\mathbf{n}}dldt^{\prime}.$$

$$dx(ix)=dx_{0}+dx_{0}\frac{f_{x}}{2}\Bigl{[}2-\tanh\frac{ix-c_{1}}{w}+\tanh\frac{ix-c_{2}}{w}\Bigr{]}$$

$$dy(iy)=dy_{0}+dy_{0}f_{y}\Bigl{[}1+\tanh\frac{iy-c}{w}\Bigr{]},$$

Full text

††software: LARE2D v.2.11 (Arber et al., 2001)

1 Introduction

In our previous paper (Tarr et al., 2017, hereafter Paper I) we reported on a numerical simulation of magnetohydrodynamic (MHD) wave propagation through a model solar atmosphere. The study of MHD waves in stratified atmospheres stretches back many decades (see for example Ferraro & Plumpton (1958) and the extensive bibliography given in Paper I). Many observations have demonstrated the ubiquity of oscillatory phenomena in basically every feature of the solar atmosphere that researchers have decided to study: sunspot umbra and penumbra, pores, active region filaments, large prominences, quiet sun regions, extreme ultraviolet (EUV) loops above active regions, and the more extended corona (see the Introduction in Paper I for specific examples, or reviews by Nakariakov & Verwichte, 2005; Khomenko & Calvo Santamaria, 2013; Khomenko & Collados, 2015; Warmuth, 2015; Arregui et al., 2018). These observations indicate that a substantial amount of wave energy is present throughout the atmosphere and therefore this energy is likely important for solving the long–standing coronal and chromospheric heating problems. This general topic raises the important question: “How do the waves get to the locations where they are observed?” Convective motions in the photosphere are the obvious source, but early models found that convectively induced compressive waves should reflect from relatively low in the atmosphere, or even below the photosphere for lower frequency waves. An important step forward was the recognition that mode conversion—the transfer of wave energy between different types of MHD waves—allows MHD waves to propagate to locations that are inaccessible to their hydrodynamic cousins. This topic was the focus of Paper I.

A distinguishing feature of our previous work was that the magnetic field contained a strong inhomogeneity, owing to the presence of a single magnetic null point. At the null, the magnetic field is zero and has a discontinuous first derivative. Although MHD waves in the presence of nulls have been studied repeatedly in the past (see the Introduction and Discussion in Paper I, and reviews by McLaughlin et al., 2011; Pontin, 2012), few had done so in the context of a stratified atmosphere. The presence of the null turned out to greatly affect the propagation of waves—and the energy they carry—throughout the domain, compared to the case of a more smoothly varying field. The null created a location in the upper chromosphere where mode conversion could take place, in addition to the lower–lying conversion region typically considered, for example, in studies of mode conversion in sunspots (Rijs et al., 2016). Essentially, the null allows wave energy to tunnel through a region in which it would otherwise be evanescent, although this process happens in a roundabout way involving mode conversion.

Wave mode conversion (the exchange of energy between the fast, Alfvén, and slow MHD waves, though we focus here only on the fast and slow waves), is an important process in our simulations, and we have discussed it at length in Paper I. Here, we summarize several important concepts needed to understand the present work. Mode converstion may occur when the eigenvectors describing the MHD modes become nearly degenerate (Tracy et al., 2003). A necessary condition for this is that the sound and Alfvén speed are equal, $c_{s}=v_{A}$. At such locations, there is near equipartition between pressure and magnetic forces, which, for a $\gamma=5/3$ plasma, give plasma $\beta=2c_{s}^{2}/\gamma v_{A}^{2}=1.2$. Hence, regions where $\beta\approx 1$ allow for mode conversion. Magnetic null points have $\mathbf{B}=0$, meaning that (i) $v_{A}\rightarrow 0$ at nulls, (ii) a high–$\beta$ region must exist around a null, and therefore, (iii) any wave approaching a null is susceptible to mode conversion (except in the strictly $\beta=0$ limit). Further, this means the nature of an MHD wave may change considerably as it crosses an equipartition region and transitions from a high–$\beta$ to low–$\beta$ environment, or vice versa. At the same time, the properties of fast and slow waves are only truly distinct far from equipartition layers, which complicates the discussion of mode conversion and wave propagation near the equipartition layer, and around nulls in particular. See McLaughlin & Hood (2006) for an early study of, and attempt to quantify, mode conversion near a magnetic null. We refer the reader to Paper I, especially their §5.1 and §7.2, for further discussion on mode conversion near nulls and how the properties of fast and slow MHD waves change across an equipartition boundary.

In Paper I we found that the conversion process at the null coincided with two important dynamic effects: (i) the null point collapsed into a current sheet, and (ii) a set of propagating slow waves arose at the conversion region surrounding the null. The formation, oscillation, and dissipation of the current sheet suggested that the incoming wave induced magnetic reconnection at the null point. Indeed, this appears to be a general feature of the interaction of waves with magnetic nulls: waves refract toward and deposit energy at null points, thereby forming current sheets. In the other causal direction, whenever reconnection happens, the unbalancing of the Lorentz force generates MHD waves, and the newly launched waves may steepen to form shocks as they propagate away from the reconnection site, heating the plasma away from the null in the process. In this sense, the study of MHD waves and reconnection are not separate topics, and will often need to be treated simultaneously. In our previous work we only noted the existence and oscillation of the current sheets and slow mode waves. In the present work we will discuss how the slow modes form shocks, how the current sheet forms, and the ultimate fate of both.

In Paper I we did not investigate the physical process by which the current sheet arose and then oscillated. However, Murray et al. (2009) explored what appears to be a similar phenomenon in their 2.5D simulation of an emerging flux tube in a background field. A null point naturally formed between the emerging and surrounding flux, and as emergence drove the system, the null periodically collapsed into a current sheet and underwent reconnection. It is unclear if the current sheet forms in the same way in our wave driving experiment, and we answer the question of how our current sheet forms, below. The salient feature of Murray et al. (2009) is that the current sheet retracts because force balance across separatrices near the outflow causes the separatrices to pull apart.

We describe the present simulations in §2, including, in §2.1, how they differ from our previous work in Paper I, and the primary results in §2.2. Details for one representative simulation are presented in §3, which explains the null point collapse in terms of force balance in §3.1, quantifies the reconnection rate and compares it to Ohmic heating in §3.2, and explains slow shock formation and dissipation in §3.3. Part one of the parameter study, in §4, compares the long term properties for simulations with varying transitions region heights, $y_{tr}$. Part two of the parameter study, in §5, determines the mode conversion efficiency as a wavepacket is introduced at different locations throughout the simulation, and confirms the results of a WKB analysis given in Paper I. Section §6 applies our results to the question of radiative losses in coronal bright points, and compares our results with other relevant work. Finally, we briefly summarize our results in §7.

2 Simulation Overview

2.1 Setup and differences from Paper I

The simulations studied in this work are updated versions of the one described in Paper I; we refer the reader to §2 of that paper for a complete description, and here note only changes made to better isolate the behavior of the shocks. The full resistive MHD equations are solved using the LARE2D code (version 2.11; Arber et al., 2001) in terms of the primitive variables mass density, specific internal energy, plasma velocity, and magnetic field, $\rho,\epsilon,\mathbf{v},\mathbf{B}$, respectively. Our simulation consists of an initially magnetohydrostatic atmosphere that is subsequently driven by injecting a compressive wave packet through the lower boundary. The solid black lines with arrows in Figure 1 show the magnetic configuration near the null, and the solid contours show several cases of the plasma stratification, as described below. The null is located at $(x,y)=(0,3.75)\,\mathrm{Mm}$, and the four field lines connected to it form the separatrix dome. Compared to Paper I, we introduce three main changes to the simulation: two to the initial condition and one change to the wave driver.

First, we now use a stretched numerical grid that has high resolution in regions of interest near the null and lower resolution approaching the side and top boundaries. The high resolution region has uniform grid spacings $\Delta_{x}=\Delta_{y}\equiv\Delta=L_{N}/32\approx 4.7\,\mathrm{km}$ and spans $\lvert x\rvert\lesssim 9\,\mathrm{Mm}$ horizontally and $0<y\lesssim 6\,\mathrm{Mm}$ vertically. We use the subscript “N” to denote numerical normalization factors, which are unchanged from Paper I. The factor $L_{N}=150\,\mathrm{km}$ normalizes the numerical lengths, and is the density scale height at the model photosphere. The high resolution is a factor of four (in each direction) greater than that used in our previous simulation. The other numerical normalizations are $\rho_{N}=3.03\times 10^{-4}\,\mathrm{kg}\,\mathrm{m}^{-1}$, which normalizes the density, and $B_{N}=0.12\,\mathrm{T}$ which, normalizes the magnetic field. All remaining normalizations can be derived in terms of $B_{N},L_{N},\rho_{N}$. For quick reference, velocity is normalized to $v_{N}=B_{n}/\sqrt{\mu_{0}\rho_{N}}\approx 6.177\,\mathrm{km/s}$, time to $t_{N}=L_{N}/v_{N}\approx 24.28\,\mathrm{s}$, and current to $j_{N}=B_{N}/(\mu_{0}L_{N})\approx 0.63\,\mathrm{A}\,\mathrm{m}^{-2}$. Resistivity is normalized to $\mu_{0}L_{N}v_{N}$, so that $S=\mu_{0}L_{N}v_{N}/\eta$ is the Lundquist number defined using the pressure scale height at the lower boundary. The uniform resistivity is $\eta=0.0333\approx\Delta/L_{N}$, compared to a value of $0.1$ in Paper I111We followed the same guidelines in setting the uniform resistivity $\eta$ as we did in Paper I. See their §2.1 for a discussion of how numerical resistivity scales with the grid resolution and typical lengthscales of dynamic evolution.. For the remained of the paper, we have written equations in normalized units, yet restored physical units for the figures.

We stretch the grid in both the vertical and horizontal directions using $\tanh$ profiles as described in Appendix A. The uniform, high resolution region is described above. The lowest resolution regions (furthest from the dome) have $\Delta_{x}\approx 0.526L_{N}\approx 79\,\mathrm{km}$ and $\Delta_{y}\approx 0.47L_{n}\approx 70\,\mathrm{km}$. The stretched grid ensures that the side and top boundaries are far from the null point and injection locations, mitigating possible boundary effects, such as reflections. The side and top regions also include the damping term described in Paper I, which reduces reflections by removing kinetic energy near the boundaries. The full simulation has $6144\times 4096$ computation cells and spans $\lvert x\rvert\lesssim 52.5\,\mathrm{Mm}$ and $0<y<120\,\mathrm{Mm}$ in the horizontal and vertical directions, respectively.

For the second main difference, we modify the height of the transition region $y_{tr}$, varying from the value $y_{tr}=3.0\,\mathrm{Mm}$ used in Paper I (see their Equation 7) to $y_{tr}=3.45\,\mathrm{Mm}$ in steps of $0.15\,\mathrm{Mm}$. As the transition region height increases the stratification results in a lower (and more realistic) coronal density: the plasma temperature, and through it the density scale height, remains lower to a greater height above the photosphere, and therefore the density decreases. The primary simulation we focus on has $y_{tr}=3.3\,\mathrm{Mm}$; compared to Paper I, the density decreases an additional $\approx 2$ $e$–foldings before stabilizing to an isothermal corona with nearly uniform density.

In all simulations we have kept the null point at a height of $y_{null}=3.75\,\mathrm{Mm}$ so that the density is equal to (for $y_{tr}=3.0\,\mathrm{Mm}$) or lower at the null compared to the simulation in Paper I. For increasing $y_{tr}$, the sound speed is correspondingly lower near the null while the Alfvén speed is higher, so that the the equipartition layer, where $c_{s}=v_{A}$, is situated closer to the null, with an approximate radius of $r_{E}\approx 0.25\,\mathrm{Mm}$ for the $y_{tr}=3.3\,\mathrm{Mm}$ simulation. The solid green, magenta, red, and blue contours in Figure 1 show (for, respectively, the $y_{tr}=3.45,3.30,3.15,3.0\,\mathrm{Mm}$ simulations) the two equipartition regions in the simulation, one around the null and one running near the lower boundary (where all the curves nearly overlap). When used without qualification, “equipartition region” refers to the one surrounding the null. Figure 1 shows how the equipartition region shrinks as the transition region height increases and density at the null decreases. The magenta curve corresponds to the primary simulation in this work, with $y_{tr}=3.3\,\mathrm{Mm}$. The blue contour corresponds to the Paper I case, where $r_{E}\approx 0.75\,\mathrm{Mm}$. We will briefly compare the results of simulations using all four values $y_{tr}$ in §4.

For the third main difference, we modify the introduced “wave” to contain just a single Gaussian–shaped pulse (c.f. Paper I Equation 21). This was done to better study a single interaction of the pulse and the null point. When a wave interacts with a null it generates a complicated response that includes current sheet formation, mode conversion, and shock generation. If subsequent amplitude peaks in a wave train are introduced they further interact with each type of response, complicating the analysis. For this reason, in this paper, we have chosen to focus on the system’s response to just a single pulse.

We introduce a wave pulse using the method described in Paper I: we specify a spatio–temporal vertical velocity perturbation in the simulation’s lower–boundary ghost cells and use the velocity perturbation to define adiabatic energy and density perturbations. The velocity pulse is described in space and time by

[TABLE]

with associated specific internal energy and density terms

[TABLE]

Subscripts [math] and $1$ refer to background and perturbation quantities, respectively. For the primary simulation, the pulse has spatial parameters $x_{d}=1\,\mathrm{Mm}$, $w_{x}=1.977L_{N}=0.297\,\mathrm{Mm}$, and $w_{y}=0.15L_{N}=0.0225\,\mathrm{Mm}$, and is advected upward through the lower boundary at the local sound speed $c_{s}=\sqrt{\gamma P_{0}/\rho_{0}}=\sqrt{\gamma(\gamma-1)\epsilon_{0}(y=0)}$. The terms $x_{d}$ and $w_{x}$ are the pulse’s centroid and width in $x$. Because of the vertical advection at the sound speed, the vertical width $w_{y}$ sets the driving duration $T_{d}\approx w_{y}/c_{s}$. In the parameter study discussed in §5, we vary the wavepacket injection location $x_{d}$ between $-6.0\,\mathrm{Mm}$ and $3.0\,\mathrm{Mm}$ for the $y_{\text{tr}}=\,\mathrm{Mm}$ atmosphere.

To better isolate the waves from the background, we also run a baseline simulation, with the exact same initial condition, but no injected wave. The difference between the simulations with and without the wave isolates the wave–perturbed quantities $\rho_{1}$, $\mathbf{v}_{1}$, etc. Other methods of background subtraction, such as differencing the wave–injected simulation with itself at an initial time or a running average of previous times, all show nearly the same dynamics, but also include minor contributions from currents and velocities that arise due to secular drifts in the simulation caused by imperfect numerical representation of derivatives (see Tarr et al., 2017, §.2.2). These dynamics are unimportant compared to those of the injected wave packet: amplitudes in the baseline simulation are typically at least 2 orders of magnitude less than those in the wave injection simulation. The background subtraction we use primarily produces visualizations with fewer distracting features when the color scales are saturated to best visualize the dynamics near the null, but does not otherwise affect the analysis.

2.2 Primary results

Figure 2 shows a time series of the $y_{tr}=3.3\,\mathrm{Mm}$ simulation output in terms of the out–of–plane current density. This view highlights the essential elements of the numerical experiment: (i) the pulse, injected at $\mathbf{x}_{0}=(1,0)\,\mathrm{Mm}$, initially rises and expands ($t<4.5t_{N}$); (ii) a portion of the pulse refracts towards the null while the majority refracts away and back to the photosphere ($t\approx 4.8t_{N}$); (iii) the null point collapses into a current sheet ($t\approx 5.3t_{N}$); (iv) the current sheet eventually retracts and reforms oriented $90^{\circ}$ to the original sheet, and with opposite current density ($t\approx 5.8t_{N}$); (v) outward propagating disturbances from the null are largely confined to be near the separatrices ($t\geq 5.3t_{N}$); and finally, (vi) some additional oscillations at the null take place as the system seeks a (new) equilibrium state.

Other physical parameters (velocity, density perturbation) show essentially the same information as the current density in Figure 2. A more useful representation of the data for understanding the dynamics is shown in Figure 3. Here, each color channel represents a different energy density term in the wave conservation relation derived in Paper I,

[TABLE]

The three terms within the time derivative are the kinetic ($E_{K}$), acoustic $(E_{A})$, and magnetic ($E_{M}$) energy densities, and the two terms within the divergence are acoustic $(\mathbf{F}_{A})$ and magnetic ($\mathbf{F}_{M}$, i.e. Poynting) fluxes. Note that each term is a second order quantity that describes an energy density or energy flux of waves, excluding the background. The upper left panel contains three color wheels that show how the terms combine. Cyan (yellow) colors show equipartition between kinetic and magnetic (acoustic) energy densities. Therefore, outside the magenta contours in Figure 3, where $\beta<1.2$, cyan (yellow) represent fast (slow) waves, while inside the magenta contours ($\beta>1.2$) the reverse holds. The color scales are saturated to best represent the wave energy densities in the upper chromosphere to low corona, basically around the height of the null222The relations between the acoustic versus magnetic partition of wave potential energy, the type of MHD wave, and the local value of $\beta$ is described in detail in Paper I: see their §3 and §5.1..

Figure 3 and its animation (available in the online material) illustrate the dynamics of the simulation in terms of the partition of wave energy ($E_{K},\ E_{A},\ E_{M}$) at each location. Initially the pulse travels upward ($t=3.5-4.7t_{N}$) and then splits ($\approx 4.9t_{N}$) into multiple pieces. The sides refract back toward the photosphere while the central region refracts toward the null. All wave sections are predominately fast mode waves in the $\beta\lesssim 1$ regions, having near equipartition between magnetic and kinectic energies (cyan). The splitting of the initial wave packet into multiple subpackets was described in detail in Paper I §3 using a ray–tracing WKB method; we briefly summarize that method in our §3, below. See especially Paper I Figures 3-5, which show how different sections of an initial wave packet refract either toward or away from the null due to local gradients in the phase speed of the wave, and their §5.1-5.3, which discusses how the partition of wave energy changes for each section of the wave packet.

In the following work we focus only on the central wavepacket region that refracts toward the null: we will speak of the “central portion” of the pulse front and the “wings,” both referring to the just central region of the initial pulse, after the split. Figure 4 shows three example fast wave trajectories calculated as described in §3 that illustrate the three basic ways that the central portion of the wave packet interacts with the null.

The initial pulse steepens as it propagates, and eventually gives rise to multiple shocks. The plasma $\beta$ varies throughout the simulation, so the pulse shocks in distinct ways along different portions of the front. The central portion propagates directly toward the null, approximately perpendicular to the magnetic field, and causes the null to collapse into a current sheet. The solid line in Figure 4 illustrates this path, while the energy densities are seen in Figure 3.

The wings of the pulse refract and wrap around the null, eventually reaching locations where they propagate nearly parallel to the field. In Figure 4, this occurs where the dotted trajectory passes the equipartition contour. At those locations, the plasma displacement vector and wave speeds are such that mode conversion is an efficient process (Tracy et al., 2003). For an extensive discussion of the mode conversion process near the null and quantification of the energy exchanged between modes, see Paper I, especially their §5. In Figure 3, the mode conversion sequence can be seen by following the cyan regions near the equipartition contour at time $4.9t_{N}$, which convert to the magenta/white regions at $5.1t_{N}$, and finally the four strongest yellow regions at $5.3t_{N}$. These latter are the primary converted slow shocks that propagate away from the null along the separatrices. Behind each slow shock is a rarefaction wave followed by a reverse shock that brings the plasma nearly, but not quite, back to its original state.

Finally, at the far outer edges of the wave front (but still in the portion that refracts toward the null) the fast and slow wave speeds remain distinct so that the incoming fast waves do not immediately mode convert, but instead wrap around the null: see the dashed trajectory in Figure 4. Eventually, these parts of the wings also reach locations where mode conversion is an efficient process, at which point they mode convert to generate secondary, weaker sets of slow shocks on the opposite sides of each separatrix: see, for example, the two yellow regions above and leftward of the current sheet at $t=5.3t_{N}$.

Several more details are worth highlighting now. As the central region of the pulse impinges on the null around $t=4.9t_{N}$, the central portion concentrates in the region where the magnetic field orientation is roughly perpendicular to the pulse front (Figure 3). This fast pulse causes the null to collapse and form a current sheet oriented $+45^{\circ}$ counterclockwise (CCW) from the $\hat{\mathbf{x}}$ direction (Figure 2 at $5.3t_{N}$): the sheet is elongated parallel to the front. The sheet then reconnects, shrinks, and eventually reforms in the $-45^{\circ}$ direction around $t=5.7t_{N}$, as shown in Figure 5. At that level of detail, we see that the separatrices connect to the current sheet in a cusp shape ($\prec$) instead of a Petschek “Y” shape, as described by Vekstein & Priest (1992, 1993); Uzdensky & Kulsrud (1997); Kulsrud (2011). The cusp shape is required by the presence of currents along the separatrices, which connect to the current sheet of the collapsed null. This is true whenever the current sheet is present (e.g., also at $5.3t_{N}$, and during each later episode).

Subsequently, the sheet undergoes several more oscillations, damping all the while, as seen in Figure 2. These oscillations coincide with reconnection, a process that has been described in numerous other places, for instance, for a simulation of a flux tube emerging into a uniform vertical field (Murray et al., 2009), and for a linear null point (McLaughlin et al., 2012a). We discuss the oscillatory reconnection further in §3.2.

Returning once more to the wings of the pulse, the mode–converted slow waves subsequently shock and the newly formed slow mode shocks propagate away from the null along each separatrix. These shocks are distinct from the slow mode shocks of reconnection outflow, which do not propagate and instead remain attached to the collapsed–null current sheet.

The slow shocks lose energy as they propagate, heating the plasma in their wake. Because a given shock that forms due to mode conversion is localized to one side of the separatrix, so is the heating. As we show below, the shocks that form underneath the separatrix dome are stronger than those on the outside, so the heating is asymmetric across each separatrix. The heated plasma expands to restore force balance, thus pushing the magnetic field away from a potential configuration and generating currents localized along each separatrix in the process.

The final state of the system therefore includes: (i) a tangential discontinuity across each separatrix that divides the more strongly shocked plasma on one side from the less shocked plasma on the other; (ii) an hourglass shaped region of increased internal energy density at the null related to reconnection outflows and dissipation within the collapsed–null current sheet; and (iii) persistent current density at the null and along each separatrix, as required to balance the new pressure gradients. We detail each of these features in the following.

3 Formation of Shocks and current sheets

In this section, we provide an analysis of the current densities which localize to the separatrix surfaces, as reported in Paper I. The formation and evolution of these current sheets are interwoven with the formation and propagation of various shocks in the simulation: the two studies cannot be decoupled. A recent overview of the types of shocks that can arise in resistive MHD can be found in Goedbloed et al. (2010) Chapter 20.

Since the (acoustically dominated) fast wave packet is injected away from $x=0$, as was the case in Tarr et al. (2017), most of the initial wave energy efficiently transmits across the lower equipartition layer to become a magnetically dominated fast wave in the low $\beta$ portion of the simulation. As described above, a section of the fast wave refracts toward the null. Along the wave front’s trajectory, the fast mode wave speed first increases upward from the photosphere but then decreases again near the null. That latter decrease (even as the density drops) causes the wave to shock once it nears the null point; that is, it would not shock in the absence of a null, or at least a magnetic minimum. Identifying the type of shock is difficult in this case: the pulse has some properties of a fast shock, but the transient dynamics of reconnection and mode conversion confuse the analysis. We choose instead to analyze the forces near the null in the following. On the other hand, the waves leaving the null behave like quasi–steady–state slow mode shocks.

To characterize the evolution of the shocks we start with the generalized Wentzel–Kramers–Brillouin (WKB: Weinberg, 1962) approximation to find the trajectory of wave packets. We present just the necessary definitions for the present analysis and refer the reader to Paper I §3.1 and Appendix A for details. To lowest order, a shock will follow the same trajectory through an undisturbed background as a small–amplitude wave, but propagate at a faster speed (Uralova & Uralov, 1994). Determining the trajectory amounts to integrating Hamilton’s equations for a wave packet characterized by an initial location $\mathbf{x}(t=0)$ and wave vector $\mathbf{k}(t=0)$:

[TABLE]

Here, $\mathbf{v}_{g}=\partial_{\mathbf{k}}\omega$ is the group velocity, $v_{\phi}=\omega/\lvert k\rvert$ is the phase speed, and $\omega$ is the frequency of the specific wave type: Alfvén, fast, or slow. The phase and group velocities are calculated from the spatially–varying dispersion relation (see Paper I §3.1),

[TABLE]

where $v_{A}$ and $c_{s}$ are the Alfvén and sound speeds, and $\omega_{A},\ \omega_{+}$, and $\omega_{-}$ select the Alfvén, fast, and slow branches of the dispersion relation, and each has a positive and negative solution, for a total of six solutions. We measure distance along a trajectory by $s$. For each solution to $\mathcal{D}=0$, $\omega\propto k$ so the waves are dispersionless. However, the waves refract as they propagate due to the spatial dependence of the phase and group velocities. The magnitude of the wave vector divides out of Hamilton’s equation, so hereafter we only refer to the normalized wave vector; e.g. the normal direction to a phase front has $\lvert\mathbf{k}\rvert=1;\mathbf{k}=[\cos(\phi),\sin(\phi)]$ where $\phi$ is the angle measured counterclockwise from the positive $x$ direction. Ray trajectories are found by numerically integrating equations (4) using a Runge Kutta method.

3.1 Current Sheet formation

To study how the pulse causes the null to collapse into a current sheet, we first find a fast mode trajectory connecting the injected wave packet to the region around the null: Equations (4) are solved for the $\omega_{+}$ branch of Equation (5) with the initial conditions $\mathbf{x}(0)=(0.86,0.005)\,\mathrm{Mm},\ \mathbf{k}(0)=(0,1)$; this trajectory (Figure 1), passes within a few computational cells of the initial null point at a distance $s=3.94\,\mathrm{Mm}$ along the curve. The center of the initial current sheet is shown at time $t=4.95t_{N}$ in Figure 6. Plasma parameters extracted along this curve for all output simulation times are analyzed below. The long axis of the current sheet is perpendicular to the incoming fast ray trajectory, so the current sheet forms in the plane of the front, $\approx 45^{\circ}$ CCW from $+\hat{\mathbf{x}}$ direction.

Figure 7 shows the magnetic pressure $\frac{1}{2}B^{2}$ (black) and plasma pressure $P=(\gamma-1)\rho\epsilon$ (red) in the vicinity of the current sheet along the fast trajectory. Solid curves are for time $4.95t_{N}$ and dashed curves for the background state. The abscissa shows distance $s$ along the trajectory, oriented to agree with the direction of propagation in the simulation: if time were to progress, the pulse would move from right to left. Derivatives of these curves are the magnetic and plasma pressure forces projected along the trajectory, which we discuss in depth below.

As the current sheet forms the magnetic field becomes more distorted and enhanced on the incoming side of the null (rightward in Figure 7) compared to the outgoing side, indicating that the magnetic part of the pulse is not able to propagate across the null. However, the pressure increases across the null, indicating that the pressure dominated part of the wave is, in contrast, able to propagate across the null. From this analysis we found that the initial incoming fast wave is split into a magnetically dominated and an acoustically dominated pulse, each with separate dynamics. The magnetic portion of the pulse has converted to a slow–type disturbance and is unable to cross the null, while the energy and density disturbance (the plasma pressure) transmits across the null as a fast wave, but carries substantially less energy. This is consistent with the dynamics of wave splitting and propagation determined in Paper I by numerically integrating the wave conservation relation (3) about the null. The result, then, is that the current sheet forms due to a strong gradient in the magnetic pressure as flux piles up and is unable to rapidly leave the null region.

Because $\beta>1$ inside the equipartition region, we identify the pressure pulse that is able to cross the collapsing null as a transmitted fast wave, while we identify the magnetic pulse that stagnates and forms the current sheet as a slow wave. The slow wave phase propagation vector $\mathbf{k}$ is perpendicular to the field, and so the slow phase velocity becomes 0. But recall that the slow wave’s group velocity limits to the cusp velocity $v_{c}=v_{A}c_{s}/\sqrt{v_{A}^{2}+c_{s}^{2}}$ directed along the magnetic field in this case (Goedbloed & Poedts, 2004, §5.3.2) so the newly formed slow wave packets propagate outward along field lines connected to the current sheet. Meanwhile, the sheet’s center moves along the fast mode trajectory as it forms, decelerating as it does so. This is because the transverse forces on the sheet are unbalanced. The sheet’s length and maximum value of current density increase until $t\approx 5.04t_{N}$, after which the sheet begins to retract, weaken, and migrate backwards along the same trajectory, reconnecting all the while.

We can get a better handle on the current sheet formation by considering the forces in the momentum equation, written in the normalized units used in the code as

[TABLE]

where $\frac{d}{dt}=\partial_{t}+\mathbf{v}\cdot\nabla$ is the advective derivative. Near the null, the gravitational force is several orders of magnitude smaller than the pressure and magnetic forces caused by gradients in background or the wave, so we will mostly ignore that term. The Lorentz force $\mathbf{j}\times\mathbf{B}$ may be decomposed into two terms corresponding to a pressure and a tension,

[TABLE]

where the curvature vector $\mathbf{t}$ defines the radius of curvature $R_{c}$ through $\lvert\mathbf{t}\rvert=\frac{1}{R_{c}}$, and $\mathbf{B}=B\hat{\mathbf{b}}$. Note that the parallel components of the derivative (along $\hat{\mathbf{b}}$) cancel between the two terms on the right side of (7), so that the remaining terms are perpendicular to the field line: in the 2D case, they act in the same (or opposite) direction.

By design, the initial condition is in force balance. Figure 8(a) shows several examples of the curvature magnitude (red circles) and direction (blue lines), while Figure 8(b) shows the magnetic pressure (solid) and tension (dotted) force density terms calculated separately, as well as their sum (dash–dot), along a $-45^{\circ}$ line passing through the null. Near the null, the derivative of the plasma pressure along this line is initially nearly zero. The derivative of the perpendicular magnetic pressure along the line smoothly transitions from negative below the null to positive above it, while the tension does the opposite. Although $\beta>1$ around the null, the Lorentz force is still initially self–balanced.

Figure 9 shows similar force diagrams for time $4.95t_{N}$, where each force term is calculated along the solid fast ray trajectory from Figure 4, reproduced here as the red curve. The grayscale in panel (a) is the current density $j_{z}$. Panel (b) shows that the introduced wave pulse breaks the balance between the magnetic tension and magnetic pressure gradient forces, as easily seen at distances $<1\,\mathrm{Mm}$.

A detail of Figure 9 near the null point is shown in Figure 10. Panel (b) now also includes the plasma pressure gradient force (dashed). A time series of the force balance near the null is shown in Figure 11 and the supplemental animation available in the online material. The time series includes additional force terms and their combinations: the plasma pressure gradient force (dashed line), the gravitational force (solid red), and the total force (dash-dot-dot-dot).

Figure 10 highlights the balance of forces as the initial current sheet forms at time $4.95t_{N}$. At this time, the initial pulse has split in two: there is the main pulse associated with the current sheet located at a distance $s=3.94\,\mathrm{Mm}$ along the trajectory, and a transmitted pulse around $s=4.2\,\mathrm{Mm}$. The magnetic pressure gradient force and the plasma pressure gradient force are the derivatives of the black and red curves from Figure 7, respectively.

Panel (b) shows that the tension force (dotted line) smoothly varies through zero across the collapsing current sheet, marked by the vertical line labeled $j_{max}$: the tension force does not cause the localization of currents. Conversely, the magnetic pressure gradient force (solid) shows multiple local minima and maxima, with each extrema associated with extrema of the current density (compare to panel (a)). The sum of the two magnetic terms, the total Lorentz force (dash–dot), shows that the magnetic pressure dominates the tension inside the pulses (e.g., at locations $\approx 3.9$ and $4.2\,\mathrm{Mm}$). At the same time, a plasma pressure gradient (dashed) has built up to nearly counteract the magnetic pressure gradient.

As we now detail, the decoupling of the plasma pressure force from the magnetic pressure force at the null is what allows the current sheet to develop. From the time series in Figure 11 we see that as the pulse approaches the null, the magnetic pressure gradient (solid) increases in magnitude relative to the background at the pulse’s leading edge, and it decreases in magnitude behind the front. This is most easily seen in the change in the solid line to the right of the null $(s\approx 3.8\,\mathrm{Mm})$ between times $4.825$ to $4.9t_{N}$. During this time, the tension term (dotted) simply increases in magnitude everywhere behind front, though this is a subtle effect: it is most easily seen from the leftward drift over time of the dotted line along the lower boundary of each panel.

The total Lorentz force therefore has a localized leading force directed towards the null, due mostly to magnetic pressure, and a fairly uniform trailing restoring force away from the null, due mostly to tension. This situation holds until a reflected pulse forms and the current sheet begins to retract at time $\approx 5.1t_{N}$, as seen in the online animation.

The dynamics of the pulse rapidly changes once it crosses the equipartition region. When the magnetic pressure pulse crosses the equipartition region, the plasma compression associated with it begins to cause an energetically important increase in plasma pressure: see the orange dashed line in Figure 11 starting at time $4.775t_{N}$. The leading and trailing edges of the plasma compression give rise to pressure gradient forces oriented in the forward and backward directions along the fast trajectory (dashed line at $\approx 3.7\,\mathrm{Mm},4.8t_{N}$). As the pulse reaches the location where the magnetic field switches direction (the null), two things happen at once: (i) the forward plasma pressure pulse is able to cross the null as a (still high–$\beta$) transmitted fast wave and becomes the dominant force on the other side of the null: this can be seen by comparing the total force (dash–dot-dot-dot) to just the magnetic force (dash–dot) at $t=4.9t_{N}$; (ii) the magnetic pressure pulse stalls on the incoming side of the null, and is mostly balanced by the trailing plasma pressure gradient pulse: compare the dashed and solid lines near $x=3.8\,\mathrm{Mm}$.

Once the transmitted fast wave crosses the equipartition region on the exiting side (leftward) of the null at $4.95t_{N}$, the plasma pressure gradient begins to lose its dominance while the magnetic pressure again starts to dominate the pulse dynamics: this pulse remains a transmitted fast wave, now in the low–$\beta$ region. The online animation of this figure shows that this transition is complete by $5.05t_{N}$, when the transmitted fast wave is at $s=4.6\,\mathrm{Mm}$. According to our calculations in Paper I, the transmitted fast wave carries away $\approx 7\%$ of the energy to impinge upon the null.

3.2 Reconnection and null point heating

The strong currents that localize to the null lead to both reconnection and Ohmic heating near the null. Reconnection involves an extended area around the null where field lines are swept into the current sheet. Because the field strength drops with height, these field lines map back to a narrow region at the lower boundary which straddles the separatrix near $x=0$. To study the reconnection, we traced 1000 field lines from $y=0$ spanning $x=\pm 0.47\,\mathrm{km}$ (one tenth of a simulation voxel to either side of $x=0$) and tracked through time the locations of their opposite footpoints. During reconnection, these jump discontinuously between the locations of the other two photospheric separatrices at $x=\pm 3.75\,\mathrm{Mm}$.

Figure 12 describes the reconnection over the course of the simulation. Panel (a) shows $L(x)$, the position of the central separatrix along the lower boundary. This separatrix moves laterally as reconnection transfers field lines through the current sheet. Several reconnection episodes are evident, but the secular trend is to reconnect field lines between the right domain inside the dome and the left domain outside the dome, which shifts the separatrix in the $+\hat{\mathbf{x}}$ direction. This is consistent with reconnection across the initial current sheet. Note that the lower boundary is line–tied with $v_{x}(x,0)=0$, so the shift of the separatrix is due to a change of topology.

The rate at which flux $\Phi$ is swept out as the separatrix location changes with time gives the reconnection rate $\dot{\Phi}$, which is shown in black in Figure 12(b); the finite number of traced field lines cause the discrete steps. The reconnection rate should be related to the current density at the null, where reconnection occurs, by $\dot{\Phi}=\eta j_{z}$, and $\eta=0.0333$ is the uniform simulation resistivity. We determine the location of the null by finding the minimum of $B^{2}$ using a Newton–Raphson method, and then perform a bilinear interpolation of the current density $j_{z}$ to that point from the closest grid locations.

The red curve in Figure 12(b) shows the null current density as a function of time, given in the normalized units $j_{N}$ on the right axis. The blue and green curves show the maximum amplitude current density in the current sheet and the current density at $y=4.7\,\mathrm{km}$, one gridpoint above the lower boundary along the $x=0$ separatrix. The null location is typically (but not always) the location of maximum magnitude current density in the current sheet. The scaling between the right and left axes is indeed $\eta$. Comparing the current density to the reconnection rate, we see that there appear to be two distinct processes at work: (i) one where the current density at the null scales with the reconnection rate like $\eta$ (e.g., times $\approx 5t_{N},\ 12t_{N},18t_{N}$), and a (ii) second process where the reconnection rate scales (rather roughly) with the current density near the lower boundary ($\approx 8t_{N},\ 15t_{N}$). The first set appear to be genuine reconnection at the null point. The second set appear to be diffusion caused by strong gradients in the current density near the lower boundary just as the fast wave reflects off that boundary. The side and top boundaries are far removed from the null so that, combined with the damping terms near these simulation boundaries, no reflections from these locations impinge upon the null.

The strong currents at the null in process (i) occur every time a fast wavefront sweeps across the null. The first of these ($5t_{N}$) is due to the initial wavepacket reaching the null. Note from Figure 12 that this episode has two reconnection phases, where reconnection happens first in one direction (peak at $5t_{N}$) and is immediately followed by reconnection in the opposite sense (peak at $5.5t_{N}$). There does not appear to be further oscillatory reconnection associated with this episode. This is also the case for the other two process (i) episodes at $12t_{N}$ and $18t_{N}$. These latter two episodes correspond to the fast wave traveling back from the null to the lower boundary, reflecting, then finally refracting back to the null point. Again, they do not seem to be a form of oscillatory reconnection inherent to the dynamics of the null. Each time process (i) occurs it generates a substantial reconnection episode, influencing a large area around the null. This implies that the processes of fast wave refraction into the null region, flux pile–up, and mode conversion are important for Ohmic heating within the null region.

In contrast, process (ii) does not coincide with strong currents at or around the null point, meaning its contribution to null point heating is minimal. Instead, it appears to be due to diffusion of magnetic field caused by strong currents near the lower boundary, cotemporal with the fast wave reflection off the lower boundary. It is an integral effect over the separatrix field line, leading to the lesser correlation between the green and black curves in Figure 12(b) compared to the red, blue, and black curves. Its amplitude also diminishes rapidly, as less of the fast wave energy returns to the lower central separatrix. The similarity in the amplitudes between process (i) and process (ii) is due to the uniform resistivity, which was set to just exceed the numerical resistivity for this simulation.

The grayscale in Figure 13 depicts the change in temperature near the null point at the final simulation time, $t=20t_{N}$, relative to the baseline simulation where no wave packet was injected. A dumbbell shaped region of enhanced temperature is located near the null, $\approx 0.3\,\mathrm{MK}$ or $40\%$ above the temperature of the background stratified atmosphere. The plasma temperature along each separatrix is enhanced by $\approx 3-5\%$. The separatrix dome is highlighted in green, and red dashed curves mark the extent of field lines underneath the dome that were involved in reconnection, as determined in Figure 12(a). The field line associated with the initial separatrix surface is shown as the dashed blue line and demonstrates the direction of flux transfer. The plasma on opposite sides of the initial current sheet has been heated, and this heating is bounded by the field lines involved in reconnection. The distribution of plasma density shows essentially the same features as temperature distribution, but in the opposite sense: the final density is reduced by $\approx 30\%$ relative to the background within the dumbbell region, and by $\approx 1-2\%$ along the separatrices.

To further determine if the reconnection is associated with heating around the null we integrated the Ohmic heating term $H_{O}=\eta j^{2}$ in time for a set of traced plasma elements. Because we did not include conduction in our model, the heating integrates with each fluid element over time. We therefore calculated the integrated displacement of a $(51\times 51)$ grid of points surrounding the initial location of the null: $(-0.094<x<0.140,3.66<y<3.89)\,\mathrm{Mm}$. The location of each of point as a function of time is given by $\mathbf{r}(t)=\mathbf{r}(0)+\int_{0}^{t}\mathbf{v}(t^{\prime})dt^{\prime}$. This method allows only a rough comparison because we used the simulation output at a cadence of $dt\approx 150\delta t$, where $\delta t$ is the fundamental simulation timestep given by the Courant condition. At each time we found the current density at the displaced locations by bilinear interpolation from the closest Eulerian grid points. We next mapped the total Ohmic heating $\int H_{O}(t)dt$ at each Lagrangian point over time back onto the Eulerian simulation grid at the final time, and compared the integrated Ohmic heating with the final temperature of the plasma. Contours of the Ohmic heating are displayed as solid blue lines in Figure 13. The temperature and Ohmic heating distributions do overlap, which indicates that Ohmic heating within the current sheet is responsible for a stable, dumbbell–shaped high energy region around the null at late times.

Interestingly, the dumbbell shaped enhanced energy region is oriented perpendicular to the plane of the original current sheet, which forms parallel to the incoming wave front. Comparison between the Lagrangian traces, the plasma density, and the current density reveals why. Reconnection outflows from the ends of first the current sheet cause a mass density depletion within the sheet. Most of the heating takes place inside the first current sheet while it is at its maximum extent, but this happens after the density depletion. Ohmic heating is thus able to heat the remaining plasma to high temperatures so that the overall pressure near the null is enhanced relative to the initial state. As the high pressure, high magnetic energy current sheet collapses and begins to reconnect in the opposite sense, the previous strongly heated regions get pushed out along the original inflow directions. The result is the high temperature, low density, dumbbell shaped region oriented perpendicular to the original wave front.

After the first reconnection episode ends around $6.0t_{N}$, the plasma density remains low and the temperature high relative to the initial state. In each later reconnection episode, the plasma temperature continues to increase by Ohmic heating so the pressure becomes locally enhanced and has strong gradients. In between each episode the plasma relaxes until the pressure has minimal gradients.

3.3 Slow shocks along the separatrices

In this section we study the mode–converted slow waves that form along each separatrix. As we will see, these waves eventually form slow mode shocks. We study the area around the left separatrix in detail, but the same analysis holds around all four separatrices. Note that, in particular, it holds around the upward separatrix that funnels wave energy higher into the corona. This brings to mind the observations of spicules, which we discuss briefly in §6.

The dashed curve in Figure 1 labeled “Slow Ray” shows a solution to Equations (4) for the slow mode branch of Equation (5) with initial conditions $\mathbf{x}(0)=(-0.563,3.436)\,\mathrm{Mm},\ \mathbf{k}(0)=(0,1)$: it traces out the trajectory of a mode–converted slow mode wave. The initial condition was selected to describe the leading edge of the wave pulse in the mode conversion region around the time when the initial pulse converts, at $t=4.925t_{N}$. The intersection of the wave front with the region of substantial mode conversion creates an extended region where the slow shocks form, all propagating on similar trajectories, essentially creating a shock tube. The case we analyze is roughly in the middle of this tube.

The shock itself travels faster than a slow wave (it is a shock) but follows the same trajectory. Afanasyev & Uralov (2012) discuss a ray tracing method that captures this behavior, where the only essential difference is the replacement of the phase velocity $a$ by $(a+\frac{1}{2}\kappa U_{sh})$ in the $\hat{\mathbf{k}}$ term of the ray equation, where $U_{sh}$ is the shock velocity, $\kappa$ is a numerical coefficient dependent on local plasma properties and propagation angle to the magnetic field, and the shock velocity depends on the shape of the waveform that generated the shock. We do not perform that analysis here, but simply determine the slow trajectory and then extract all simulation output along that trajectory at every snapshot in time.

We have performed the extraction along many trajectories, and as expected the plasma properties do not vary significantly for nearby trajectories within the shock tube. The shock strength varies from path to path, decreasing away from the separatrix. Inspecting the various plasma parameters along the slow trajectory from Figure 1 as a function of time shows that the shock fully forms after traversing $\approx 1\,\mathrm{Mm}$ along this path: the $*$ symbols in Figure 1 are spaced at $1\,\mathrm{Mm}$ intervals.

To define important features of the shock we introduce the fractional density change along the shock path relative to the background,

[TABLE]

where $s$ is the distance along the shock trajectory. Figure 14 shows an example of $\delta\rho$ at time $t=6.0t_{N}$. The form of the wave pulse resembles a standard N–wave, with a forward shock front, rarefaction wave, and a reverse shock at the tail (Courant & Friedrichs, 1948). Both the rarefaction wave and the tail, however, display more complicated behavior than for the strictly hydrodynamic N–wave (Ro & Matzner, 2017), pointing to the unique behavior of MHD shocks and the possible influence of reconnection outflows on the tail.

In order to automatically and robustly track the shock evolution we require a consistent definition of the preshock value and the shock peak. To do so, we fit the shock front with a combination of a hyperbolic tangent and an exponential decay behind the peak. The resulting functional form is

[TABLE]

The coefficients $\{s_{c},s_{w},a,b,c\}$ are determined using lmfit.pro, IDL’s built–in implementation of the Levenberg–Marquardt algorithm. The first two coefficients set the location and width of the shock front, while the latter three set the shock strength, decay rate, and depth of the rarefaction wave behind the shock, respectively. We define the preshock location as the point where the normalized $\tanh$ function reaches $0.001$ of its maximum and the postshock location (the shock peak) as the zero of the second derivative of the fitted profile; the former is calculated directly from the analytic function while the latter is determined numerically using a Newton–Raphson method.

The shock’s fractional density profile (black), analytic fit to the shock front (green), and locations of the preshock (black vertical) and postshock (blue vertical) locations are shown in Figure 14 at time $t=6.0t_{N}$. For reference, the red vertical line in Figure 14 demarks the location where the shock peak would be if it traveled at the slow mode speed after fully forming at $s=1\,\mathrm{Mm}$ at time $t=4.925t_{N}$, when the fully formed shock peak passed that location.

Figure 15 shows the temporal evolution of the fractional density (top) and current density (bottom) profiles along the slow shock trajectory; the color scale indicates time. The amplitude of the shock increases until $t\approx 6t_{N}$ due to the convergence of the surrounding slow mode trajectories, which are guided by the convergence of field lines away from the null. After that time, even though the field lines and slow trajectories continue to converge, shock dissipation becomes the dominant effect, and the shock amplitude decreases. As the shock dissipates and its amplitude decreases, so does its propagation speed. After $t=6.5t_{N}$ it essentially travels at the slow wave speed (the distance between the red and blue vertical lines in the time evolution of Figure 14 remains fixed). A sharp spike in current density forms at the leading edge of the shock (see Figure 15(b) at $t\approx 6t_{N},\ s\approx 1.75\,\mathrm{Mm}$), whose origin is currently unexplained.

Figure 16 compares the fractional pressure change along two slow mode trajectories, $\delta P$ defined analogously to $\delta\rho$ in Equation (10). Panel (a) shows the wave energy densities using the same scaling as in Figure 3, at $t=5.7t_{N}$, when disturbances are present on either side of the separatrix. Slow Path A, which is the same as that in Figure 1, traces one of the primary mode-converted disturbances, while Slow Path B follows an energy pathway that has first circled over the top of null, similar to the dashed path in Figure 4. More energy follows path A than path B. Because all of this energy eventually dissipates along each trajectory (and all corresponding trajectories near each of the other separatrices), and there is more energy along path A than path B, the regions of the simulation that experience the first and strongest mode conversions are ultimately heated more. Figure 16a also shows the asymmetry across separatrix at $x=0\,\mathrm{Mm}$.

The strongest shock heating is thus localized to one side of each separatrix. The resulting pressure imbalance causes the more heated side to expand and the less heated side to contract until an equilibrium is reached between the pressure gradient and Lorentz forces. The initial condition had a zero current (potential) magnetic field, so the new equilibrium contains currents, localized to each separatrix surface. Thus, a tangential discontinuity forms across each separatrix. The next two sections detail the consequences of asymmetric heating across the dome.

4 Long term properties

In order to better understand the interplay between current sheet formation, shock heating, and the size of the equipartition contour, we repeated a version of the simulation using the four different values for the transition region height: $y_{tr}=3.0,\ 3.15,\ 3.3$, and $3.45\,\mathrm{Mm}$. These values generate the four different stratification profiles (and corresponding changes in the size of the equipartition region) shown in Figure 1. For each simulation, we launched a high amplitude ($0.5v_{N}$) wave packet from $x_{0}=-1\,\mathrm{Mm}$, then allowed the system to relax. In each case, we find that the temperature increased throughout the flux domains lying underneath the separatrix dome. However, the enhanced temperature region becomes more strongly localized to the separatrix as the equipartition region shrinks.

Figure 17 shows a map of the pressure variance $\delta P(x,y)$. Continuing oscillations are visible in each simulation, and the phase of these oscillations varies from map to map because of the strong dependence of the Alfvén speed on $y_{tr}$. These oscillations are mostly visible outside the dome.

In contrast, the variations in plasma pressure surrounding the separatrix are basically stationary long after the pulse is introduced. As the height of the transition region increases the pressure increase becomes more strongly localized to the separatrix, as seen in the line plots of Figure 18. This is due to the shrinking of the equipartition region as the transition region height is increased, with the narrower conversion region producing more tightly localized shocks near the separatrix. Additionally, in the $y_{tr}=3.0\,\mathrm{Mm}$ case, the conversion primarily produces simple waves instead of shocks because of shallower gradients throughout the simulation, leading to comparatively less shock heating.

In all cases, energy conversion is asymmetric across each separatrix, as in the discussion of Figure 16. The asymmetric shock dissipation across the separatrices gives the plasma a higher pressure, due mostly to enhanced internal energy, on inside of the dome. The pressure gradient across the dome is balanced by a Lorentz force, which accounts for the enhanced current density near the separatrices. Thus, both the magnetic and plasma properties change discontinuously across the separatrices: these are tangential discontinuities, which are allowed normal to topological boundaries (see, e.g., Goedbloed et al., 2010, §20.2.2).

5 Parameter study

\floattable

In the previous section we demonstrated how the development of shocks and current sheets depend on the atmospheric parameters of the system. To study this topic further, we ran a set of simulations to form a parameter study on two variables: the height of the transition region $y_{tr}$ and the centroid location of an injected wave packet, $x_{0}$. The initial condition for each of these new simulations is the near–equilibrium state reached in the four simulations with varying $y_{tr}$ described in the previous section. To each of these simulations we inject a new wave packet with the properties described by equations (1) and (2). We perform 14 simulations for the $y_{tr}=3.3\,\mathrm{Mm}$ case, varying the injection location $x_{0}$ between $-6$ and $3\,\mathrm{Mm}$, and 3 additional simulations, injecting the wave packet at $x_{0}=-1\,\mathrm{Mm}$ for the cases $y_{tr}=3.0,\ 3.15,$ and $3.45\,\mathrm{Mm}$.

For each simulation, we measured the magnetic and acoustic fluxes, defined in Equation (3), entering and exiting the null region, by repeating the analysis given in §5.3 of Paper I, e.g., using their Equations (23) for $W(t)$ and (24) for $W_{\text{flux}}(t)$, which are reproduced below. The total wave energy inside a region $\mathcal{A}$ is given by the spatial integral of the wave energy density over the region at each time,

[TABLE]

while total flow of energy into or out of the region is given by the integral of the acoustic and magnetic energy fluxes across the boundary $\partial\mathcal{A}$:

[TABLE]

As in Paper I, we take $\mathcal{A}$ to be the area within one equipartition scaleheight of the equipartition layer. The equipartition layers are approximately circular with average radii $r_{E}$. The equipartition scaleheight, defined in Schunker & Cally (2006), is $H_{E}^{-1}\equiv\partial_{\tau}((c_{s}^{2}/v_{A}^{2})\rvert_{c_{s}=v_{A}})$ where we use the boundary normal for the derivative direction $\tau$ and then average over the $c_{s}=v_{A}$ surface. The values for $r_{E}$ and $H_{E}$ for each value of $y_{tr}$ are given in Table 1.

The total energy of the injected wave packet, $E_{\text{input}}$, is known from the simulation setup. The total energy density flux through the equipartition region around the null is then used to measure the fraction of the initial injected wave energy that makes it to the null. The ratio of the net acoustic flux to net Poynting flux which passes through the equipartition contour surrounding the null measures the efficiency of mode conversion in that region, thus giving the fraction of the initial injected energy that ends up as the slow mode shocks along the separatrices. The results for this analysis for the 17 simulations are summarized in Table 1 and plotted in Figures 19 and 20. We found the mode–conversion efficiency varied only by $\approx 20\%$ across all simulations, much less than the other factors, and do not consider it further.

Figure 19 shows that the fraction of injected wave energy to reach the null decreases as the height of the transition region increases. This quantifies the discussion in §4. The result for $y_{tr}=3.0\,\mathrm{Mm}$, with $E_{\text{null}}/E_{\text{{input}}}\approxeq 0.17$, repeats the result of Paper I333The fraction given here of $0.17$ differs slightly from the $0.155$ reported in the previous work, due to the slightly different background conditions.. A ratio of $0.04$ was found for the primary simulation discussed in the present work, where $y_{tr}=3.3\,\mathrm{Mm}$.

Figure 20 shows the effect of varying the wave packet injection location $x_{0}$ for the $y_{tr}=3.3\,\mathrm{Mm}$ simulations. These results are shown as asterisks connected by solid black lines. In Paper I §6 we calculated the effect of varying the wave packet injection location using a WKB method, (see Figure 11, dot-dashed line of that paper). The WKB estimate agrees with the simulation results when scaled by a factor $\alpha\approx 0.12$. That factor, which minimizes the square–difference between the simulation results and the WKB estimate, is expected: it accounts for both the result of Paper I, that the WKB method overestimates the amount of energy to reach the null by a factor of 1.8, and for the results of Figure 19, that a factor 4.25 less energy makes it to the null for $y_{tr}=3.3\,\mathrm{Mm}$ compared to the $y_{tr}=3.0\,\mathrm{Mm}$ case.

The strength of the outgoing shocks depend on the size of the equipartition region near the null and the amount of the wave’s initial energy that makes it to that region (in addition to the initial injection amplitude). In our simulations, as the height of the transition region increases and the equipartition region shrinks, the amount of energy that makes it to the null also decreases. At the same time, the null collapse and conversion process itself becomes more violent and localized, leading to stronger shock formation. This results in the localization of the pressure disturbance seen in Figures 17 and 18.

6 Discussion

The results of the previous two sections have important implications for heating regions of the Sun that share this type of dome topology, especially coronal bright points (CBP: Golub et al., 1974; Alipour & Safari, 2015). CBPs are typically (though not always) associated with a magnetic null and dome topology (Galsgaard et al., 2017). From our study, it is clear that fast mode waves introduced at different locations throughout a dome (and even slightly outside a dome) will refract towards the null, undergo mode conversion, and form slow shocks localized to the separatrices. If we assume that the incident wavefield is stochastic—i.e., not a single wave packet localized to one lobe of the magnetic dome, but a continuous, uniform spectrum across the entire lower boundary, then the average of the curve in Figure 20 is the expected energy to be pumped from the wave field into shocks and Ohmic dissipation, and ultimately heating the plasma in the dome. That value is 0.016, meaning that of order $1\%$ of the high frequency444(Say, a few times the Brünt–Väisälä or acoustic cutoff frequency, both of which are $\approx 5\,\mathrm{mHz}$ in the photosphere or low chromosphere Mihalas & Toomre, 1981; Vigeesh et al., 2017; Cally, 2006). convective wave flux through the photosphere could heat the solar plasma in these regions.

It is interesting to consider whether the processes described in this paper may help explain the association of CBPs with magnetic dome topologies. If the magnetic concentrations associated with them do not substantially suppress convection, then the convective wave field will continually pump energy into the plasma near the separatrices. For a survey of 285 CBPs, Longcope et al. (2001) determined their radiated power using SOHO EIT and MDI data (Delaboudinière et al., 1995; Scherrer et al., 1995). From the distribution of radiated power versus CBP size they found a typical heat flux over the CBP area of $1000\,\mathrm{W}\,\mathrm{m}^{-2}$. A typical radius for these structures is $5\,\mathrm{Mm}$, similar to the size of our simulated dome. Bello González et al. (2010) found that the total upward acoustic flux through the photosphere to be $\approx 7000\,\mathrm{W}\,\mathrm{m}^{-2}$; this is roughly double previously reported values, owing to better spatial and temporal coverage. Those previous estimates (Bello González et al., 2009) further found that $\approx 1/3$ of the total power is contained in the high frequency component $(>10\,\mathrm{mHz})$, or $\approx 2000\,\mathrm{W}\,\mathrm{m}^{-2}$ using the newer total flux estimates. These measurements then indicate that if high frequency waves are responsible for CBP emission, $50\%$ of the high frequency energy incident upon a CBP dome must be converted to thermal energy within the dome to account for the observed CBP heat flux of $1000\,\mathrm{W}\,\mathrm{m}^{-2}$. Instead, we found of order $1\%$ efficiency, so the process we have outlined above cannot be the main source of power necessary to replenish radiative losses by CBPs that have a dome–type magnetic topology. Under the most favorable conditions from our parameter study, the focusing of wave energy towards the null to generate shocks and current sheets might account for $5-10\%$ of the the radiative loss budget of CBPs. Apparently, the bulk of the energy must come from other processes, such as reconnection due emergence or continual shearing, as discussed in Galsgaard et al. (2017).

The type of wave dissipation we have studied is also unlikely to play a dominant role in larger structures. Tarr et al. (2014) determined the total radiated power from an emerging ephemeral region, again with a dome topology similar to the one we have simulated here. That region contained a coronal null point $\approx 15\,\mathrm{Mm}$ above the photosphere whose fan surface intersected the photosphere in a ring roughly $30-50\,\mathrm{Mm}$ across, covering an area $\approx 10^{15}\,\mathrm{m}^{2}$. They found a dome–shape structure in many EUV and X–Ray wavelengths that clearly coincides with the null’s separatrix surface (c.f. their Figures 2,3,5, and 7). The dimensions for this structure are rather larger than our simulated dome rotated about its central axis. The total radiated power from this region varied between $10^{18}-10^{19}\,\mathrm{W}$ during emergence, giving a radiated flux of some $2000-4000\,\mathrm{W}\,\mathrm{m}^{-2}$ from the dome, essentially the same as for CBPs. Once again, clearly another process, such as the reconnection between emerging and external fields proposed by Tarr et al. (2014), must play a more dominate role in that case.

On the other hand, the formation of the upward propagating slow shocks in our simulation may help explain what drives type–II spicules (de Pontieu et al., 2007). Martínez-Sykora et al. (2017), particularly in their supplementary material, presented a simulation of the chromosphere that includes sets of transverse waves that mode convert into slow mode waves and associated flows, which they identify as type–II spicules (see their Figure S6 and its caption). At locations where they identify the conversion process, both the field line structure and the existence of closed–contours of the $\beta=1$ surface point to the presence of magnetic X–lines or null points: see, for example, their Figure S3 at location $\approx(91.5,2\,\mathrm{Mm})$, or Figure S6 at location $(22,2\,\mathrm{Mm})$. Furthermore, out–of–plane currents localize at the approximate locations of the separatrices attached to the conversion points, similar to the concentration of currents we find in our simulations (see, e.g., their Figure S3, especially panel D). This suggests that a similar mechanism leads to both their results and ours, even though their simulations include the additional physics of ambipolar diffusion.

The slow mode shocks we find at first seem similar to those reported in McLaughlin et al. (2009). They, too found slow mode shocks traveling outward from a null after a fast mode wave passed through. However, these shocks arise from a different process. Their initial condition is zero $\beta$ so only fast modes are initially supported by the plasma. An incoming initially symmetric pulse eventually creates a ring of fast shocks. The shocks collide, forming a cusp, and the heating is much greater in the cusps. The pressure gradient set up in the cusps then gives rise to a set of outward propagating slow shocks. The inward propagating fast shocks then collapse the null into a current sheet and drive reconnection outflows (jets) into the cusp, further heating the plasma and driving the slow shocks.

Our slow shocks form instead by mode conversion of the fast wave, a process not immediately available when the initial condition is zero $\beta$. This difference is illustrated by considering the low amplitude limit. When McLaughlin et al. (2009) reduced their initial amplitude to keep the simulation in the linear regime (see their Appendix A), the slow shocks disappear from the solution while the fast mode wave simply propagates symmetrically inward with increasing energy density. In contrast, in our simulations if we reduce the amplitude of the perturbation or change the system parameters $(y_{tr})$ to make the gradient of the phase speed more gradual, the strong slow mode shocks smoothly transition into weak shocks and then slow mode waves, still localized to the separatrices, which are what we reported in Paper I.

Does this analysis carry over into 3D? Thurgood et al. (2017) shows that yes, it can for an idealized 3D null. They perform a 3D upgrade to McLaughlin et al. (2009), where the initial pulse is a donut–shaped fast wave, and they find very similar propagating density shocks to those we discussed in §3.3 (see their Figure 8). This demonstrates that the conversion of incoming waves to compressive modes remains a substantial effect in 3D. On the other hand, they also find evidence that the incoming fast wave generates outward propagating Alfvén waves, which we deliberately avoided in our 2D setup. Clearly in 3D, in order to understand how the energy of propagating waves is ultimately distributed, one must treat the conversion between all modes simultaneously. This will be the basis of a future study.

Another further step into the 3D realm has already been taken by McLaughlin et al. (2019), who performed a 3D WKB study of fast wave refraction in a 3D dome topology. They could not address mode conversion due to their cold plasma approximation, by they did determine that the ability of a null to directly capture wave energy decreases as the null height increases above the photosphere. At the same time, it is apparent that any inhomogeneity in the magnetic field affects the propagation of MHD waves. Relaxing the $\beta=0$ assumption and adding a stratified atmosphere are important next steps in determining how waves propagate in solar–like environments.

Santamaria and coauthors have discussed a similar set of simulations in a series of papers (Santamaria et al., 2015, 2016, 2017; Santamaria & Van Doorsselaere, 2018). Santamaria et al. (2017) report on a set of jet–like slow shocks leaving a null, analogous to those we analyzed in detail. They focus mostly on properties of the wave field induced by the wave–null interaction, finding that the null may act like a resonant cavity, producing higher frequencies compared to those originally introduced in their wave packet. As they further report in Santamaria & Van Doorsselaere (2018), the high frequency disturbances appear to be localized to the separatrices, so that the high frequencies they find are related to the properties of mode-converted slow mode shocks: higher amplitude, faster traveling shocks produce higher frequency dynamics.

We do not appear to find direct evidence of a resonant cavity at the null in our simulations. The apparent oscillatory reconnection we find seems to be related to the properties of the incoming wave (reconnection once in each direction corresponding the leading and trailing edges of the pulse, then repeated for each reflection off the lower boundary). Oscillatory reconnection has been considered in other studies (Craig & McClymont, 1991; Murray et al., 2009; McLaughlin et al., 2012a, b; Threlfall et al., 2012; Thurgood et al., 2019). Notably, in these simulations the region about the null is either continuously driven, say by ongoing flux emergence (Murray et al., 2009; McLaughlin et al., 2012b), or the initial perturbation has a strong symmetry about the null while the resulting dynamics include reflections off exterior boundaries (Craig & McClymont, 1991; McLaughlin et al., 2012a; Threlfall et al., 2012; Thurgood et al., 2019). Oscillatory reconnection should be a process inherent to the dynamics of the null, but appears difficult to isolate in all of the above studies. Still, the nonlinear dynamics reported in Thurgood et al. (2019) belie a simple explaination due to reflections off a boundary. That we fail to see oscillatory reconnection in our simulation could be because: (i) our resistivity of $0.0333$ is too high and therefore creates too much damping (the inherent length scale of a stratified atmosphere leads to rather high numerical resistivity when attempting to resolve both a transition-region null and a photosphere, see Tarr et al. (2017) §2.1); (ii) the efficiency of mode conversion allows the slow modes to carry away the energy that would otherwise act as the potential energy needed to cause rebound; or (iii) the wave packets that reach the null do so rather asymmetrically, while a quadrupolar (cartesian) mode about the null may be necessary to excite oscillatory reconnection (Craig & McClymont, 1991; Longcope & Tarr, 2012).

Our simulations neglected radiation, but the general effect has long been understood (Bray & Loughhead, 1974, §6.5.3). Shocks cause nonadiabatic changes: they increase the plasma entropy, while radiation removes entropy from the local plasma (indeed, plots of the change in plasma entropy $S=P\rho^{-\gamma}$ for our simulations show similar structure to Figure 17). A full shock cycle including radiation will return the plasma to its initial state. If we included radiative losses, the shocks in our simulation would lose energy faster by radiative damping of the shock front, while the tangential discontinuities in the shock wake will reduce with time as the heated plasma radiatively cools and the balance between plasma pressure gradients and Lorentz forces shifts back towards a more magnetic–potential state. Thus, for solar observations, we would expect the mode–conversion/shock process to involve a more localized intensity enhancement that quickly fades. The upward and downward propagating shocks would generate some Doppler shifts of spectral lines, but only as a component of the background emission. This behavior is similar to the observations reported in Reardon et al. (2013), who found localized, spatially spreading intensity enhancements, $10-20\%$ above the background, and perhaps $1-1.5\,\mathrm{Mm}$ at maximum extent. However, further simulations that include radiation are required for more detailed comparison.

7 Conclusion

We have studied the dynamics that arise when compressive waves are introduced in stratified, solar–like atmospheres that include magnetic nulls. A surprising result of our previous work (Paper I) was that purely compressive waves introduced in the photosphere lead to currents localized to the magnetic skeleton of the overlying field. Here, we report the mechanism by which that localization occurs. We found that the initial current sheet formation is dominated by magnetic compression localized to one side of the initial null point. Magnetic reconnection across this current sheet locally heats the plasma near the null through Ohmic dissipation. There are two distinct reconnection scenarios at work. The first is caused by fast mode refraction and subsequent conversion around null region, and is responsible for the majority of the Ohmic heating. The second appears to be due to diffusion as the fast wave reflects off the lower boundary, and thus does not substantially contribute to the heating at the null. The resulting plasma pressure gradients are balanced by the magnetic pressure term in the Lorentz force, leading to a persistent current density localized to the null.

The currents along the separatrices, on the other hand, are caused by mode–converted slow mode waves that shock as they propagate away from the null. The strongest shocks form underneath and adjacent to the separatrix dome, so the plasma is heated more in those locations compared to the opposite sides of the separatrices. As the higher pressure regions expand they compress the magnetic field on the opposite side of the separatrix, generating the long–lived currents localized to each separatrix.

The fraction of the injected wave packet’s energy which eventually heats the plasma depends on the region of the dome where the packet is introduced and the properties of the stratified atmosphere. We found that our previous (Paper I) WKB estimate of this dependence matched our full MHD simulation results quite well. A greater decrease in density with height, caused by a higher transition region, allows less wave energy to reach the null, while a more compact equipartition region surrounding the null generates more focused mode conversion and stronger slow mode shocks along the separatrices, better localizing energy dissipation to the separatrix surface.

Finally, we applied our results to the question whether the trapping and dissipation of high frequency waves in magnetic dome topologies association with coronal bright points is sufficient to balance observed radiative losses. We found that the processes we described may only account for $\lesssim 10\%$ of the required enegies under the most favorable circumstances.

This work is supported by the the Chief of Naval Research and the NASA HSR and LWS programs. The simulations were performed under a grant of computer time from the DoD HPC program. This research has made use of NASA’s Astrophysics Data System. A portion of this work was carried out while L.A.T. was a visiting scientist at the University of Hawaii Institute for Astronomy, Maui, HI, and as an NRC Research Associate hosted by the Naval Research Laboratory in Washington, DC. We thank an anonymous referee for many helpful comments.

Appendix A Grid stretching

The stretched grids are defined in cell coordinates and then scaled to the total simulation length. In the $x$ direction, $L_{\text{total}}=700L_{N}$, and the width of a cell $-N_{x}/2\leq ix\leq N_{x}/2$ is given by

[TABLE]

where we use $dx_{0}=L_{N}/32,\ c_{1}=-N_{x}/2+512,\ c_{2}=N_{x}/2-512$, and $w=160$. The above function smoothly varies the cell width between $dx_{0}$, or $\Delta_{x}=L_{n}/32$, at the center of the domain (near the magnetic dome) to $\Delta_{x}=dx_{0}(1+2f_{x})$, or $\Delta_{x}=0.526L_{N}$, at the simulation edges. The geometric factor is $f_{x}/2=\frac{L_{\text{total}}/dx_{0}-N_{x}}{N_{x}+A_{x}}$, where A_{x}=2w\ln\Bigl{(}\cosh\frac{N_{x}/2-c_{1}}{w}\operatorname{sech}\frac{-c_{2}}{w}\Bigr{)}. The factor $f_{x}$ is related to the amplitude of the stretching, while $A_{x}$ is related to the fraction of stretched versus unstretched cells.

The $y$ direction is defined similarly, where now $L_{\text{total}}=800L_{N}$, $0<iy<N_{y}$, and

[TABLE]

in which $dy_{0}=L_{N}/32,\ c=1920$, and $w=160$. The $y$ direction has only one stretching center, so the geometric factors are slightly different: $f_{y}=\frac{L_{\text{total}}/dy_{0}-N_{y}}{N_{y}+A_{y}}$ and A_{y}=w\ln\Bigl{(}\cosh\frac{N_{y}-c}{w}\operatorname{sech}\frac{-c}{w}\Bigr{)}. The $y$ grid varies smoothly between $\Delta_{y}=dy_{0}=L_{N}/32$ and $\Delta_{y}=dy_{0}(1+2f)=0.34L_{N}$.