Magnetoacoustic Waves in a Stratified Atmosphere with a Magnetic Null Point

We solve the resistive magnetohydrodynamic equations in 2.5D in Cartesian coordinates using the LARE2D code (Arber et al. 2001). The horizontal and vertical directions are x and y, and the out-of-plane direction is z. We notate spatial points by ${\boldsymbol{x}}=(x,y)$ and velocities by ${\boldsymbol{v}}=({v}_{x},{v}_{y},{v}_{z})$. Individual components may also be referenced in the standard way by a subscript i or j: i = 1, 2, 3; j = 1, 2, 3. In order, the equations describe mass, momentum, and energy conservation, and the magnetic induction:

$$\begin{eqnarray}&&\displaystyle \frac{D\rho }{{Dt}}=-\rho {\rm{\nabla }}\cdot {\boldsymbol{v}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}}=-\displaystyle \frac{1}{\rho }{\rm{\nabla }}P+\displaystyle \frac{1}{\rho }{\boldsymbol{j}}\times {\boldsymbol{B}}+{\boldsymbol{g}}+\displaystyle \frac{1}{\rho }{\rm{\nabla }}\cdot S,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D\epsilon }{{Dt}}=-\displaystyle \frac{P}{\rho }{\rm{\nabla }}\cdot {\boldsymbol{v}}+{\zeta }_{{ij}}{S}_{{ij}}+\eta {j}^{2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{B}}}{{Dt}}=({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{v}}-{\boldsymbol{B}}({\rm{\nabla }}\cdot {\boldsymbol{v}})-{\rm{\nabla }}\times (\eta {\boldsymbol{j}}).\end{eqnarray} \tag{ 5 }$$

$D/{Dt}=\partial /\partial t+{\boldsymbol{v}}\cdot {\rm{\nabla }}$ is the advective derivative. Owing to the translational invariance in z, all derivatives in the out-of-plane direction are zero. Our primitive variables are mass density, specific internal energy, plasma velocity, and magnetic field, $\rho ,\epsilon ,{\boldsymbol{v}},{\boldsymbol{B}}$, respectively. The current density is defined through ${\boldsymbol{j}}={\rm{\nabla }}\times {\boldsymbol{B}}/{\mu }_{0}$, where μ₀ is the permeability of free space. The stress tensor has components ${S}_{{ij}}=\nu ({\zeta }_{{ij}}-(1/3){\delta }_{{ij}}{\rm{\nabla }}\cdot {\boldsymbol{v}})$ with ${\zeta }_{{ij}}=(1/2)(\partial {v}_{i}/\partial {x}_{j}+\partial {v}_{j}/\partial {x}_{i})$, and we use a uniform viscosity $\nu =2.8\,\mathrm{kg}\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$. We set the resistivity to $\eta =116{\rm{\Omega }}\,{\rm{m}}.$ The viscosity and resistivity are such that the Reynolds number R = 10⁵ and the Lundquist number S = 10, which we discuss in more depth below.

Temperature, internal energy, and pressure are related through the ideal gas law

$$\begin{eqnarray}&&P=\displaystyle \frac{\rho {k}_{{\rm{B}}}T}{{\mu }_{m}}\quad \mathrm{and}\quad \epsilon =\displaystyle \frac{{k}_{{\rm{B}}}T}{{\mu }_{m}(\gamma -1)},\end{eqnarray} \tag{ 6 }$$

where k_B is Boltzmann's constant and γ = 5/3 is the ratio of specific heats. Gravitational acceleration is set to the solar value of g = 274 m s⁻², with ${\boldsymbol{g}}=-g\hat{y}$. We use the fully ionized limit of MHD, but set the reduced mass to the average ion mass, ${\mu }_{m}={m}_{i}=1.25{m}_{p}$. This choice more correctly reproduces coronal densities (see discussion in Leake & Linton 2013). The factor 1.25 models the effect of heavy ions.

We normalize the governing Equations (2)–(5) by writing each variable as a constant multiplying a normalized coordinate: $x={L}_{N}{x}^{* },\rho ={\rho }_{N}{\rho }^{* },$ and so on. We set L_N = 1.50 × 10⁵ m, ρ_N = 3.03 × 10⁻⁴ kg m⁻³, and B_N = 0.12 T. ρ_N is the photospheric density, and L_N is approximately the scale height at the photosphere. All other variables may be defined through these three. The velocity normalization, for instance, is ${V}_{N}={B}_{N}/\sqrt{{\mu }_{0}{\rho }_{N}}=6177\,{\rm{m}}\,{{\rm{s}}}^{-1},$ the photospheric Alfvén speed for unit magnetic field. Time normalization is t_N = L_N/V_N = 24.28 s. The viscosity is normalized to ν_N =ρ_NV_NL_N. We define the Reynolds number R = ρ_NV_NL_N/ν using the photospheric Alfvén speed and scale height. Using this, we set the normalized viscosity to ${\nu }^{* }=\nu /{\nu }_{N}=\tfrac{1}{R}={10}^{-5},$ resulting in the value ν quoted above. From here on we take all variables to be normalized unless explicitly stated, and we suppress the star notation.

Our domain extends from the photosphere at the lower boundary up to the low corona. The initial thermodynamic equilibrium is set to be invariant in x, so that T(x, y, t = 0) = T₀(y) (a subscript 0 will be used to refer to the initial state for all variables). We model a low-temperature photosphere, steep transition region, and isothermal corona as a hyperbolic tangent function:

$$\begin{eqnarray}&&{T}_{0}(y\gt 0)={T}_{\mathrm{ph}}+\displaystyle \frac{{T}_{{cor}}-{T}_{\mathrm{ph}}}{2}\left({\mathsf{\tanh }}\left(\displaystyle \frac{y-{y}_{\mathrm{tr}}}{{w}_{\mathrm{tr}}}\right)+1\right).\end{eqnarray} \tag{ 7 }$$

The parameters are T₀(0) = T_ph = 5778 K, T_cor = 150T_ph, y_tr = 3.0 Mm, and ${w}_{\mathrm{tr}}=.75\,\mathrm{Mm}$. The initial density profile is determined by numerical integration of the hydrostatic equation, ${\rm{\nabla }}{P}_{0}=-{\rho }_{0}g\hat{{\boldsymbol{y}}}$, with ρ₀(0) = 1.0ρ_N as a boundary condition. Figure 1(a) shows the initial hydrostatic equilibrium for temperature (left axis) and density (right axis) on a log scale.

We use a uniform staggered grid of 1024 × 1024 cells, with ρ and defined at cell centers, ${\boldsymbol{v}}$ at cell vertices, and ${\boldsymbol{B}}$ at cell faces (see Arber et al. 2001, for details). The cell width is Δ_x = Δ_y = Δ = L_N/8, which sets the size of our domain to x = (−9.6, 9.6) Mm and y = (0, 19.2) Mm.

Finally, we estimate the magnitude of numerical diffusion caused by the finite-difference scheme using the method of Arber et al. (2007). We consider the 1D analog to Equation (5), ${\partial }_{t}B+C{\partial }_{x}B=0$, Taylor-expand the second-order finite-difference equation, and isolate the error term, whose coefficient is the effective numerical resistivity, η_num. We find that ${\eta }_{\mathrm{num}}=C{{\rm{\Delta }}}^{2}/6L$, where L is a typical length of the dynamic evolution and C the fastest wave speed. For a worst-case scenario corresponding to a shock across three cells, L = 3Δ. If $C\approx 20{V}_{N}$ at locations where shocks form, we find that ${\eta }_{\mathrm{num}}\approx {\rm{\Delta }}/{L}_{N}=0.15$ in the normalized units defined below. We have run multiple simulations varying only the explicit resistivity and checking the solution in regions of strong gradients (there is little discernible effect outside these regions). We found that explicit resistivity begins to dominate over the numerical resistivity at the expected value η ≈ 0.1, confirming our approximate calculations. We therefore use this value of resistivity for our simulations. This ensures that resistive effects are due to the explicit term in Equations (4) and (5), while simultaneously keeping the resistive effects as small as possible, for the chosen grid resolution. 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 results are the values of resistivity and Lundquist number stated above.

To the hydrostatic background, we add an initial magnetic field derived from a flux function, ${\boldsymbol{B}}={\rm{\nabla }}\times {\boldsymbol{A}}=({\rm{\nabla }}{A}_{z})\times \hat{z}$,

$$\begin{eqnarray}&&{A}_{z}(x,y)=\displaystyle \sum _{p=-1}^{1}\displaystyle \frac{{\psi }_{p}}{\pi }\,\mathrm{atan}\,\displaystyle \frac{y-{y}_{s}}{x-{x}_{p}},\end{eqnarray} \tag{ 8 }$$

while A_x and A_y are both zero. As Figure 1(b) shows, the flux function includes three sources of flux (per ignorable length) ψ_p, each located at depth y_s = −0.825 Mm, and having horizontal locations ${x}_{0}=0\,\mathrm{Mm},{x}_{\pm 1}=\pm 3.75\,\mathrm{Mm}.$ The inner source ψ₀ has the opposite sign of flux of the outer two, ${\psi }_{\pm 1}$. The resulting magnetic field contains one null point, located at

$$\begin{eqnarray}&&{y}_{\times }\,=\,-{y}_{s}+\sqrt{\displaystyle \frac{{x}_{1}^{2}}{| 2{\psi }_{1}/{\psi }_{0}| -1}}.\end{eqnarray} \tag{ 9 }$$

We choose y_× = 3.75 Mm, which sets $| {\psi }_{1}/{\psi }_{0}| \approx 0.836,$ and ${\psi }_{0}=-5{B}_{N}{L}_{N}=-9\times {10}^{10}\,\mathrm{Mx}/{L}_{\mathrm{ignorable}}$, creating a maximum vertical field at the lower boundary of ≈1 kG for each source. The field strength along the lower boundary is shown in the bottom panel of Figure 1(b).

The magnetic field generated by Equation (8) is a potential field, whose z component analytically satisfies Laplace's equation ${\mu }_{0}{j}_{z}{(x,y)=({\rm{\nabla }}\times {\boldsymbol{B}})}_{z}=-{{\rm{\nabla }}}^{2}{A}_{z}=0$, and j_x = j_y = 0 as well. However, the second-order finite-difference scheme used by LARE2D results in spurious currents, as may be found by substituting the Taylor series expansion of Equation (8) into the finite-difference scheme. All odd-order derivative terms cancel numerically, the second derivative term correctly reproduces the Laplacian, but even-order derivative terms of order 4 and higher do not cancel in general. The initial condition is therefore slightly out of force balance. We allow the system to come to equilibrium during an initial relaxation period for each simulation. The resulting changes are minimal, though the largest changes are near the X-point. The change in B²(x, y) is everywhere <1.1% and less than 0.1% for more than 99.5% of the domain; the change in internal energy density is everywhere less than 0.1%. For the remainder of this paper, we refer to t = 0 as the end of this initial relaxation and only discuss dynamics for t > 0. Subscripted variables ${\rho }_{0},{{\boldsymbol{B}}}_{0},\ldots$ therefore refer to this background initial condition.

LARE2D requires two ghost cells surrounding the domain to implement the boundary conditions. In the convection zone, the increasing sound speed will eventually cause waves to refract upward at a depth that depends on the horizontal wavelengths of the waves. We therefore use a reflective lower boundary and add to it a time-dependent driver, as described in Section 4. The phase of the waves reflected from our lower boundary will not be accurate, but this does not affect our analysis in Section 5, which focuses on the region near the null before any reflections reach it.

The side boundaries are line tied, and the top boundary is zero gradient. Linear damping regions are implemented for both, and these reduce the velocities in a cell based on distance to the boundary. This has the effect of removing kinetic energy from the system. The amount of reduction starts at 0 at $y=102.4{L}_{N}\ (15.36\,\mathrm{Mm})$ or $| x| =51.2{L}_{N}\ (7.68\,\mathrm{Mm})$ and increases toward the appropriate boundary (left, right, or top). If L_damp is the size of the damping region $(25.6{L}_{N}$ in y, $12.8{\L}_{N}$ in x) and ΔX is the distance from a cell within the damping region to the start of the region, then at each time step the velocities in that cell are reduced by a factor $1+0.3{\delta }_{t}\tfrac{{\rm{\Delta }}X}{{L}_{\mathrm{damp}}}$, where δ_t is the numerical time step. Our analysis will focus on dynamics near the null point and magnetic dome, and testing has shown that the damping regions prevent the majority of reflections, leading to a small contribution to the dynamics at locations of interest.

An important aspect for wave propagation is the structure of the sound and Alfvén speed throughout the domain. Hydrostatic stratification creates a sound speed that varies vertically only: ${c}_{s}(y)=\sqrt{\gamma {k}_{b}{T}_{0}(y)/{\mu }_{m}}=\sqrt{\gamma {P}_{0}(y)/{\rho }_{0}(y)}$. This vertical sound speed stratification is shown as gray scale in Figure 2(a) and as the solid black line in Figure 2(b). The sound speed ranges from ≈8 to 100 km s⁻¹ within the domain.

Zoom In Zoom Out Reset image size

The Alfvén speed is structured horizontally as well as vertically: ${v}_{{\rm{A}}}(x,y)=\tfrac{{B}_{0}(x,y)}{\sqrt{{\mu }_{0}{\rho }_{0}(y)}}$. Thin solid blue lines in Figure 2(a) show contours of the Alfvén speed, in steps of 35 km s⁻¹, with the value indicated at several levels. The function has a minimum at the null point where v_A(0, 3.75 Mm) = 0, and a maximum of v_A ≈ 400 km s⁻¹ in the two lobes on either side at ${\boldsymbol{x}}\approx (\pm 5,+3.75)\,\mathrm{Mm}$. There are two saddle points along x = 0, at y ≈ 2.4 and 9.1 Mm, though these are not on the same contour level. Figure 2(b) shows the Alfvén speed along two vertical slices, one passing through the null at x = 0 Mm (red), the other passing near the left maximum at x = −5 Mm (dashed line).

We will use LARE2D to solve the full, nonlinear set of Equations (2)–(5). However, we will also use linear theory to determine the form of a wavepacket injected into the system, approximate its propagation through the (stationary) background, and help interpret the simulation output. For this linear analysis, we set the viscous and resistive terms to zero and add a perturbation to each background quantity: $\rho ={\rho }_{0}+{\rho }_{1},P={P}_{0}+{P}_{1},{\boldsymbol{B}}={{\boldsymbol{B}}}_{0}+{\boldsymbol{b}}$, and ${\boldsymbol{v}}={{\boldsymbol{v}}}_{1}({\rm{with}}\,{\rm{}}{{\boldsymbol{v}}}_{0}=0)$. The linearized ideal MHD equations are then

$$\begin{eqnarray}&&{\partial }_{t}{\rho }_{1}+{{\boldsymbol{v}}}_{1}\cdot {\rm{\nabla }}{\rho }_{0}+{\rho }_{0}{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{1}=0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rho }_{0}{\partial }_{t}{{\boldsymbol{v}}}_{1}={\rho }_{1}{\boldsymbol{g}}-{\rm{\nabla }}{P}_{1}+\displaystyle \frac{1}{{\mu }_{0}}({\rm{\nabla }}\times {\boldsymbol{b}})\times {{\boldsymbol{B}}}_{0},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\partial }_{t}{P}_{1}={c}_{s}^{2}{\partial }_{t}{\rho }_{1}+{c}_{s}^{2}{{\boldsymbol{v}}}_{1}\cdot {\rm{\nabla }}{\rho }_{0}-{{\boldsymbol{v}}}_{1}\cdot {\rm{\nabla }}{P}_{0},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{b}}={\rm{\nabla }}\times ({{\boldsymbol{v}}}_{1}\times {{\boldsymbol{B}}}_{0}).\end{eqnarray} \tag{ 13 }$$

