On the construction of general large-amplitude spherically polarised Alfvén waves

Abstract

In a magnetised plasma on scales well above ion kinetic scales, any constant-magnitude magnetic field, accompanied by parallel Alfvénic flows, forms a nonlinear solution in an isobaric, constant-density background. These structures, which are also known as spherically polarised Alfvén waves, are observed ubiquitously in the solar wind, presumably created by the growth of small-amplitude fluctuations as they propagate outwards in the corona. Here, we present a computational method to construct such solutions of arbitrary amplitude with general multidimensional structure, and explore some of their properties. The difficulty lies in computing a zero-divergence, constant-magnitude magnetic field, which leaves a single, quasi-free function to define the solution, while requiring strong constraints on any individual component of the field. Motivated by the physical process of wave growth in the solar wind, our method circumvents this issue by starting from low-amplitude Alfvénic fluctuations dominated by a strong mean field, then ‘growing’ magnetic perturbations into the large-amplitude regime. We present example solutions with non-trivial structure in one, two and three dimensions, demonstrating a clear tendency to form very sharp gradients or discontinuities, unless the solution is one-dimensional. As well as being useful as an input for other calculations, particularly the study of parametric decay, our results provide a natural explanation for the extremely sharp field discontinuities observed across magnetic field switchbacks in the low solar wind.

1. Introduction

The existence of the incompressible Alfvén wave is perhaps the most important distinguishing feature between the dynamics of neutral fluids and plasmas, underlying key physics of turbulent energy dissipation and enabling non-trivial dynamics well below the scale of the mean free path (e.g. Schekochihin et al. 2009). Interestingly – and unlike other wave-like perturbations of fluids or plasmas – Alfvén waves have an unambiguous nonlinear generalisation (Barnes & Hollweg 1974): defining the plasma density $\rho$, pressure $P$, magnetic field $\boldsymbol {B}$, flow velocity $\boldsymbol {u}$ and field strength $B=|\boldsymbol {B}|$, any perturbation that satisfies

(1.1a–d)\begin{equation} P = {\rm const.},\quad\rho= {\rm const.},\quad B^{2} = {\rm const.},\quad \delta\boldsymbol{u} ={\pm} \delta \boldsymbol{B}/\sqrt{4{\rm \pi} \rho} \end{equation}