We have assumed an adiabatic relation between the perturbed pressure and density, $\delta {P}_{1}={c}_{s}^{2}\delta {\rho }_{1}$. These equations may be combined to yield a conservation relation for wave energy density and energy flux (Bray & Loughhead 1974; see Appendix B):

$$\begin{eqnarray}&&{\partial }_{t}\left[\displaystyle \frac{1}{2}{\rho }_{0}{v}_{1}^{2}+\displaystyle \frac{{P}_{1}^{2}}{2{\rho }_{0}{c}_{s}^{2}}+\displaystyle \frac{| {{\boldsymbol{b}}}^{2}| }{2{\mu }_{0}}-\left(\displaystyle \frac{1}{2}{\partial }_{y}{\rho }_{0}g+\displaystyle \frac{{\rho }_{0}{g}^{2}}{2{c}_{s}^{2}}\right){Y}^{2}\right]\\ &&\quad +{\rm{\nabla }}\cdot \left[{P}_{1}{{\boldsymbol{v}}}_{1}+\displaystyle \frac{1}{{\mu }_{0}}({{\boldsymbol{B}}}_{0}\times {{\boldsymbol{v}}}_{1})\times {\boldsymbol{b}}\right]=0.\end{eqnarray} \tag{ 14 }$$

We can identify each term in Equation (14) in order: the energy densities are kinetic (E_K), acoustic (E_A), magnetic (E_M), and gravitational (E_G), while the energy fluxes are acoustic $({{\boldsymbol{F}}}_{{\rm{A}}})$ and magnetic (Poynting, ${{\boldsymbol{F}}}_{M}$):

$$\begin{eqnarray}{\partial }_{t}E & = & {\partial }_{t}({E}_{K}+{E}_{A}+{E}_{M}+{E}_{G})=-{\rm{\nabla }}\cdot {\boldsymbol{F}}\\ & = & -\,{\rm{\nabla }}\cdot ({{\boldsymbol{F}}}_{A}+{{\boldsymbol{F}}}_{M}).\end{eqnarray} \tag{ 15 }$$

We can further group the energy density terms into kinetic (E_K) and potential (E_A, E_M, E_G) parts. In the gravitational term, $Y(t)={\int }_{0}^{t}{v}_{y}({t}^{\prime }){{dt}}^{\prime }$ is the vertical displacement of a fluid element. For frequencies greater than a few times the Brünt–Väisälä frequency (∼4.5 mHz), the gravitational term becomes increasingly less important than the other terms (Hansen & Cally 2009). We will focus on frequencies of ∼40 mHz, and therefore we drop this term from our analysis.

The above equation involves real quantities that may be directly calculated from the output of the simulation at each time t: ${P}_{1}(t)=P(t)-{P}_{0},{\boldsymbol{b}}={\boldsymbol{B}}(t)-{{\boldsymbol{B}}}_{0}$, and so on. This will be useful for determining the amount and type of wave energy at various locations in the simulation. On the other hand, how the wave energy propagates through the system should be controlled by the properties of the background state. To determine wave propagation, we first combine the linearized Equations (10)–(13) into a single wave equation for the velocity perturbations by taking the time derivative of the momentum equation and substituting in the continuity, energy, and induction equations. The result is

$$\begin{eqnarray}{\partial }_{{tt}}^{2}{{\boldsymbol{v}}}_{1} & = & {c}_{s}^{2}{\rm{\nabla }}({\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{1})+\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}\{{\rm{\nabla }}\times [{\rm{\nabla }}\times ({{\boldsymbol{v}}}_{1}\times {{\boldsymbol{B}}}_{0})]\}\\ & & \times \,{{\boldsymbol{B}}}_{0}+{\boldsymbol{g}}(\gamma -1){\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{1}+{\rm{\nabla }}({\boldsymbol{g}}\cdot {{\boldsymbol{v}}}_{1})\\ & & +({\boldsymbol{g}}\cdot {{\boldsymbol{v}}}_{1})\left[\displaystyle \frac{{\rm{\nabla }}{c}_{s}^{2}}{{c}_{s}^{2}}+\displaystyle \frac{{c}_{s}^{2}}{{\rho }_{0}}{\rm{\nabla }}\left(\displaystyle \frac{{\rho }_{0}}{{c}_{s}^{2}}\right)-\displaystyle \frac{{\boldsymbol{g}}}{{c}_{s}^{2}}\right].\end{eqnarray} \tag{ 16 }$$

We assume that each perturbed quantity varies in space and time only by a common phase term: ${{\boldsymbol{v}}}_{1}({\boldsymbol{x}},t)={\boldsymbol{a}}{e}^{i{\rm{\Phi }}},{\rm{\Phi }}={\boldsymbol{k}}\cdot {\boldsymbol{x}}-\omega t$. Next, we apply the WKB approximation (Weinberg 1962) that the phase function varies much more rapidly than any background quantity: $k=| {\rm{\nabla }}{\rm{\Phi }}| \gg 1/h$, where h represents any spatial scale of the background. We again drop the explicit gravitational terms from consideration, though note that part of the stratification's effect is implicit through the spatial dependence of c_s and ρ₀. After applying these assumptions, we can write the wave equation in dyadic notation (see Thomas 1982; Campos 1983, who kept the explicit gravitational term):

$$\begin{eqnarray}{\omega }^{2}{\boldsymbol{I}}\cdot {{\boldsymbol{v}}}_{1} & = & [{c}_{s}^{2}{\boldsymbol{kk}}-{v}_{{\rm{A}}}^{2}(\hat{{\boldsymbol{b}}}\cdot {\boldsymbol{k}})\hat{{\boldsymbol{b}}}{\boldsymbol{k}}+{v}_{{\rm{A}}}^{2}{(\hat{{\boldsymbol{b}}}\cdot {\boldsymbol{k}})}^{2}{\boldsymbol{I}}\\ & & +{v}_{{\rm{A}}}^{2}{\boldsymbol{kk}}-{v}_{{\rm{A}}}^{2}(\hat{{\boldsymbol{b}}}\cdot {\boldsymbol{k}}){\boldsymbol{k}}\hat{{\boldsymbol{b}}}]\cdot {{\boldsymbol{v}}}_{1},\end{eqnarray} \tag{ 17 }$$

where ${\boldsymbol{I}}$ is the unit dyad. Setting the determinant of this equation to zero, we find the dispersion relations for Alfvén (ω_A) and fast and slow magnetoacoustic waves (ω_±):

$$\begin{eqnarray}{ \mathcal D }(\omega ,{\boldsymbol{k}},{\boldsymbol{x}}) & = & ({\omega }_{{\rm{A}}}^{2}-{k}^{2}{v}_{{\rm{A}}}^{2}{\cos }^{2}\theta )\left({\omega }_{\pm }^{2}-{k}^{2}\left[\displaystyle \frac{1}{2}({v}_{{\rm{A}}}^{2}+{c}_{s}^{2})\right.\right.\\ & & \pm \left.\left.\displaystyle \frac{1}{2}\sqrt{{v}_{{\rm{A}}}^{4}+{c}_{s}^{4}-2{v}_{{\rm{A}}}^{2}{c}_{s}^{2}\cos 2\theta }\right]\right)=0.\end{eqnarray} \tag{ 18 }$$

The angle between the propagation direction and the magnetic field (the attack angle or phase angle) is defined through $\hat{{\boldsymbol{k}}}\cdot \hat{{\boldsymbol{b}}}=\cos (\theta )$. Equation (18) has the same form as the standard relation for a homogeneous compressible plasma (Kulsrud 2005). The WKB approximation we have used simply takes ${v}_{{\rm{A}}}({\boldsymbol{x}})$ and c_s(y) to be spatially varying functions rather than uniform. The determinant factors into ${ \mathcal D }={{ \mathcal D }}_{{\rm{A}}}{{ \mathcal D }}_{+}{{ \mathcal D }}_{-}$, which describe separate conditions on ${\boldsymbol{k}}$ for a given frequency ω for the Alfvén, fast, and slow mode, respectively (we will generically apply the subscripts +, −, A to indicate solutions for each branch). Note that in each case the relation between frequency and wavenumber is linear, of the form ω ∝ k. This is the relation for dispersionless waves, where each frequency wave propagates in the same way, just as for MHD waves in a homogeneous plasma. Retaining the gravitational terms in Equation (16) introduces dispersion at low frequencies near the Brünt–Väisälä frequency. Our first-order WKB approximation therefore includes the effect of refraction but excludes dispersion, which is unimportant for high-frequency waves.

Each mode will propagate through the system along a different path. As explained in more detail in Appendix A, we can trace the path of a ray (say, the fast ray) ${\boldsymbol{x}}(\tau )$, where τ parameterizes the distance along the ray, by picking an initial condition ${\boldsymbol{k}}(\tau =0),{\boldsymbol{x}}(\tau =0)$ and solving for ${\boldsymbol{k}}(\tau )$ and ${\boldsymbol{x}}(\tau )$ subject to the constraint that the correct condition is satisfied, i.e., D₊ = 0 for the fast ray). The result is that the ray satisfies Hamilton's equations

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{k}}}{{dt}}={-\displaystyle \frac{\partial \omega }{\partial {\boldsymbol{x}}}\right|}_{{\boldsymbol{k}}}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\displaystyle \frac{\partial \omega }{\partial {\boldsymbol{k}}}\right|}_{{\boldsymbol{x}}},\end{eqnarray} \tag{ 20 }$$

where ω comes from the condition ${ \mathcal D }=0$. We give explicit expressions for these equations for the fast and slow ray in Appendix A; the Alfvén waves simply follow field lines. Note that ${\partial }_{{\boldsymbol{k}}}\omega ={{\boldsymbol{v}}}_{g}$, the wave's group velocity. This is the velocity at which energy is transferred along the ray, and is equal to ${\boldsymbol{F}}/E$ from Equation (15). The energy relation is true up to the assumption that there is no dissipation and that the wave stays on a single branch of this dispersion relation, a point that we will return to repeatedly in the following.

The fast and slow branches of the dispersion relation change their characters as the plasma shifts between pressure and magnetically dominated regions. The fast mode is increasingly acoustic (potential energy dominated by E_A) for β > 1 and magnetic (dominated by E_M) for β < 1, while the opposite is true for the slow mode. Mode coupling is allowed under certain conditions (Tracy et al. 2003), where the essential requirement is that the gradient of the phase function Φ for each mode be similar in a region of space, so that, for example, ${{\boldsymbol{k}}}_{+}(\omega )\approx {{\boldsymbol{k}}}_{-}(\omega )$. This can occur where the phase speeds (${v}_{\phi }=\omega /k$) of the two modes are approximately equal, along the equipartition curves where c_s = v_A.

Cally (2007) and Hansen & Cally (2009) explain the typical terminology used in the helioseismic literature for mode conversion and transmission. Physically this depends on which term dominates Equation (15), E_K, E_A, or E_M, and if that changes as a wave propagates. Conversion refers to waves whose energy shifts from the magnetic to acoustic term (or vice versa) while staying on the same branch of the dispersion relation; transmission refers to a wave whose dominant potential energy term remains nearly the same as it propagates. Thus, an acoustic fast wave originating from below the equipartition height may continue propagating as a magnetic fast wave above the equipartition height, and mode conversion is said to have taken place. Because the ray theory is a solution along a single branch of the dispersion relation, it assumes perfect conversion. These definitions depend only on the properties of the wave, as can be determined directly from the full MHD simulation presented in Section 5, and are therefore to be preferred.

Figure 3 shows the trajectories for bundles of fast magnetoacoustic rays, which are solutions for the ${{ \mathcal D }}_{+}$ branch of the dispersion relation. The slow waves (solutions to D₋) simply follow the field lines to a high degree of accuracy, and we do not show them here. The thick red and magenta lines again show the separatrices and the two contours where c_s = v_A, while select contours of the Alfvén speed from Figure 2 are shown as thin blue lines, for easy comparison. Over this background, panel (a) shows the paths for a bundle of fast rays initialized from the lower boundary as vertically propagating ($\hat{{\boldsymbol{k}}}=\hat{{\boldsymbol{y}}}$) fronts between x ≈ [−4.5, −2.0], with each ray depicted in a different color. Clearly visible is the importance of refraction due to the inhomogeneous background: rays refract away from (toward) regions of high (low) phase speed in an amount that depends on the angle between their propagation direction and the local magnetic field direction. From left to right in the figure, the rays refract back down toward the photosphere (purple to yellow rays), escape through the side or upper boundaries without passing near the c_s = v_A region surrounding the null (yellow to orange), and pass near the conversion region surrounding the null point (orange to red).

Zoom In Zoom Out Reset image size

An important feature of the ray solution is that some rays propagate directly through the c_s = v_A region surrounding the null point. This is because $\beta \ne 0$ and in the WKB approximation the fast ray describes a wave that transitions smoothly from magnetically dominated to acoustically dominated perturbations. These rays cross each other multiple times, forming multiple sets of caustics. Afanasyev & Uralov (2012) have studied similar caustics both analytically and visually by determining ray trajectories near a linear 2D null point for a $\beta \ne 0$ plasma with uniform density and temperature. In our case, tracing the paths of many rays that pass near the null point (not shown here) results in a very similar pattern of caustics to those in Figure 3 of Afanasyev & Uralov (2012). The differences arise because of the stratification: near the null, our magnetic field is fairly linear in the horizontal direction, but not in the vertical direction.

In contrast, for a β = 0 plasma the slow-mode solution vanishes from the dispersion relation, and the wave speed of the fast rays decreases to zero at the null. In that case, the ray paths form logarithmic spirals focused on the null point (McLaughlin & Hood 2006a; Longcope & Tarr 2012; in particular, see Figure 7 of the former), accumulating strong currents at the null, which can then dissipate. As has been noted before, the existence of a finite plasma pressure term thus makes focusing of wave energy on a null point more difficult than it otherwise would be. We will explore this in more detail below by direct numerical simulation.