is a nonlinear solution that propagates in the direction $\hat {\boldsymbol {b}} = \boldsymbol {B}/B$ at speed ${{v}_{A}} \equiv |\overline {\boldsymbol {B}}|/\sqrt {4\pi \rho }$ (where the overline signifies a spatial average and $\delta \boldsymbol {B} =\boldsymbol {B}-\overline {\boldsymbol {B}}$). Such solutions are necessarily maximally ‘imbalanced’ (the magnitude of the fluctuation's cross-helicity equals its energy), and, in the wave frame at speed ${{v}_{A}}$, involve constant plasma kinetic energy (constant $|\boldsymbol {u}|$) and zero motional electric field $\boldsymbol {u}\times \boldsymbol {B}$ (Matteini et al. 2015). Moreover, unlike large-amplitude sound or magnetosonic waves, spherically polarised perturbations do not steepen into shocks even if $|\delta \boldsymbol {B}|\gg |\overline {\boldsymbol {B}}|$, although they can be unstable (Sagdeev & Galeev 1969; Cohen & Dewar 1974). They also share many properties of linear (small-amplitude) Alfvénic fluctuations, such as how they change in amplitude and/or refract in a non-homogenous background medium (e.g. Barnes & Hollweg 1974; Hollweg 1974). The solution (1.1a–d) is even valid in a collisionless plasma as well as in collisional magnetohydrodynamics (MHD), so long as $\delta \boldsymbol {B}$ varies over scales much larger than the ion gyroradius or skin depth (Barnes & Suffolk 1971; Kulsrud 1983). Given these properties, it is perhaps not surprising that perturbations close to (1.1a–d) are observed ubiquitously in our best-studied example of a natural plasma – the solar wind (Belcher & Davis 1971). In this context, they are known as spherically polarised and their properties likely underlie a number of key aspects of solar-wind physics, including heating and acceleration (e.g. Bale et al. 2021; Shoda, Chandran & Cranmer 2021), turbulent spectra (e.g. Matteini et al. 2018; Bowen et al. 2021) and properties of ‘switchback’ field reversals (Kasper et al. 2019; Squire, Chandran & Meyrand 2020; Johnston et al. 2022).

Despite these important features, to our knowledge, there does not exist any general method to construct and study such nonlinear Alfvénic solutions when $\delta \boldsymbol {B}\gtrsim |\boldsymbol {B}|$, or even any general results regarding their existence. Previous methods (Valentini et al. 2019; see also Roberts 2012), though promising, may lead to unnecessary discontinuities in the solutions, the causes of which are discussed below. The difficulty arises from the twin constraints of $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = 0$ and $B=|\boldsymbol {B}|={\rm const.}$, which leave only one degree of freedom from which to construct the three-dimensional (3-D) vector field $\boldsymbol {B}$. Moreover, as we show below, this degree of freedom cannot in general be freely chosen when $\delta \boldsymbol {B}$ approaches $\overline {\boldsymbol {B}}$ in magnitude, even though simple, smooth solutions do exist (they just require non-trivial constraints on all components of $\boldsymbol {B}$ simultaneously). It is the purpose of this contribution to present such a method and briefly explore the properties of the solutions (1.1a–d). These are of interest for two reasons. First, the method can be used as input or initial conditions for other numerical calculations, such as the study of parametric decay (Del Zanna 2001; Primavera et al. 2019) or large-amplitude reflection-driven turbulence (Squire et al. 2020; Johnston et al. 2022). Second, the solutions are interesting in and of themselves – our method shares strong similarities with the plasma-expansion process that generates large-amplitude Alfvénic states in the solar wind, so the structures and properties that arise may be of direct physical relevance. Of particular interest is the development of sharp field discontinuities, which we show form much more readily two-dimensional (2-D) or 3-D solutions than one-dimensional (1-D) ones; this is promising, since highly discontinuous fields seem to be observed in switchbacks by Parker Solar Probe and other spacecraft (Bale et al. 2019; Akhavan-Tafti et al. 2021). More generally, the method provides an extremely broad class of nonlinear solutions to the MHD equations that complements other solutions such as force-free magnetostatic equilibria (Marsh 1996). Seen differently, it provides a way to construct a general divergence-free, unit vector field.

Below, we first discuss the basic idea of the method and equations involved, then present its application to magnetic field configurations that vary in one, two and three dimensions. The 1-D version, although idealised and explored previously in Mallet et al. (2021) (hereafter M$+$21) and Squire et al. (2022), provides a helpful illustration of the method and the challenges involved in higher dimensions. We finish with a brief discussion of some interesting features of the solutions, focused particularly on the generation of sharp gradient and/or discontinuities, as well as some possible applications of our results.

2. Method

The method we propose is based on splitting the magnetic field into its mean ($\overline {\boldsymbol {B}}$) and fluctuating ($\delta \boldsymbol {B}$) parts, then devising an equation to grow $\delta \boldsymbol {B}$ in amplitude while maintaining $\overline {\delta \boldsymbol {B}}=0$, $\boldsymbol {\nabla }\boldsymbol {\cdot }\delta \boldsymbol {B}=0$, fixed $\overline {\boldsymbol {B}}$ and constant $B^{2} = |\overline {\boldsymbol {B}} + \delta \boldsymbol {B}|^{2}$. In this way, the method ‘grows’ a small-amplitude spherically polarised wave, which can be relatively easily constructed from linear Alfvénic perturbations via standard optimisation methods, to any desired amplitude $\mathcal {A}\equiv (\overline {\delta \boldsymbol {B}^{2}}/\overline {\boldsymbol {B}}^{2})^{1/2}$. The obvious candidate to evolve $\delta \boldsymbol {B}$ is simply the induction equation supplemented by exponential growth,

(2.1)\begin{equation} \frac{\partial}{\partial t} \delta \boldsymbol{B} = \delta \boldsymbol{B} + \boldsymbol{\nabla}\times (\tilde{\boldsymbol{u}}\times \delta \boldsymbol{B}),\end{equation}

which clearly maintains $\overline {\delta \boldsymbol {B}}=0$ and $\boldsymbol {\nabla }\boldsymbol {\cdot }\delta \boldsymbol {B}=0$ by construction. The flow $\tilde {\boldsymbol {u}}$ should be non-Alfvénic: this is not the flow associated with the Alfvén wave itself ($\delta \boldsymbol {u}$ in (1.1a–d)), but rather the flow needed to change the shape of $\delta \boldsymbol {B}$ as it grows with $t$, in order to maintain constant $B$. The clear choice is the potential/compressive flow $\tilde {\boldsymbol {u}} = \boldsymbol {\nabla }\phi$, which cannot contribute to the Alfvénic flow on a periodic domain because $\boldsymbol {\nabla }\boldsymbol {\cdot }\tilde {\boldsymbol {u}}=\nabla ^{2}\phi$. To proceed, we form ${\partial }/{\partial t}|\overline {\boldsymbol {B}} + \delta \boldsymbol {B}|^{2} = 2(\overline {\boldsymbol {B}} + \delta \boldsymbol {B})\boldsymbol {\cdot } {\partial \delta \boldsymbol {B}}/{\partial t}$ and assert that its fluctuating part must be zero, i.e. that $\phi$ should be chosen so that $\partial /\partial t\,\delta ( B^{2})= 0$. Using the fact that $\boldsymbol {\nabla } B^{2}=0$ and $\boldsymbol {\nabla } \overline {\boldsymbol {B}} = 0$, this yields

(2.2)\begin{equation} B^{2} \boldsymbol{\nabla}_{{\perp}}^{2}\phi = B^{2}\nabla^{2}\phi - \sum_{ij}B_{i}B_{j}\boldsymbol{\nabla}_{i}\boldsymbol{\nabla}_{j}\phi ={-}\delta \boldsymbol{B}\boldsymbol{\cdot}\overline{\boldsymbol{B}}, \end{equation}

where $B_{i} = \overline{B}_{i} + \delta B_{i}$. If (2.2) can be solved at each step, with the result used to evolve $\delta \boldsymbol {B}$ through (2.1), the total field will maintain spatially constant $B^{2}$ (and $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = 0$) as $\delta \boldsymbol {B}$ grows to $\mathcal {A}\gg 1$.

This method was motivated by the reduced equations of M$+$21, which describe the evolution of 1-D spherically polarised Alfvénic states in the expanding solar-wind plasma. In this case, the amplitudes of the fluctuating field and (radial) mean field evolve as $\propto a^{-1/2}$ and $\propto a^{-1}$, respectively, where $a$ is the plasma's expansion factor. This leads to a similar situation where $\delta \boldsymbol {B}$ grows compared to $\overline {\boldsymbol {B}}$. Expansion also causes the gradient operators in (2.1) to differ in the perpendicular and parallel directions and change with $a$, as well as rotating the mean field if it has non-radial components. These effects necessitate extra terms in (2.1) and (2.2) (see § 3 for details), but do not fundamentally modify the method. Such terms cause intriguingly different nonlinear solutions to develop (see, e.g. figure 4 of Squire et al. (2022)) pointing to interesting generalisations of our method, such as making the direction of $\overline {\boldsymbol {B}}$ change in time or modifying the gradient operators to capture the physical effects of super-radial expansion and wind acceleration relevant to the inner heliosphere (Tenerani & Velli 2017). The correspondence with M$+$21 also shows that, while (2.1) and (2.2) are formulated on a purely mathematical basis, the process of $\delta \boldsymbol {B}$ growth does share important similarities with the real physical processes that generate large-amplitude spherically polarised states in the solar wind.

3. One dimension

In 1-D solutions, $\delta \boldsymbol {B}$ varies only in the $\hat {\boldsymbol {p}}$ direction, which can be at an arbitrary angle to $\overline {\boldsymbol {B}}$. This case is significantly more straightforward than 2-D or 3-D solutions and serves to illustrate some useful points. Denoting the coordinate in the $\hat {\boldsymbol {p}}$ direction as $\lambda$, such that $\boldsymbol {\nabla }\phi = \hat {\boldsymbol {p}} \,{\rm d}\phi /{\rm d}\lambda = \hat {\boldsymbol {p}}\phi '$, where $\hat {\boldsymbol {p}}$ is a unit vector, (2.1) and (2.2) simplify to

(3.1a,b)\begin{equation} \frac{\partial\delta \boldsymbol{B}}{\partial t}= \delta \boldsymbol{B} - \frac{\partial}{\partial \lambda}\left[\phi' (\delta \boldsymbol{B} + \overline{\boldsymbol{B}}_{{T}})\right],\quad |\delta \boldsymbol{B} + \overline{\boldsymbol{B}}_{{T}}|^{2} \phi'' ={-} \delta \boldsymbol{B}\boldsymbol{\cdot}\overline{\boldsymbol{B}}.\end{equation}

Here $\overline {\boldsymbol {B}}_{{T}} = \overline {\boldsymbol {B}} - \hat {\boldsymbol {p}} (\hat {\boldsymbol {p}}\boldsymbol {\cdot }\overline {\boldsymbol {B}})$ is the part of $\overline {\boldsymbol {B}}$ that is transverse to the mean field. We see a clear correspondence with (59) and (61) of M$+$21 (with $\overline {\boldsymbol {B}}_{{T}}$ denoted $\boldsymbol {v}_{\rm AT}$), which can actually be equivalently derived by asserting that $\partial /\partial t\,\delta (B^{2})=0$ (i.e. in the same way as (2.1) and (2.2)), as opposed to the asymptotic expansion in slow expansion rate used therein. The extra terms in M$+$21 result from expansion-induced rotation of $\hat {\boldsymbol {p}}$ and $\overline {\boldsymbol {B}}$, as mentioned above.

Solving (3.1a,b) first requires an initial condition with constant $B^{2}$. This is easily constructed by arbitrarily specifying the component of $\delta \boldsymbol {B}$ in the $\hat {\boldsymbol {n}} \equiv \hat {\boldsymbol {p}}\times \overline {\boldsymbol {B}}/|\hat {\boldsymbol {p}}\times \overline {\boldsymbol {B}}|$ direction – i.e. the $\delta \boldsymbol {B}$ direction for a linear Alfvén wave – then solving for the $\delta \boldsymbol {B}$ component in the $\hat {\boldsymbol {m}}\equiv \hat {\boldsymbol {p}}\times \hat {\boldsymbol {n}}/|\hat {\boldsymbol {p}}\times \hat {\boldsymbol {n}}|$ direction, noting also that $\hat {\boldsymbol {p}}\boldsymbol {\cdot }\delta \boldsymbol {B}=0$ due to periodicity in $\lambda$. This gives

(3.2)\begin{equation} \hat{\boldsymbol{m}}\boldsymbol{\cdot}\delta \boldsymbol{B} ={-} \hat{\boldsymbol{m}}\boldsymbol{\cdot}\overline{\boldsymbol{B}} + \sqrt{B^{2} - (\hat{\boldsymbol{p}}\boldsymbol{\cdot}\overline{\boldsymbol{B}})^{2} - (\hat{\boldsymbol{n}}\boldsymbol{\cdot}\delta \boldsymbol{B})^{2}},\end{equation}

with $\hat {\boldsymbol {n}}\boldsymbol {\cdot }\delta \boldsymbol {B}$ an arbitrary function of $\lambda$. This must be combined with the condition $\overline {\hat {\boldsymbol {m}}\boldsymbol {\cdot }\delta \boldsymbol {B}}=0$, which determines $B^{2}$. The solution (3.2) illustrates two difficulties that also manifest in higher dimensions: first, $B^{2}$ must be computed self-consistently for a particular form of $\delta \boldsymbol {B}$; second, there may not exist real solutions for arbitrary choices of $\hat {\boldsymbol {n}}\boldsymbol {\cdot }\delta \boldsymbol {B}$ once $|\delta \boldsymbol {B}|$ approaches $|\overline {\boldsymbol {B}}|$. If one is not careful, either of these difficulties will almost inevitably cause discontinuities in $\delta \boldsymbol {B}$ as its solution is constructed. But, we reiterate that this does not signify that spherically polarised solutions do not exist, or are non-smooth; rather, they require constraining both components of $\delta \boldsymbol {B}$ simultaneously, as opposed to specifying one component and solving for the other. As a concrete example, Barnes & Hollweg (1974) show that for $\hat {\boldsymbol {n}}\boldsymbol {\cdot }\delta \boldsymbol {B} = \mathcal {A}_{n}\sin (k\lambda )$, solutions exist only for $\mathcal {A}_{n}<\mathcal {A}_{n,{\rm crit}}=({\rm \pi} /2)\sin \vartheta$, where $\vartheta = \cos ^{-1}(\hat {\boldsymbol {p}}\boldsymbol {\cdot }\overline {\boldsymbol {B}}/|\overline {\boldsymbol {B}}|)$. However, as we show below, if one starts with such a solution and evolves it according to (3.1a,b), there is no singular or otherwise interesting behaviour as $\mathcal {A}$ passes through $\mathcal {A}_{n,{\rm crit}}$; $\hat {\boldsymbol {n}}\boldsymbol {\cdot }\delta \boldsymbol {B}$ simply changes shape to avoid the amplitude limit. This example illustrates why ‘growing’ the $\delta \boldsymbol {B}$ solution from small amplitudes is necessary – without this, one is limited by the inability to choose an appropriate free function that will enable the constant-$B$ constraint to be satisfied.

With a small-amplitude $\delta \boldsymbol {B}$ specified, (3.1a,b) is easily solved by standard numerical methods. Here we use sixth-order finite differences for consistency with the 2-D and 3-D calculations, but a Fourier pseudospectral method is equally suitable in one dimension (M$+$21; Squire et al. 2022). An example is shown in figure 1 at several times as $\delta \boldsymbol {B}$ grows in amplitude. By the final snapshot, with $\mathcal {A}\approx 400$, the solution is approaching the zero-mean-field limit, which would give a purely stationary Alfvénic solution, since the propagation velocity $\boldsymbol {v}_{A}$ becomes much smaller than the flows $\delta \boldsymbol {u}$ in the nonlinear solution. Note that once $\delta \boldsymbol {B}\gg \overline {\boldsymbol {B}}$, $\phi$ becomes small and the shape of $\delta \boldsymbol {B}$ remains nearly constant with growing amplitude.

Figure 1. One-dimensional spherically polarised solutions from (3.1a,b), starting from a solution of (3.2) with $\mathcal {A}\approx 0.2$ constructed from a random collection of the first six Fourier modes. The wavevector $\hat {\boldsymbol {p}}$ is at an angle of $30^{\circ }$ to $\overline {\boldsymbol {B}}$, which lies in the $\hat {x}$ direction. Colours show $B_{x}$ (blue), $B_{y}$ (red), $B_{z}$ (yellow) and $B$ (black). From top to bottom, the panels show $\mathcal {A}\approx 0.2$ (the initial conditions), $\mathcal {A}\approx 0.65$, $\mathcal {A}\approx 5$ and $\mathcal {A}\approx 400$, illustrating how the solution approaches the zero-mean-field limit with $\delta \boldsymbol {B}\gg \overline {\boldsymbol {B}}$.

4. Two dimensions

To explore general 2-D solutions, we stipulate that $\delta \boldsymbol {B}$ is a function only of $\ell \equiv \hat {q}_{x} x + \hat {q}_{y}y$ and $z$, taking $\overline {\boldsymbol {B}} = {B_{0}}\hat {\boldsymbol {x}}$. This geometry is the obvious generalisation of the 1-D wave described above, with $\delta \boldsymbol {B}$ varying only in the plane angled at $\theta _{\rm 2D} \equiv \tan ^{-1}(\hat {q}_{y}/\hat {q}_{x})$ to the mean field. As in one dimension, the 2-D calculation proceeds in two steps by first constructing a low-amplitude near-linear solution, then growing this using (2.1) and (2.2). Both steps are significantly more complex than in one dimension.

Low-amplitude solution. First, we note that the clear generalisation of the $\hat {\boldsymbol {n}}$-directed 1-D linear Alfvénic field to two (or three) dimensions is the field $\delta \boldsymbol {B} = \boldsymbol {\nabla }\times (A_{x}\hat {\boldsymbol {x}})$, with $A_{x}$ chosen arbitrarily. The goal, then, is to add additional field components to enforce constant $B^{2}$, which is best done using the vector potential ($A_{y}$ and/or $A_{z}$) so as to maintain $\boldsymbol {\nabla }\boldsymbol {\cdot }\delta \boldsymbol {B}=0$ (the vector potential is not needed in the 1-D case because enforcing $\hat {\boldsymbol {p}}\boldsymbol {\cdot }\delta \boldsymbol {B}=0$ ensures $\boldsymbol {\nabla }\boldsymbol {\cdot }\delta \boldsymbol {B}=0$). The equation for $B^{2}$ becomes

(4.1)\begin{equation} B^{2} = {\rm const.} = |{B_{0}}-\partial_{z}A_{y} + \partial_{y} A_{z}|^{2} + |\partial_{z}A_{x} - \partial_{x}A_{z}|^{2} + |\partial_{x}A_{y} - \partial_{y}A_{x}|^{2},\end{equation}

where $\partial _{x}=\hat {q}_{x}\partial _{\ell }$ and $\partial _{y}=\hat {q}_{y}\partial _{\ell }$. If one of $A_{y}$ or $A_{z}$ is fixed, (4.1) is a first-order nonlinear partial differential equation in $A_{z}$ or $A_{y}$, which can in principle be solved using the method of characteristics. However, such a method, which is effectively that used by Valentini et al. (2019), is problematic on a periodic domain because the value of the constant $B^{2}$ is not known a priori, but itself depends on the solution ($A_{z}$ or $A_{y}$). An incorrect choice of $B^{2}$ manifests as an additional component of the mean field (as occurs for $\hat {\boldsymbol {m}}\boldsymbol {\cdot }\delta \boldsymbol {B}$ in (3.2)), which will require $\boldsymbol {A}$ to contain linear gradients, and thus lead to spurious discontinuities on a periodic domain. To overcome this, we instead stipulate that $\boldsymbol {\nabla }B^{2}=0$ in (4.1), which removes the issue of $B^{2}$ being undetermined at the cost of increasing the order of the partial differential equation. We also use the ‘mean-field Coulomb gauge’ $A_{y} = -\partial _{z}\alpha$, $A_{z}=\partial _{y}\alpha$ so that (4.1) involves just one free function without causing an artificial difference between the $\ell$ and $z$ directions. We then solve the resulting nonlinear partial differential equation for $\alpha (\ell,z)$ in Fourier space using MATLAB's fsolve function with the trust-region dogleg method. The two equations of $\boldsymbol {\nabla } B^{2}=0$ are easily combined into one by minimising only non-zero Fourier modes with fsolve, and the low-order approximate solution $(\partial _{y}^{2}+\partial _{z}^{2})\alpha = -(|\partial _{z}A_{x}|^{2} + |\partial _{z}A_{x}|^{2})/(2B_{0})$ provides a good initial guess for the optimisation. So long as $A_{x}$ is chosen to be relatively smooth (e.g. the first one or two Fourier modes in each direction), the procedure is rapid and straightforward because only a small number of Fourier modes are needed to represent $\alpha$. We have found no evidence for the development of discontinuities in such solutions, and a solution for $\alpha$ can be found for most choices of smooth, low-amplitude $A_{x}$ (at least in two dimensions; see below).

Large-amplitude solutions. Starting from a constant-$B$ initial condition, (2.1) and (2.2) can be solved using standard methods. For simplicity, we use an Euler timestepper with the timestep chosen based on the maximum of $\phi$ over the domain. The only complication arises from the non-homogenous derivative operator in the Poisson equation (2.2), which must be recomputed and solved at every timestep as $\delta \boldsymbol {B}$, and thus $\boldsymbol {\nabla }_{\perp }$, change. We choose to use a periodic sixth-order finite-difference representation,¹ forming $\boldsymbol {\nabla }_{\perp }^{2}$ as a sparse matrix through left-multiplication of the relevant gradient operators by $B_{i}$. A difficulty in the interpretation and solution of (2.2) is that the $\boldsymbol {\nabla }_{\perp }^{2}$ matrix has a rather high-dimensional nullspace, which means that, depending on the source $-\delta \boldsymbol {B}\boldsymbol {\cdot }\overline {\boldsymbol {B}}$, (2.2) may not have a solution.² We circumvent the issue by interpreting (2.2) in the least-squares sense, thus finding the best approximation to the flow $\tilde {\boldsymbol {u}}=\boldsymbol {\nabla }\phi$ that maintains constant $B$. We use the least-squares conjugate-gradient method (Paige & Saunders 1982) with an incomplete LU preconditioner. As we show below, the solution is generally very good, but the method certainly does not guarantee constant $B$, at least at finite resolution, unlike in one dimension. Presumably, if $-\delta \boldsymbol {B}\boldsymbol {\cdot }\overline {\boldsymbol {B}}$ has a large component in the nullspace of $\boldsymbol {\nabla }_{\perp }^{2}$, this implies that $\delta \boldsymbol {B}$ has developed structures that cannot continue to grow in amplitude while maintaining constant $B$. Understanding the conditions under which this occurs requires further study.

An example solution is shown in figure 2. We use a periodic domain of size $L_{\ell }=1$, $L_{z}=1$, with $384$ grid points in each direction and $\theta _{\rm 2D}=30^{\circ }$. The initial conditions in $A_{x}$ are constructed from a random combination of $|k_{\ell }|\leqslant 2{\rm \pi} /L_{\ell }$ and $|k_{z}|\leqslant 2{\rm \pi} /L_{z}$ modes, scaled to give $\mathcal {A}\approx 0.2$. In the top three panels, we show the field structure by taking an angled trace through the domain, accounting for the periodicity in order to show most of the solution (but note that some larger structures appear twice; the bottom-left panel illustrates the correspondence to the 2-D domain). The smooth, small-$\delta B_{x}$ initial conditions are shown in the top panel of figure 2 and have a very small variation in $B$, with $({\overline {B^{2}}})^{1/2}/\overline{B}\approx 10^{-5}$. As $\delta \boldsymbol {B}$ grows, it develops quite sharp gradients by modest amplitudes ($\mathcal {A}\gtrsim 0.7$; see second panel), which are numerically unresolved at this resolution. The large-amplitude $\mathcal {A}\approx 5$ solution is shown both in the bottom trace and in the image plots in the bottom row, illustrating the discontinuities that have developed at several locations in the domain. These sharp gradients represent a distinct difference compared to 1-D solutions, which only ever steepen modestly from the initial waveform as they grow to large amplitudes (see figure 1; this is discussed in more detail below). The method clearly maintains extremely constant $B$, aside from numerical ringing near some of the sharp gradient discontinuities that form (excluding these regions, $({\overline {B^{2}}})^{1/2}/\overline{B}\approx 0.8\,\%$).

Figure 2. Two-dimensional spherically polarised solution on an $\ell,z$ plane angled at $\theta _{\rm 2D}=30^{\circ }$ from the $\hat {\boldsymbol {x}}$ (mean-field) direction. The top three panels show periodic traces of $B$ (black), $B_{x}$ (blue), $B_{y}$ (red) and $B_{z}$ (yellow) along the white line plotted on the bottom-left panel (this lies at angle $\alpha \approx 11.3^{\circ }$ from $\hat {\boldsymbol {\ell }}$; integer $l$ values are marked to illustrate the correspondence between the trace and the 2-D solution). The initial condition is constructed from a random superposition of modes in $A_{x}$, scaled to give amplitude $\mathcal {A}\approx 0.2$ (top panel). It then grows in time according to (2.1) and (2.2). The bottom panels show the 2-D structure of the components of $\boldsymbol {B}$ at the time corresponding the bottom trace, when $\mathcal {A}\approx 5$. At least to the precision of the $384^{2}$ resolution used here, discontinuities develop in the field structure, unlike the 1-D solutions (the most prominent is near $l=1$ on the trace plots). Aside from numerical ringing caused by the development of these discontinuities, however, $B$ remains very constant throughout the domain (the colourbar of $B$ on the bottom left is scaled to ${\pm }2\,\%$).

We have also explored solutions on planes at other angles $\theta _{\rm 2D}$ (not shown) and with differing initial conditions. As we demonstrate in figure 3, which shows another example large-amplitude solution, the propensity to form sharp gradients depends on the structure of the initial low-amplitude $A_{x}$ (see below for further discussion). Unlike figure 2, the solution in figure 3 is well resolved and relatively smooth, constituting a practical demonstration that 2-D, smooth, spherically polarised solutions exist (which, as far as we are aware, was not previously known). We have not found any simple heuristic to determine how the sharpness of the final solution relates to the initial low-amplitude one. In general, solutions with larger $\theta _{\rm 2D}$ – i.e. those which are more elongated along the magnetic field – develop somewhat larger variation in $B$. This seems to be due to larger inaccuracies in the solutions of (2.2), although whether this relates to discontinuities remains unclear.

Figure 3. Same as figure 2 but starting from a different initial random collection of Fourier modes in the low-amplitude $A_{x}$ initial conditions. We show only the solution with $\mathcal {A}\approx 5$. In this case, discontinuous structures do not develop and the solution is well resolved at this resolution of $256^{2}$ (which is lower than that of figure 2).

5. Three dimensions

Three-dimensional solutions are constructed through a method that is almost identical to that of the 2-D case. The only additional complication is that the method to construct the low-amplitude solution using the mean-field Coulomb gauge ($A_{y} = -\partial _{z}\alpha$, $A_{z}=\partial _{y}\alpha$) cannot be used to construct a field that varies in $x$ but not in $y$ or $z$ (i.e. a field with power in $k_{\perp } = \sqrt {k_{y}^{2}+k_{z}^{2}}=0$ modes). Although a different method of constructing a low-amplitude solution may alleviate this, we opt instead to adjust the chosen $A_{x}$ so that its associated $B^{2} = |\partial _{y}A_{x}|^{2} + |\partial _{z}A_{x}|^{2}$ has very little power in $k_{\perp }=0$ modes. This is easily done using MATLAB's lsqnonlin function after constructing $A_{x}$ from a collection of random Fourier modes. The method seems to work well; for example, it generates smooth initial conditions with $\mathcal {A}\approx 0.2$ and $({\overline {B^{2}}})^{1/2}/\overline{B}\approx 3\times 10^{-5}$ (see figure 4). In solving (2.1) and (2.2), the only extra challenge compared with the 2-D case is the computational expense, and we are limited to constructing solutions with resolutions ${\lesssim }48^{3}$ due to our non-parallelised implementation using MATLAB. However, the only significant computational difficulty – inverting $\boldsymbol {\nabla }_{\perp }^{2}$ in (2.2) – involves standard iterative matrix-solve methods, which have robust and efficient parallel implementations that could allow for much higher resolutions if desired.

Figure 4. Three-dimensional spherically polarised solution in a cubic box with a resolution of $48^{3}$. As in figure 2, the top three panels show periodic traces of $B$ (black), $B_{x}$ (blue), $B_{y}$ (red) and $B_{z}$ (yellow) along a line in the direction $(\cos \theta _{\rm 3D}, \sin \theta _{\rm 3D}\cos \varphi _{\rm 3D},\sin \theta _{\rm 3D}\sin \varphi _{\rm 3D})$, with $\theta _{\rm 3D}\approx 30^{\circ }$ and $\varphi _{\rm 3D}\approx 11.3^{\circ }$, with $l=0$ at the centre of the box (the units are scaled to the size of the box). The initial condition is constructed from random modes in $A_{x}$ with an amplitude such that $\mathcal {A}\approx 0.2$ (top panel), then growing in time according to (2.1) and (2.2). The bottom panels show the 3-D structure of the components of $\boldsymbol {B}$ at the time corresponding the bottom trace, when $\mathcal {A}\approx 5$.

An example 3-D solution in a cubic box is shown in figure 4. As in figure 2, we show the solution structure using a trace along an angled, periodic line, at $\mathcal {A}\approx 0.2$ (the initial conditions), $\mathcal {A}\approx 0.7$ and $\mathcal {A}\approx 5$. Like in the 2-D case, the solutions appear to become discontinuous, although it is more difficult to diagnose in detail because of the limited resolution. The solutions have somewhat larger variation in $B^{2}$ compared with figure 2, with $({\overline {B^{2}}})^{1/2}/\overline{B}\approx 3\,\%$ by $\mathcal {A}\approx 5$, although this is at least partially a result of the lower resolution. We have also explored solutions in boxes that are elongated along $\overline {\boldsymbol {B}}$ with $L_{x}=4$, finding similar structures and properties, at least to the accuracy achievable here. We also note that we have confirmed that $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}$ remains zero within the tolerances of the finite-difference representation (with sixth-order finite differences at $48^{3}$ the root-mean-square of $\boldsymbol {\nabla }\boldsymbol {\cdot }\delta \boldsymbol {B}$ is ${\approx } 1.5\times 10^{-6}$).

6. Discussion and conclusions

Above, we have presented a method to construct large-amplitude spherically polarised Alfvénic structures, which form a broad class of nonlinear solutions to the compressible MHD (or kinetic MHD) equations. The method works by ‘growing’ a low-amplitude, nearly linear ($\delta \boldsymbol {B}\ll \overline {\boldsymbol {B}}$) wave, thus allowing the exploration of nonlinear solutions in the large-amplitude $\delta \boldsymbol {B}\sim \overline {\boldsymbol {B}}$ regime applicable to the solar wind, or even zero-mean-field limiting solutions with $\delta \boldsymbol {B}\gg \overline {\boldsymbol {B}}$. We have presented some examples of such solutions in one, two and three dimensions, albeit with somewhat limited resolution due to computational challenges. Others – including, in principle, solutions with turbulent-like spectra³ – could be generated using our method by choosing a different form for the initial, small-amplitude wave. These solutions exhibit some important features observed in the solar wind, particularly the asymmetry (one-sidedness) of $\delta B_{\|}=\delta \boldsymbol {B}\boldsymbol {\cdot }\overline {\boldsymbol {B}}/|\overline {\boldsymbol {B}}|=\delta B_{x}$ fluctuations at modest $\mathcal {A}$, which is a consequence of maintaining constant $B$ through large changes in $\delta \boldsymbol {B}$ (Gosling et al. 2009). At larger $\mathcal {A}$, all solutions cause magnetic field reversals ($|\delta B_{\|}|>|\overline {\boldsymbol {B}}|$, a commonly used definition of switchbacks), although these effects would be stronger if we chose more oblique solutions (e.g. larger $\theta _{\rm 2D}$ or a more elongated box in three dimensions; Mallet et al. 2021). However, it is also worth noting that despite its apparent success in nearly all cases we have explored, the method is not guaranteed to produce perfectly constant-$B$ structures because the Poisson-like equation (2.2), which determines how $\delta \boldsymbol {B}$ changes shape as it grows, lacks exact solutions in general, at least at finite resolution. These mathematical issues should be explored in more detail in future work. Other possible ideas for future studies include using a time-dependent mean field or gradient operators, which would change the large-amplitude $\delta \boldsymbol {B}$ that results from chosen small-amplitude initial conditions (M$+$21; Squire et al. 2022), or the extension to the relativistic regime (Mallet & Chandran 2021), which may have interesting applications to a number of high-energy processes such as disk coronae (Chandran, Foucart & Tchekhovskoy 2018) or pulsar magnetospheres (Kumar & Bošnjak 2020; Zhang 2020).

One interesting use case for these solutions is as input for nonlinear MHD or kinetic simulations to study processes such as parametric decay (instability of Alfvén waves; Sagdeev & Galeev 1969) or large-amplitude reflection-driven turbulence (Johnston et al. 2022). For example, a 3-D MHD simulation could be initialised with a 1-D, 2-D or 3-D solution $\delta \boldsymbol {B}$ from our method by first extending it via symmetry to three dimensions (e.g. $\delta \boldsymbol {B}_{\rm 3D}(x,y,z) = \delta \boldsymbol {B}(\hat {\boldsymbol {p}}\boldsymbol {\cdot } \boldsymbol {x})$ for a 1-D solution), then setting $\delta \boldsymbol {u} = \pm \delta \boldsymbol {B}/\sqrt {4{\rm \pi} \rho }$. Study of parametric decay especially has been somewhat limited by a lack of general, large-amplitude Alfvénic solutions other than circularly polarised waves that can be evolved numerically, so our method may enable important progress in this area (the modest-amplitude case with $\delta \boldsymbol {B}\lesssim \overline {\boldsymbol {B}}$ has been studied in one dimension by Del Zanna (2001); Del Zanna et al. (2015) and Tenerani et al. (2020) considered a larger-amplitude 2-D solution). A particularly interesting regime, which could in principle be studied for solutions of any dimensionality, would be the zero-mean-field limit, where the nonlinear Alfvénic solutions become a stationary, non-propagating tangle of magnetic field lines and Alfvénic flows. The system is reminiscent of the tangled fields possible in magnetostatic force-free equilibria (those with $\boldsymbol {u}=0$ and $\boldsymbol {B}\times (\boldsymbol {\nabla } \times \boldsymbol {B})=0$, see e.g. Chandrasekhar & Woltjer 1958; Marsh 1996; Hosking, Schekochihin & Balbus 2020), but with a different class of equilibria that includes flows and has seen comparatively little study.

A second reason for interest in these solutions relates directly to their properties and structure. As far as we are aware, it was previously not known whether smooth, spherically polarised solutions with 2-D or 3-D structure even existed; our method has provided near-perfect examples in two dimensions, one of which is shown in figure 3. In three dimensions, we are resolution-limited due to our non-parallelised numerical implementation (the example in figure 4 is not smooth at this resolution), but it seems unlikely that it would differ fundamentally from the 2-D case. Interestingly, in some other solutions, an example of which is shown in figure 2, $\boldsymbol {B}$ develops extremely sharp gradients that remain discontinuous at our highest resolution of $384^{2}$. We demonstrate this in figure 5, which quantifies the development of sharp structures by plotting how the maximum gradient of $\boldsymbol {B}$ evolves with solution amplitude, over a scan in resolution. For reference, we plot with dotted lines the empirically determined power-law behaviour of an unconverged solution ($||\boldsymbol {\nabla } \boldsymbol {B}||^{\infty }/||\boldsymbol {B}||^{\infty } \propto N^{\chi }$; the exponent $\chi$ differs between 1-D and 2-D cases). Convergence occurs at low resolution in 1-D solutions (figure 5a), which only ever develop modestly sharper structures compared with the initial wave, and are converged by $N\approx 48$ for this example ($b_{n}$ initialised with $k\leqslant 4{\rm \pi} /L$ modes). In two dimensions, both solutions become significantly sharper, but figure 5(b) (that from figure 2) shows no sign of convergence at all by $384^{2}$, while figure 5(c) (the case in figure 3) converges at around $128^{2}$. Given there does not seem to be any simple distinguishing feature(s) that determine whether discontinuities develop,⁴ the most likely scenario seems to be that continuous solutions do generally exist, but they require extremely sharp gradients in some cases. If true, this suggests 1-D waves are simply a particular special case of more general 2-D or 3-D fields that happens to allow large-amplitude solutions with an especially smooth structure.

Figure 5. Measurement of discontinuity formation in 1-D solutions (a) and 2-D solutions starting from two different initial conditions (b,c). We plot the normalised infinity norm of the gradient of the solutions, $||\boldsymbol {\nabla } \boldsymbol {B}||^{\infty }/||\boldsymbol {B}||^{\infty } = (1/3)\sum _{i}\sum _{j}\max (\boldsymbol {\nabla }_{j} B_{i})/\max (B_{i})$ (with $\boldsymbol {\nabla }_{j}$ taken along either $\lambda$ or $\ell$ and $z$), as a function of $\mathcal {A}$ for a scan in resolution in each case (an $N\times N$ grid is used in two dimensions; we list only every second $N$ in the legend for clarity). In one dimension (a), we initialise with a linear combination of modes with $k\leqslant 4{\rm \pi} /L$; in two dimensions, we initialise with a linear combination of modes with $k_{\ell }\leqslant 2{\rm \pi} /L_{\ell }$ and $k_{z}\leqslant 2{\rm \pi} / L_{z}$ (the $N=384$ case of b is that shown in figure 2). Dotted lines in each case show a scaling $N^{\chi }$, with $\chi$ chosen to match the scaling of unconverged solutions ($\chi \approx 0.7$ in one dimension and $\chi \approx 0.8$ in two dimensions). Clearly, the 1-D solution converges at very low resolution ($N\approx 64$), showing that (2.1) does not lead to particularly small-scale features. In contrast, in the 2-D case, sharp field structures form much more readily: in the first example in (b), which is that from figure 2, there is no convergence even at the highest resolution that is feasible using our current computational implementation ($N=384$); but, the second example in (c), which is that from figure 3 and simply starts from a different random initial condition, achieves convergences around $N=128$.

The conclusions of the previous paragraph and figure 5 may have interesting consequences for switchback formation in the solar wind. Specifically, our method based on (2.1) and (2.2) shares clear similarities with the physical processes that occur in the solar wind, where Alfvénic structures grow in normalised amplitude due to plasma expansion. Although the direct effect of expansion is not included in (2.1) and (2.2) (this would cause time dependence of the gradient operators), the general idea – whereby, in order to maintain constant $B^{2}$ as it grows, $\delta \boldsymbol {B}$ changes shape by means of a compressive flow – is very similar to physical expansion, and would presumably produce similar $\delta \boldsymbol {B}$ evolution.⁵ In support of this idea, in one dimension with expansion effects added, (2.1) and (2.2) become exactly the system of M$+$21, which was derived asymptotically from the expanding MHD equations and produces similar large-amplitude structures to (2.1) and (2.2). If this correspondence holds in the 3-D case, our results imply that small-amplitude Alfvénic perturbations released from the low corona would often develop very sharp gradients or discontinuities in the process of growing, unless they are in some sense 1-D (unlikely if they are created by turbulence in the chromosphere; van Ballegooijen et al. 2011). This is quite promising for the in situ Alfvénic scenario of switchback formation (Squire et al. 2020; Shoda et al. 2021), in which switchbacks are simply Alfvénic fluctuations that have grown to large amplitudes by expansion. While M$+$21 suggested, based on 1-D solutions, that the model may struggle to explain the extremely sharp switchback structures seen in some observations (Bale et al. 2019; Kasper et al. 2019; Akhavan-Tafti et al. 2021), the answer may simply be that near-discontinuities naturally develop when starting from low-amplitude fluctuations with more general (non-1-D) structure. Whether this idea has direct observable consequences – for example, in the relative populations of rotational versus tangential discontinuities (e.g. Neugebauer et al. 1984) – requires better understanding of why and how sharp gradients develop, so is left to future work.