Figure 3(b) depicts a connectivity graph via a set of bounding rays for initially vertically propagating fast rays initialized in the left half of the domain, between x ≈ −3 and −0.5 Mm. The two leftmost green rays are traced from the centers of adjacent computational cells and illustrate how quickly ray paths may diverge. Rays traced farther to the left of these (initial location x ≲ −3 Mm) all have turning points and close back down at the photosphere. This is the typical behavior for low-β fast modes when the Alfvén speed increases with height. Rays traced to the right of x ≈ −3 Mm display a much different behavior and increasingly refract toward the null point. The two yellow lines are bounding curves for rays that refract strongly toward the null point. Rays launched between the rightmost green ray and the center of the domain (between x = −0.5 and 0 Mm) again exhibit a turning point and refract back downward to the photosphere. The pattern repeats in mirror image for rays launched from the positive x side of the domain, due to the reflectional symmetry of the system.

The path of a ray depends on its initial propagation direction, $\hat{{\boldsymbol{k}}}$, in addition to its initial location. Figure 4 shows the paths of fast rays launched from a single initial position near the lower boundary, ${{\boldsymbol{x}}}_{0}=[-1.0,0.08]\,\mathrm{Mm}$, and initial propagation directions ${\chi }_{0}=\mathrm{atan}({k}_{x}/{k}_{y})$ up to ±15° to the vertical (χ = 0°), in steps of 1°. These values correspond to the center of the wavepacket we use to drive the numerical simulation, described in Section 4, and its initial range of propagation angles, as determined in Section 6. The plus signs along each ray are equally spaced in units of travel time, at intervals of $0.5{t}_{N}$, or, equivalently, through the phase distance τ. Each set of equal time points thus traces out a front of constant phase ${\rm{\Phi }}({\boldsymbol{k}}(\tau )\cdot {\boldsymbol{x}}(\tau )-\omega t(\tau ))$. The phase fronts are more closely spaced in the lower portion of the figure, where both the sound and Alfvén speeds are small, and more widely spaced higher up, where the phase speeds are greater. The set of all phase fronts fills in the phase function throughout the domain. Rays passing close to the null equipartition region again display complex trajectories. A steady prograde change in initial angle of propagation does not lead to a steady prograde change in the direction of the outgoing ray. Instead, outgoing rays jump from prograde to retrograde change several times. The rays cross, and so the phase function is, in general, multivalued throughout the domain.

Zoom In Zoom Out Reset image size

We can combine the ray-tracing results to try to predict how an initial, localized disturbance will propagate through the system. A spatially localized wavepacket launched from the lower boundary will exhibit a combination of the behaviors illustrated by the ray paths of Figures 3 and 4. At each spatial location, the packet will be best described by a range of propagation directions. In Section 6 we will need to know the range of initial propagation angles that pass through the null point's equipartition layer. To determine this, we initialized fast rays between x = −5 Mm and x = 0.15 Mm with propagation directions χ =  ±15° in steps of 0.1° We recorded each ray's closest approach, d, to the centroid of the equipartition layer, located at ${{\boldsymbol{x}}}_{c}=[0,3.9]\,\mathrm{Mm}$ (note that the equipartition centroid is shifted slightly above the null point, due to the stratification).

The null's equipartition curve is approximately circular, with a radius r_e ≈ 0.75 Mm. Figure 5(a) shows the result for angles between ±10°, with the closest approach distance for each ray displayed in logarithmic scale as a function of the initial ray location (abscissa) and initial propagation angle (ordinate). The distribution is nearly flat outside of this angle range, as seen in Figure 5(b). Rays with positive (negative) angles initially propagate in the positive (negative) x-direction. The red, green, and blue contours correspond to 1, 1.3, and 3 times r_e. We take ${R}_{n}=1.3{r}_{e}=0.96\,\mathrm{Mm}$ to define a region of influence for the null. This region is based on the equipartition scale height, ${H}_{e}^{-1}={{\partial }_{\tau }({c}_{s}^{2}/{v}_{{\rm{A}}}^{2})| }_{{c}_{s}={v}_{{\rm{A}}}}$, so that ${R}_{n}={r}_{e}+{H}_{e}$. Therefore, rays initialized at the lower boundary with initial propagation angles falling between the green curves are strongly influenced by the null point.

Zoom In Zoom Out Reset image size

The three vertical lines in Figure 5(a) show the location of slices through the contour plot and are displayed as line plots in Figure 5(b). The range of angles for which rays pass near the null reaches a maximum for rays initialized near x = −1 Mm and a minimum around x = −3.8 and 0 Mm, where the separatrices intersect the lower boundary (squares in panel (a)). The actual range of angles undergoes significant variation, is typically not centered on 0°, and has multiple inflection points. Only two sets of vertically propagating rays, one between x = −3.0 and −0.5 Mm and the other in a narrow range around x₀ = 0, pass close to the null. The first of these corresponds to the region between the yellow rays of Figure 3(b). Several dark bands run through the entire ${x}_{0}-{\chi }_{0}$ plot and indicate rays that pass very close to the center of the interaction region, almost directly through the null. There are typically multiple distinct dark bands for a given x₀. This is a manifestation of the prograde-retrograde-prograde behavior of outgoing rays discussed for Figure 4(b). The angles at which the behavior changes for a given x₀ are located at the minima along vertical slices through the plot (see, e.g., the slices plotted in Figure 5(b)). Interestingly, the dark band that arises around x₀ = −1.5 Mm bifurcates at x₀ = −0.55 Mm. Moving farther to the right, the three distinct bands at that point remain distinct across x₀ = 0, which can be seen as the three small dips in the solid curve of Figure 5(b).

We have found the ray-tracing analysis to be a useful tool for understanding the properties of our model atmosphere and, as will become apparent, for analyzing and interpreting the MHD simulation presented in the next two sections. Here, we summarize our results from the ray theory and point out several important differences from previous ray-tracing investigations.

We used a WKB method to estimate how a fast-mode wavepacket crossing the lower boundary will propagate through the background field. The WKB theory assumes that a ray stays on a single branch of the dispersion relation. A wavepacket is a spatially localized disturbance with a distribution of wavevectors ${\boldsymbol{k}}$. Thus, to see how a wavepacket will propagate, we traced rays from multiple initial positions and using multiple initial phase angles, with several examples shown in Figures 3 and 4. We found that some rays will refract toward the null point, while most will refract back downward toward the photosphere. A small segment between these two regions appears able to propagate to the top of our domain (and hence to infinity), as shown by Figure 3(b) for initially vertically propagating rays. For rays that refract back downward, the height of the turning point depends on the gradient of the phase velocity along each ray path (we have not attempted to analytically derive the turning point for our horizontally inhomogeneous background). Rays reaching the interaction region surrounding the null exhibit complex behavior, with some eventually refracting back toward the lower boundary, some reaching the top and side boundaries, and many adjacent rays crossing paths, some at multiple locations. The set of fast rays that pass near the null is a function of both the initial ray location and the propagation direction.

We did not discuss solutions for the slow waves using ${{ \mathcal D }}_{-}$ in any detail. The properties of slow MHD waves are well known for homogeneous plasmas and can be found, for instance, in Section 5.3.2 of Goedbloed & Poedts (2004). In that case, the maximum departure of the group velocity from the magnetic field direction, the return angle θ_R, is θ_R = −30°. It occurs for c_s/v_A = 1 and for a phase angle θ = 0°. For different phase angles and different values of the ratio c_s/v_A, the return angle is much smaller, typically ≲5°. To check that this carries over to the inhomogeneous case, we initialized slow rays at several hundred initial locations and wavevectors. For a given ray, the maximum departure from the magnetic field direction is typically less than 5°, and less than 30° for all rays, as in the homogeneous case. For most of the distance along each ray we find departures of <.5° Rays that do make a substantial angle to the magnetic field typically do so only for small portions of their trajectories near the c_s = v_A boundaries, after which they closely follow the field.

Our analysis differs from ray tracing used in studies of sunspots with slowly varying magnetic fields (Cally 2007; Khomenko & Collados 2009; Khomenko et al. 2009). They typically focus on the low-frequency dispersion and the variation of the cutoff frequency due to the magnetic field, which we ignore. These are very important effects for the interpretation of phase relations determined from observations (see Felipe 2012, and references therein).

Other authors have performed ray tracing through magnetic configurations containing nulls (McLaughlin & Hood 2004, 2006a; McLaughlin et al. 2008; Afanasyev & Uralov 2012; Longcope & Tarr 2012; McLaughlin et al. 2016). Here, the focus is typically quite different compared to the sunspot studies and addresses whether wave energy is reflected by the null or accumulated at the null, how to determine characteristic damping timescales, etc. We will consider these studies more in Section 7.

It is important to keep in mind how altering our initial condition would modify the ray behavior (the ease with which this is accomplished is one of the main advantages of the WKB method). Suppose that the magnetic field is everywhere reduced to zero, but the stratification is kept the same. Then upward-propagating fast rays, now degenerate to purely hydrodynamic sound waves, still refract off the increasing sound speed and eventually reflect from a height where their frequency divided by wavenumber matches the local sound speed. If instead we keep the magnetic field the same but set the density and temperature to uniform values, we find that a set of rays spiral in toward the null, as has been found many times before for both linear and more complicated nulls (McLaughlin & Hood 2004, 2006a; Longcope & Tarr 2012). Comparison with Figures 3 and 4 demonstrates that the vertical stratification breaks the radial symmetry close to the null and prevents the logarithmic inspiral. These differences arise just from the linear (WKB) theory and are in addition to any shock formation or mode conversion that may arise when solving the nonlinear MHD equations, as we will do in Section 5. The combination of stratified atmosphere, compressive waves, and nontrivial topology has not been well studied and represents one significant advance of this present work.

Not included in our WKB analysis is the rate of mode conversion between fast and slow modes. So, although many fast-ray paths do pass through the interaction region around the null, as shown in Figure 5, it remains to determine how much energy stays on the ${{ \mathcal D }}_{+}$ branch and how much converts to the ${{ \mathcal D }}_{-}$ branch. Slow waves thus generated will ultimately have different trajectories upon exiting the conversion region compared to the fast waves. We now turn our attention to the numerical solution of Equations (2)–(5) to answer this question.

Our analysis of the simulation will focus on the energy density and flux terms of Equation (14). Figure 6 contains 12 frames from a movie showing the interaction of the introduced perturbation with the null point. The images consist of three color channels—red, green, and blue—whose intensities correspond to the acoustic, kinetic, and magnetic energy density terms from Equation (14). The color channels are additive, e.g., green+red = yellow, as shown by the color wheels in the upper left panel. Each color wheel is scaled linearly between zero (black) at the center and 1 × 10⁻² J m⁻² at the radii where the disks overlap. The outer ∼30% of each disk is of uniform intensity to most effectively show the overlap of each energy term with the others. Overlap of all three terms shows up as white in the figure.

The energy density terms correspond to the energy carried by the waves, and together with the fluxes $({{\boldsymbol{F}}}_{A},{{\boldsymbol{F}}}_{M})$, they satisfy a conservation relation for wave energy, independent of the system's total energy. Because the wavepacket greatly expands and breaks into many individual pieces, no level of saturation for depicting the energy densities is entirely satisfactory. The saturation level of 1 × 10⁻² J m⁻² in Figure 6 was chosen to fairly accurately depict the energy densities after the wavepacket has rapidly expanded. Early in the simulation, say, at t = 2.1t_N, the peak energy density of the wavepacket is around 42 J m⁻², and the figure is highly saturated. The energy density at the peak of the wavepacket (${\boldsymbol{x}}=(-1.,0.18)\,\mathrm{Mm}$) is proportioned 32%, 50%, and 18% between acoustic, kinetic, and magnetic terms, respectively, and the kinetic and potential terms are in equipartition. Later on, the pulse expands into the low-β regions and breaks into multiple packets. The peak energy densities of the yellow pulses at 6.3t_N, near ${\boldsymbol{x}}=(-2.8,2.4)\,\mathrm{Mm}$, are ∼0.22 J m⁻², and here the energy density is divided 54%, 45%, and 1% between the acoustic, kinetic, and magnetic terms, respectively. Thus, the saturation level does a fairly accurate job at depicting the ratios of the energy terms at later times, after around 3.5t_N.

Cyan portions of the figure correspond to regions with roughly equal parts kinetic and magnetic energy densities (blue plus green), the partition of energy associated with fast magnetoacoustic waves in a low-β plasma and slow waves in a high-β plasma. In the same way, yellow regions have equal parts kinetic and acoustic energy densities and correspond to slow magnetoacoustic waves in the low-β case and fast waves in the high-β case. In the following, we will use "magnetic" or "acoustic" waves to indicate which of these energy density terms from Equation (14) dominates the other. In many circumstances, each type of wave for an ideal magnetohydrodynamic fluid exhibits equipartition of energy between its kinetic and potential terms (Zweibel & McKee 1995). Depending on the properties of the wave (standing or traveling, short or long wavelength) and the background medium (static, slowly varying, or turbulent fluctuations about an average equilibrium), equipartition may hold only weakly (in a spatially and/or temporally averaged sense), strongly (at each spacetime point), or not at all. Our simulation uses a high-frequency wave, and the ohmic dissipation term is only important in regions of strong gradients. As a result, for most locations in our simulation, we should find that the kinetic energy is in equipartition with the sum of acoustic and magnetic terms. We can think of allotting fractions of the kinetic term between the fast and slow wave in proportion to the magnetic and acoustic terms to get the total energy for each type of wave. This approach is only an approximation and will work best at locations where the two waves are strongly distinct (very high or low β regions). However, it will still prove useful in understanding the propagation of energy through the system. (For a discussion of important cases where strong equipartition does not hold, and implications for the interpretation of observations, see Goossens et al. 2013.)

Below the lower equipartition height, the wavepacket is predominantly acoustic (recall the color saturation, discussed above). The first several panels of Figure 6, up to ∼3.3t_N, show the slow upward propagation of the wavepacket. As the pulse crosses the lower equipartition layer, the wave energy mostly remains on the fast branch of the dispersion relation: it converts from an acoustically dominated (yellow) high-β fast wave to a magnetically dominated (cyan) low-β fast wave. The ongoing conversion is apparent in the division of energy at time 2.1t_N quoted above, resulting in 20% of the energy in magnetic perturbations.

The pulse rapidly expands in size as it reaches the higher phase speed portions of the domain, starting around 3.3t_N. Between 3.3t_N and 5.0t_N, the outer portions of the wavefront turn over and become downward propagating. The rapid expansion and the turnover are in agreement with the properties of the constant phase fronts from Figure 4. In the region surrounding the null, but still outside the magenta contour (e.g., still low β), a portion of the wavepacket refracts toward the null, and a portion refracts back downward toward the photosphere. Some of the initial wavefront passes around the null to reach the upper boundary and therefore escapes the system. These results again agree qualitatively with the results of the ray-tracing analysis of Section 3. Note that this is very different from what would happen if the magnetic field were uniform. Then, the phase speed always increases with height, and the fast wave will refract back toward the photosphere. The strong inhomogeneity in the field due to the presence of the null is responsible for the radically different wave behavior.

Around the null we see a continuous transfer of energy from the fast to the slow branch of the dispersion relation, and from magnetic to acoustic energy densities, starting when the wave approaches the c_s = v_A region surrounding the null, around 4.5t_N. This is visible in the animation of Figure 6 as the appearance of the yellow colored pulses propagating away from the null. They are first visible at t = 4.1t_N, in regions outside the equipartition layer where the plasma β is still less than, but close to, the equipartition value of β = 1.2. Energy transfer continues in the region spanning the equipartition contour. At least part of this behavior does not neatly fit into either of the currently used categories of transmission or mode conversion, as it involves, purely in the low-β region, a jump from one branch of the dispersion relation to the other and a change from magnetically to acoustically dominated energy density at the same time.

While the fast mode may propagate in any direction relative to the magnetic field lines, the slow mode is heavily guided along the field. In Figure 6, this shows up in the low-β regions as the yellow portions of wave energy density guided along the four separatrix field lines leaving the null. If one had a detector sitting on the separatrix at, say, ${\boldsymbol{x}}=[-2,3.5]\,\mathrm{Mm}$, one would see a fast magnetic wave headed toward the null, and later a slow acoustic wave headed away from the null. Because the magnetic field lines concentrate moving away from the null and back toward the photosphere, the energy leaving in the slow mode becomes increasingly concentrated around the separatrix field lines as it propagates downward toward the magnetic footpoints.

In summary, the simulation shows that the initial wavepacket splits up into numerous subpackets as it propagates. A small portion remains acoustic (transmits across the lower c_s = v_A boundary as a slow wave) and is confined to low-lying field lines in the low-β region. This accounts for ≈10% of the injected energy. Most of the initial packet remains a fast wave when crossing the lower equipartition height and subsequently refracts either left or right and returns to the lower boundary. However, a part of the upward-propagating fast wave is further refracted in toward the null point. This portion appears to mode-convert near the null and largely leaves the null region as slow acoustic waves confined to field lines near each of the null's four separatrices. The next two subsections will cover the mode conversion process near the null in more detail.

According to Schunker & Cally (2006), the amount of wave mode conversion depends on the angle of propagation relative to the magnetic field and the equipartition layer scale height H_e (see Equation (1)). In Section 3.2 we found that H_e ≈ 0.23 Mm for the equipartition layer surrounding the null. Because the equipartition curve is approximately circular with radius r_e = 0.73 Mm, we defined a radius of influence for the null as ${R}_{n}={r}_{e}+{H}_{e}\approx 1.3{r}_{e}=0.96\,\mathrm{Mm}$. These distances are indicated in Figure 7. Based on Equation (1), we expect mode conversion to become important as waves reach this distance from the null.

Zoom In Zoom Out Reset image size

The general description of mode conversion reported in the previous section is supported by combining the ray-tracing and numerical approaches. Figure 7 shows several curves in the vicinity of the null, with the separatrices shown in black. The blue vertical line located at x = −1 Mm passes through the initial center of the introduced wavepacket. The green curve is the path of a fast ray initialized at the lower boundary, with crosses showing equal phase distances (τ) or, equivalently, equal time differences in units of t_N, as the position of a phase point moves at the group velocity along the curve: ${\boldsymbol{x}}(\tau (t))={\int }_{0}^{t}{{\boldsymbol{v}}}_{g}({\boldsymbol{x}}({t}^{\prime }),{t}^{\prime }){{dt}}^{\prime }$, where ${{\boldsymbol{v}}}_{g}$ is defined below Equation (20). The red curve is the path of a slow ray, with crosses indicating a phase point moving at the slow speed; it closely follows a magnetic field line. The slow ray was initialized from a location along the fast-ray path, indicated by the triangle at τ ≈ 3.25t_N (${\boldsymbol{x}}(\tau )\approx [-0.49,3.11]\,\mathrm{Mm})$. This is where the fast ray reached one conversion scale height H_e away from the equipartition contour. The lengths of H_e and R_n are also indicated in the figure.

We extract the energy densities displayed in Figure 6 along the blue, green, and red curves shown in Figure 7 at each simulation output time to generate the time–distance plots shown in Figure 8. The color scheme and scaling are the same as in Figure 6. In each of the panels, the abscissa measures time in units of t_N, and the ordinate measures distance along the curve. For panel (a), distance is measured in Mm, but for panels (b) and (c), distance is measured in phase coordinates. For reference, the thick colored line at t = 0 in each panel indicates the line color of the line in Figure 7 along which the energy densities were extracted.

Zoom In Zoom Out Reset image size

Figure 8(a) is the time–distance diagram of energy density along the blue vertical line at x = −1 Mm of Figure 7. The initial shallow sloped portion at the bottom of the panel shows the upward-propagating initial pulse (magnetoacoustic fast wave), moving at the relatively slow sound speed below the equipartition height at y ∼ 0.5 Mm. Some of the energy remains acoustically dominated, switching from the fast to the slow branch of the dispersion relation. This behavior is visible as the shallow sloped green/red/yellow ridges beyond t = 5t_N, which shows that a portion of the wave refracts and becomes downward propagating around y = 1.5 Mm. However, most of the energy remains on the fast branch and converts to a magnetically dominated disturbance. This energy streams upward at a much faster speed starting around t = 3t_N, due to the increasing Alfvén speed. As the upward-propagating pulse approaches the separatrix layer, indicated by the horizontal dashed line, we again see two distinct portions. The cyan (magnetic, low-β fast-wave) portion passes across this topological barrier, while the yellow (acoustic, low-β slow-wave) portion is confined to field lines close to but underneath the separatrix.

Note that little of the behavior seen in panel (a) directly represents the physical energetics of the plasma. The velocities are mostly phase velocities caused by different portions of the wavefront refracting across the vertical slice at different times. In particular, the relation between the apparently upward-propagating magnetic energy and the acoustic energy confined near the separatrix is unclear. On the other hand, vertical integration of the different energy channels of panel (a) would most closely approximate line-of-sight observations. This illustrates the care required to interpret observations when one cannot assume that the field is nearly homogeneous.

In contrast, Figure 8(b) shows the energy densities as a function of time along the curved, fast-ray path shown in green in Figure 7. Distance is now measured as a phase through the time integral of the group velocity along the curve, the ${\boldsymbol{x}}(\tau (t))$ defined above. The white region at the top of the panel is due to this fast ray exiting the computational domain, and it is included to keep a consistent scaling between panels (b) and (c). The location of the start of the conversion region, R_N (the triangle here and in Figure 7), and the null point are each indicated by horizontal dashed lines. There are two immediately obvious and important differences compared to panel (a): (i) almost all the energy is confined to straight line paths with a slope equal to the local fast-wave group velocity (slope of 1 in the figure's units, diagonal dashed line), and (ii) almost none of the energy makes it beyond the null point. The second point is true even though energy propagation past the null is allowed as a result of the finite plasma β, because, unlike for a cold plasma, the group velocity of the fast ray does not drop to zero. Even so, it is evident that little energy actually follows this trajectory beyond the null.

Figure 8(c) shows the energy density time–distance plot for a hybrid curve of Figure 7, where the slow ray (red) has been stitched onto the fast ray (green) at the phase position where the fast ray crosses the interaction region surrounding the null (triangle). The lower portion of the plot, up to t ≈ 3.75t_N (horizontal dashed line), is the same for panels (b) and (c), with distance measured as phase position along the fast ray. Above this, distance is measured in phase position along the slow ray. Each portion is indicated by the red and green lines on the left of the plot. The magnetically dominated energy propagating along the fast ray transitions smoothly to acoustically dominated energy propagating along the slow ray, starting around t = 4t_N. In both cases the slope is 1 (compare again to the diagonal dashed line), indicating that the energy propagates at the fast group velocity in the bottom portion and at the slow group velocity in the top portion of this plot. Near the v_g label is a steep sloped region of magnetically dominated energy density, also starting around t = 4t_N. It is fully contained in the slow ray portion of the diagram and is due to a separate portion of the fast wavefront that refracts across the slow ray path. It has a high phase speed as it crosses the red path because the phase front is expanding and crosses the slow ray line at a high angle (compare with the animation of Figure 6). This is the outer portion of the fast-mode pulse that refracts back downward toward the photosphere.

The panels of Figure 8 only show the energy density along a select number of paths through the domain. However, the fact that the rate of propagation all along the fast-joined-to-slow ray is precisely the local group velocity is a strong indication that at least some energy follows this approximate path. This is, again, in sharp contrast to panel (a), which shows differences in propagation between the acoustically and magnetically dominated regions, and to panel (b), which shows energy propagating along the fast ray originally, then abruptly vanishing at the null point.

Energy densities extracted along other ray paths (not shown) all show similar results. For instance, energy densities extracted for the 31 rays traced in Figure 4 each show propagation of energy along the entire path, except those that pass within R_n of the center of the null equipartition region. Instead, the time–distance plots along those latter rays look similar to Figure 8(b). This suggests that wave energy enters the region around the null as a fast mode, but it does not exit the region as a fast mode. The energy must either convert to the slow mode or dissipate locally. We will now attempt to determine the fate of this energy.

We can estimate the total amount of mode conversion in the region surrounding the null by integrating each term of Equation (14). The energy inside a region ${ \mathcal A }$ is given by the integral of the energy density terms at each time,

$$\begin{eqnarray}&&W(t)={\int }_{{ \mathcal A }}({E}_{K}(t)+{E}_{{\rm{A}}}(t)+{E}_{M}(t))d{ \mathcal A },\end{eqnarray} \tag{ 23 }$$

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 }$

$$\begin{eqnarray}&&{W}_{\mathrm{flux}}(t)={\int }_{0}^{t}{\int }_{\partial { \mathcal A }}({{\boldsymbol{F}}}_{A}({t}^{\prime })+{{\boldsymbol{F}}}_{B}({t}^{\prime }))\cdot \hat{{\boldsymbol{n}}}{{dldt}}^{\prime }.\end{eqnarray} \tag{ 24 }$$

Plots of W(t) and W_flux(t) for two different areas are shown in Figure 9. In panel (a), ${ \mathcal A }$ is taken to be the circle of radius R_n centered on the centroid of the c_s = v_A region, resulting in W(t) displayed as the solid black line in Figure 9(a), normalized to the total wave energy introduced through the lower boundary, E_input. A fiducial at E/E_input = 0 is marked by the dotted line. We decompose the net flux into inward and outward terms, with inward magnetic (acoustic) shown as positive dashed blue (red) and outward shown as negative. About 15.5% of the injected wavepacket's total energy crosses the null's equipartition boundary (blue and red positive dashed curves). As expected, most of the energy arrives as a Poynting flux, consistent with the wave being predominantly a fast-mode wave in the β < 1 region. However, most of the energy leaving the null area is acoustic in nature and slightly lags the inward-directed Poynting flux, showing that mode conversion has taken place. The conversion is fairly efficient, with a ratio of 0.7 between the net acoustic flux and net Poynting flux across the boundary. Of the energy that passes into the conversion region, 7% appears to exit in a form consistent with how it entered. The remaining 23% is energy that enters the region but never exits: this energy heats the plasma through ohmic and shock dissipation, or is lost as a result of uncaptured numerical diffusion. We briefly describe the current accumulation at the end of this section.

Zoom In Zoom Out Reset image size

We can also calculate the energy that continues propagating upward into the corona by considering a region ${ \mathcal A }$ in the upper portion of our domain. In Figure 9(b), we show the results of computing the net flux across and energy above y = 5.6 Mm, which is above the transition region but below the upper damping region. Here, solid lines show energy density terms integrated over the entire domain above y = 5.6 Mm at each time, and dashed lines show the temporally integrated net fluxes across that boundary, up to each time. The downward fluxes are ≲1% of the upward fluxes, which is why we do not independently display them here, as we did in panel (a). Spatially integrated acoustic, kinetic, and magnetic energy densities are shown in solid red, green, and blue, respectively, while the Poynting and acoustic fluxes are shown in dashed blue and red, respectively.

As briefly discussed in Section 5.1, we do indeed find that the kinetic energy is typically in equipartition with the sum of the acoustic and magnetic energies, i.e., the solid green curve is roughly the sum of the solid red and blue curves. The exception is when the wavepackets begin to interact with the boundaries, causing the matched oscillations seen around times $t=(8-9){t}_{N}$ in Figure 9(b). The kinetic energy also appears proportionally allotted between the acoustic and magnetically dominated waves. In other words, each type of wave appears to be in equipartition, independently. For instance, the Poynting flux (dashed blue) accounts for both the magnetic energy density and its associated kinetic energy, so that the dashed blue line is essentially double the solid blue line; the same is true for the acoustic terms in red. Then, the second half of the dashed blue and red lines come from the allotted proportions of the solid green line. Because these are integrated fluxes, this plot only demonstrates equipartition between kinetic and potential energies in a spatially averaged sense.

The integrated flux across the boundary and energy above the boundary are two independent measures of the total wave energy above y = 5.6 Mm. These are shown in black and match very well until they begin to diverge around t = 6t_N. This is expected, as at that point the waves reach the damping regions at the edges of our computational domain and eventually begin to exit the domain through the upper boundary. Because the Alfvén and sound speeds are fairly uniform above this height, the waves undergo little additional refraction. The total integrated flux across the height y = 5.6 Mm is thus the total proportion of energy from the initial wavepacket that we estimate continues propagating upward into the model corona. This accounts for 11% of the initial energy of the wavepacket, with 3.3% associated with a slow magnetoacoustic wave and 7.7% with a fast magnetoacoustic wave.

Essentially all of the acoustic flux across y = 5.6 Mm originated from mode conversion at the null point. Figure 9(a) shows that ≈11% of the total energy left the null point region as an acoustic wave. It appears that, at least for this simulation, the converted energy that leaves the null is roughly equally distributed along each of the null's four separatrices, with perhaps a slight bias toward the upward leg. The animation of Figure 6 appears to qualitatively support this conclusion.

On the other hand, the Poynting flux leaving the null accounts for just ≈1% of the total energy, while 7.7% was found to cross y = 5.6 Mm. At least 6.7% must therefore be due to a portion of the fast-mode wave that refracts around the null but does not enter the null's mode conversion region. This would be the portion of the wavefront discussed in Section 3, such as rays initiated between the green and yellow curves of Figure 3(b). That figure only shows rays initiated with k_x = 0. An analysis similar to one resulting in Figure 5 could be made to determine the distribution of rays initiated at the lower boundary that exit the top of the system. We will not perform that analysis here, but simply note that a structured magnetic field can allow some fast-wave energy to propagate into the corona, whereas for a uniform field it would eventually refract back to the photosphere.

In the previous section we found that 23% of the energy that enters the null point region is never measured to exit it. This is likely due to the localization of currents to the null point and separatrices and the subsequent dissipation of those currents. Figure 10 shows a time sequence of the out-of-plane current density (gray scale) around the null. White shows positive (out of page) and black negative (into page) currents. Blue lines are magnetic field lines, red lines are the approximate locations of the separatrix field lines, and the magenta curve again marks the c_s = v_A equipartition contour. The gray scale is saturated at ±5 × 10⁻³ A m⁻². The maximum current density is $| {j}_{z}| \approx 0.016\,{\rm{A}}\,{{\rm{m}}}^{-2}$ at t = 5.5t_N, marked in the figure by the green cross (this is the greatest current density for a larger field of view than shown in the figure, as well).

Zoom In Zoom Out Reset image size

The first pulse of the wavepacket passes through the null point region around 4t_N. This stresses the field and leads to regions of current density localized along each separatrix, forming a set of current ribbons by 4.6t_N. The second pulse of the wavepacket steepens to form a shock as it approaches the null, around 5.2t_N. This appears to be the fast oblique magnetic shock reported by McLaughlin et al. (2009) and is accompanied by a strong spike in the current density, as found by Gruszecki et al. (2011) and Afanasyev & Uralov (2012) and references in those papers. We will defer details of shock formation and the development of the current sheet to a later paper.

The null point collapses to form a current sheet originally oriented at ∼45° to the separatrices. Next, a further accumulation of current density is evident in the ribbons that form along each of the separatrices. The transverse (cross-field) length scale of each ribbon decreases until it reaches the diffusion scale, ℓ = L_N/S, at which point the ribbons ohmically dissipate. This is most easily seen for t > 5.8t_N. The reduction in current density is not associated with plasma flows or energy fluxes. This dissipation accounts for the 23% of the energy that enters the equipartition region but does not leave, quoted above.

The current sheet later collapses again at −45° to the separatrices, rotated by 90° to the original collapse direction. The process repeats itself, with apparent oscillatory reconnection at the null. This type of behavior has been observed in other studies (McLaughlin et al. 2009; Murray et al. 2009; McLaughlin et al. 2012a, 2012b). A more detailed discussion of the formation of the current sheet, its dissipation, the oscillatory reconnection, and how each depends on the initial properties of the wavepacket will be presented in a future work.

A natural question is how often one expects a low-lying null, so critical to our investigation, to arise on the Sun. Current estimates for the number of nulls, and the strong gradients in the Alfvén speed that come with them, give roughly one per supergranular cell above 1.5 Mm (Close et al. 2004; Régnier et al. 2008; Longcope et al. 2009), and likely more below that height. Longcope et al. (2009), for instance, used a spectral method to determine that about 19,000 nulls with a height above 1.5 Mm are present as a result of the quiet-Sun field at any given time, about 1 null per 300 Mm². Roughly half of these occur between 1.5 and 4.5 Mm, and the statistics vary only by some 10% over the two solar minima and 550 magnetograms considered. Since high moments drop off more rapidly with height, they expect more nulls below 1.5 Mm, but the noise threshold of MDI magnetograms, the basis of their analysis, prevents estimates below that height. Freed et al. (2015) provided an observational study of the corona using AIA data to determine the distribution of nulls with height, and found rough agreement with the work cited above, for the heights they were able to observe.

The presence of many nulls between the photosphere and the transition region makes it likely that the physical processes we describe in this paper are rather common, namely, that convectively driven waves pass through multiple conversion sites as they propagate. This idea is supported by other simulations. Nutto et al. (2012) simulated a convectively unstable network region that did not contain, a priori, a null point. They found that multiple overlaying equipartition regions form self-consistently. Some even have the roughly circular shape that may indicate a null, though the authors did not specifically check for the presence of one (see their Figure 3, panel t = 200 s). They found evidence for multiple conversions between fast and slow, acoustic and magnetic pulses at each of the equipartition regions, and found that a slow-mode acoustic wave is able to propagate outward beyond the chromosphere, as did we.

We therefore expect the processes we describe in this paper to be applicable to regions of plage, network, and quiet Sun. We anticipate that the relative importance of refraction and mode conversion will vary based on the relevant length scales: the height of the null, distance between magnetic concentrations, and the size of the null's equipartition region. Future studies are required to determine which effects dominate in different atmospheric regions. The present simulation corresponds most nearly to a 2D cut through a small ephemeral region or strongly enhanced network flux.

Active regions present a different challenge. The wave fields in quiescent active regions with highly symmetric umbrae and penumbrae are likely well described by slowly varying sunspot models, as evidenced by the good agreement between observations and simulations of acoustic halos (Rajaguru et al. 2013; Rijs et al. 2016). Strong gradients in the magnetic field in general, and null points in particular, probably play a more important role in active regions with more complicated structure. Nulls are found above the photosphere in some active regions, but not the majority; at the same time, active regions with nulls are much more likely to flare than those without (Barnes 2007). Coronal nulls have separators connected to them, and photospheric driving will localize currents to these separators (Longcope 2001; Parnell et al. 2008). Focusing of wave energy at these locations may provide a way to destabilize these current sheets, accounting for the increased rate of flares, or even sympathetic eruptions (Gruszecki et al. 2011). Energy leaving the reconnection site as a slow mode will be guided and focused along the separatrices to the photospheric footpoints, the flare ribbons, which we will briefly consider below.

Our analysis is most closely related to other studies of MHD waves near nulls, and here we review those results related to our own. The key properties for comparison are that our magnetic field is nonlinear, our plasma has finite β, and our density is nonuniform. The most well-studied situation is a linear null for which β = 0 and the density is uniform. In that case, wave energy accumulates at the null because of refraction, as seen, for instance, by the inward spiral of rays (McLaughlin & Hood 2004). Focusing of the wave causes an enhancement in current density at the null, which then dissipates to locally heat the plasma. By repeated reflections between the null and an external conducting boundary, all of a wave's initial energy is dissipated at the linear null (Craig & McClymont 1991; Hassam 1992).

When the field around the null is nonlinear, saddle points and extrema in the phase speed will cause the phase front of an incoming plane wave to split, as reported by McLaughlin & Hood (2006a) and seen in our work for the rays (Figure 3) and simulation output (Figure 6). Longcope & Tarr (2012) determined the energetic consequences of this for a β = 0, uniform-density, quadrupolar field and found that only ∼40% of the energy dissipated at the null, while the rest propagated away to infinity as a fast wave. Note that those authors studied a wave initiated at the null, rather than one initially propagating toward the null; later reflections carried only part of the wave back to the null. The difference between the quadrupole and linear null is the introduction of a new length scale, the distance over which the linear term in the Taylor series expansion of the magnetic field about the quadrupole null is valid. An incoming wave that is planar on that scale will refract toward the null, and the rays describing it will form the usual logarithmic inspirals. Waves on larger scales will split and partially refract away, taking energy with them.

Finite β introduces the slow-mode wave and allows a coupling between the fast and slow modes. McLaughlin & Hood (2006b) is the only other study we are aware of that attempted to quantify conversion near a null, so it is worth discussing their work in detail. They studied a linear null and varied a parameter β₀ that set the size of the equipartition region surrounding their null. As they point out, because the field of the linear null has no inherent length scale, by changing β₀ they effectively changed the distance between equipartition layer and the initial fast plane wave launched from the boundary. When β₀ ≪ 1, the equipartition region is small compared to the extent of the wavefront (or size of the computational box), and this mimics the β = 0 case: the wave refracts and the rays spiral inward. As β₀ increases, the size of the equipartition region increases, and the wave is initiated closer to it. Once β₀ is large enough, c_s > v_A everywhere in their domain, which corresponds to the wave being initiated inside the equipartition layer.

To estimate the amount of conversion, McLaughlin & Hood (2006b) initiated a strictly fast-wave pulse in the β > 0 region and noted where it completely separated into two pulses (fast and slow) within the c_s = v_A contour. They calculated the integral of the perpendicular velocity inside the slow-wave pulse, which is a quantity related to the wave momentum (assuming that the wave is completely polarized transverse to the field). They found the ratio of this quantity to the integral of the perpendicular velocity over the initial pulse for different values of β₀. The ratio is roughly proportional to β₀. This is because increasing β₀ effectively initiates the wave closer to the null, so that a larger fraction of the wave refracts into the equipartition region where it can convert. They could therefore identify the competing effects of refraction and conversion, but not study them further.

The most conversion McLaughlin & Hood (2006b) report is ∼25% for β₀ = 10, which gives the maximum size their equipartition region can be while still fitting inside the domain. A peak conversion for this case is somewhat surprising, since conversion should begin roughly one scale height H_e outside the equipartition boundary, and their initial pulse is inside that region. It is possible that their use of the perpendicular velocity is no longer a good approximation for this case. Their estimate is difficult to directly compare to our own because they used a momentum-related quantity and only calculated it when the waves entered the c_s = v_A boundary, not when they exited. In contrast, our estimate uses wave energy and relates the total incoming to total outgoing energy of each type. We feel that our approach is more physical and allows us to quantify local dissipation. Considering their results for different β₀ suggests that we should perform a set of our own experiments in which the height of the transition region is varied while the magnetic field and wavepacket parameters are held fixed. This will change the size of the equipartition region at the null, but not the total initial wave energy, and can therefore be used to determine the relative importance of refraction and mode conversion for a given driver.

Refraction clearly dominates when β < 1, and its dominance in much of the corona has led multiple authors (McClements et al. 2004; McLaughlin et al. 2009; Threlfall et al. 2012), including McLaughlin & Hood (2006b), to argue that an azimuthally symmetric pulse will be a good approximation for waves approaching nulls. However, we did not find this to be the case in our simulation, as can be seen around t = 4t_N in Figure 6 and the animation of that figure. This is because our wavefront was of comparable size to the null's equipartition region, and because the density and temperature are not uniform around the null, so that the wave fronts do not simply spiral inward. For nulls high in the corona where the scale height is large and the equipartition contour is small, the approximation is more appropriate. On the other hand, based on the statistics reviewed in Section 7.1, we expect most nulls on the Sun to arise in locations where stratification is important.

From the discussion in this section, it seems that there are two important length scales. (i) The first is global and is caused by local extrema in the wave phase speeds. This will determine the spatial portion of an arbitrary wave that will refract in toward a null. (ii) The second is the size of the equipartition region around the null relative to the spatial extent of the wavefront that refracts toward the null. The second scale will affect mode conversion, current accumulation, and dissipation at the null. Combining our results with those of McLaughlin & Hood (2006b) shows that the combination of these two lengths will determine how much wave energy refracts toward a null, and how much of that mode converts at the null, for a given driver. They can be combined to define an overall "region of influence" for nulls.

In our simulation, we found that slow-mode energy leaving the null is concentrated near each separatrix field line. Figure 9(b) shows that ∼30% of the slow-mode energy leaving the null travels upward, so the remaining 70% propagates back to the photosphere and becomes increasingly concentrated around each of the other three separatrices. This is a focusing effect due to the converging magnetic field that acts as a waveguide for the slow waves. In total, we found that 8% of the initial upward-propagating wavepacket becomes concentrated in the three downward-propagating patches. This study therefore shows that mode conversion near the null can take a distributed upward wave flux and create localized patches of downward wave flux, focused specifically on the separatrix footpoints.

A similar effect has been reported by Russell & Stackhouse (2013) for Alfvén waves in a zero β plasma. They introduced a transverse velocity perturbation at the apex of a model arcade in order to mimic a reconnection event. The resulting pure Alfvén waves focus as they propagate downward, mostly due to the converging field lines, although phase mixing also plays a role. In contrast, when they introduced a fast wave above an arcade, they found that the energy density became more diffuse as it propagated toward their lower boundary.

In our case, the fast waves that convert near the null are not necessarily generated by the convective processes that we modeled. For instance, a reconnection-induced fast wave will be refracted toward any nulls in the chromosphere or corona. Near the nulls, fast-mode energy will partially convert to a set of slow modes, each of which will focus on the null's separatrix footpoints. The effect is especially important in the partially ionized chromosphere, which modifies the slow-mode propagation and introduces frequency-dependent damping (Soler et al. 2013). Focused Alfvén and slow waves could therefore contribute to the enhanced emission seen in flare ribbons. This has been argued before (Nakariakov & Zimovets 2011), and our simulations at least show how focused wave energy could arrive at flare footpoints. The topic seems to deserve further consideration.

We chose to isolate and analyze the dynamics of a single wave pulse. This approach makes direct comparison with most other studies on atmospheric waves difficult, as they typically rely on time-averaged properties of long-lasting wave trains or stochastic fluctuations (Fedun et al. 2011a; Santamaria et al. 2015; Rijs et al. 2016). Some more recent efforts have included either short or instantaneous pulses (Khomenko & Collados 2009; Nutto et al. 2012; Moradi et al. 2015; Rijs et al. 2015; Santamaria et al. 2015), but were analyzed by Fourier transform of, say, the velocity signal, which is still a time-averaged approach. An exception is Shelyag et al. (2016), whose primary focus was dissipation due to ambipolar diffusion in a small, 3D flux tube. They introduced a similar perturbation to ours and found that it converts to fast and slow waves, as well as Alfvénic waves due to the 3D geometry. They found strong heating due to dissipation of currents associated with the Alfvénic waves, but did not attempt to quantify the conversion process.

In a complementary simulation to our own, Santamaria et al. (2015) studied low-frequency (3.3–5 mHz) waves in a stratified atmosphere spanning the upper convection zone to low corona that also included a null point. They considered both horizontally and vertically driven waves. For the vertically driven case they found little time-averaged Poynting flux above the height of the transition region, but did find significant low-β acoustic flux confined to nearly vertical flux tubes at the edge of their domain. That mechanism is distinct from the one we describe, where a fast mode propagates toward a null, converts at the equipartition region, and exits the null point as a slow acoustic wave confined to the separatrix field lines. However, in their Section 3.3.1 they appear to mention the mode conversion process near the null that we describe in much more detail. In their case, conversion arises owing to horizontal driving at the lower boundary. Yet, by the time the wave reaches the photospheric level, it contains substantial vertical oscillations in the regions of the strongly inclined (near horizontal) field and therefore produces a similar type of driving to what we have implemented. Indeed, one can see that some acoustic power is localized near the separatrices in their Figure 9, which would indicate the presence of slow-mode waves generated by conversion near the null. Their simulations therefore suggest that our results will carry over to lower-frequency waves, closer to the Brünt–Väisälä frequency, but we caution that at those frequencies dispersion will play a strong role that we did not consider in our work.

A limitation of the present study is the 2D geometry. A 3D simulation will allow for the coupling of the Alfvén mode to the fast and slow modes discussed here, as Shelyag et al. (2016) found. Cally & Goossens (2008) and Felipe (2012) demonstrated that the Alfvén mode coupling will show a second strong dependence on the relative orientation of the wavevector and the magnetic and gravitational fields. These studies used slowly varying (or simply uniform) sunspot fields with no strong gradients.

For null point studies, McLaughlin et al. (2008) performed a WKB analysis for a linear 3D null and found that fast-wave energy accumulates at the null as a current density. Thurgood & McLaughlin (2012) performed a 3D numerical simulation and identified a coupling between the Alfvén and fast mode, though they did not attempt to quantify the amount of conversion. Both of these 3D null point studies had β = 0, so that the only form energy accumulation could take is an increase in current density as the wave contracts and the gradient across the wavefront increases. Pontin et al. (2013) performed a $\beta \ne 0$, 3D simulation (see next section) and found that currents localized to a null and separatrix surface, but did not describe the dynamics in terms of waves.

Determining the proportions of energy in each of the three modes is nontrivial in 3D. Our 2D geometry allowed us to identify with relative ease the fast and slow modes based on their energy densities and propagation with respect to the field. Distinguishing between the fast and Alfvén mode will require more care and, we anticipate, either a greater reliance on the spacetime diagrams like Figure 8 or building up a statistical version of the local dispersion relation by studying oscillations of local parameters. The latter approach has been used to analyze observational data (Tomczyk & McIntosh 2009) but could be applied in 3D to simulation output. A third approach is to project the energy density and flux terms onto characteristic directions of the magnetic field, as Felipe (2012) has done for the 3D sunspot simulation, Thurgood & McLaughlin (2012) did for a 3D numerical simulation, and McLaughlin et al. (2016) did for a 3D null in a WKB ray-tracing investigation. Mumford et al. (2015) had some success using this method to determine the energy flux associated with several modes of a 3D expanding flux tube, but the presence of multiple propagation velocities at a given height in their phase diagrams (see, e.g., their Figure 7(a)) suggests that the decomposition is not exact. It is clear that no single method can currently robustly distinguish between modes in 3D. The best approach is probably to combine or compare multiple types of analysis, as we have done in 2D in the present investigation.

In agreement with other studies, our simulations show that currents accumulate at topologically important locations: the null and the separatrices. Pontin et al. (2013) considered a 3D magnetic dome topology, similar to what would result if our initial magnetic field was spun around the axis passing through the null. In that case, the two separatrices leaving our 2D null in the vertical direction become 3D spine lines, while the other two separatrices form a single dome-like separatrix surface. The initial condition of Pontin et al. (2013) had a uniform density and temperature, as opposed to our stratified atmosphere. They applied an incompressible shear flow around the spine footpoint inside the dome and observed an accumulation of current at the separatrix and around the null (see their Figure 9). Other work on 3D nulls has also found that incompressible perturbations cause current accumulation (Pontin & Galsgaard 2007), where the details depend on whether the area around the spine or the fan is perturbed, and whether the perturbation is of rotational or shear type.

Our simulations show that current accumulation also occurs for compressive waves, even when β > 0 and waves can travel through the null. This situation is not well studied. McLaughlin et al. (2009) provide one of the most detailed analyses to date, but set β = 0 in their initial condition. They rely on current accumulation, shock formation, and the resulting ohmic and shock dissipation to heat the plasma and raise β above zero. In their case, asymmetry in the heating about the null creates a pressure gradient that must be balanced by a Lorentz force, hence the persistence of the current density at the null after the wave has passed. Mode conversion may play a role in their simulation, though they do not consider it in their paper. McLaughlin & Hood (2006b) (discussed at length above) found that most current accumulates close to a linear null, but their simulation solved the linearized MHD equations, whereas the nonlinear effects are clearly important. Afanasyev & Uralov (2012) demonstrated that nonlinear effects substantially alter wave behavior near the null and allow fast waves to pass even for the β = 0 case. They studied both β = 0 and finite β cases and used a modified WKB method that accounts for weak shocks. Their method does not apply when one shock is downstream of another, e.g., when rays cross, which prevented them from studying current sheet formation. However, we demonstrated in Figures 3–5 that this occurs frequently. Gruszecki et al. (2011) also studied the nonlinear evolution of a fast wave near a linear null and for an initially β = 0 plasma. They noted the shock formation and strong cospatial spikes in the current density, but stopped their simulations before the pulse reached the null in order to stay in the low-β regime. They were thus unable to note whether any current persisted at the null or along its separatrices. So it is clear that more work on how currents localize to nulls, particularly for finite β, is necessary.

We found that the current sheet that forms at the null oscillates between ±45°, which is evidence for oscillatory reconnection (see Figure 10). This finding is also not new (Craig & McClymont 1991), but has recently received renewed attention, especially in connection with observations of quasi-periodic pulsations associated with flares (Van Doorsselaere et al. 2016). Oscillatory reconnection is found in the null point studies (McLaughlin et al. 2009, 2012a), in simulations of flux emergence in a stratified atmosphere (Murray et al. 2009; McLaughlin et al. 2012b), and now in our simulation modeling convection-induced waves. In short, oscillatory reconnection appears to be a rather robust feature of reconnection, at least in 2D, and it remains to determine its importance in 3D.