1 Introduction↩︎

Since the introduction of Alfvén waves in 1942 [1], these fundamental magnetohydrodynamic (MHD) oscillations have been invoked to explain a wide range of phenomena in astrophysical [2], space [3]–[5], and laboratory plasmas [6]. Early solar wind observations showed that Alfvén waves are large-amplitude fluctuations characterized by quasi-constant \(|\boldsymbol{B}|\), \(p\), and \(\rho\) [5]. Recent observations from the Parker Solar Probe (PSP) [7], [8] reveal that the pristine solar wind in the upper corona [9] is permeated by large-amplitude Alfvénic fluctuations. Intriguingly, these fluctuations manifest as solitary, large-amplitude magnetic field reversals with quasi-constant \(|\boldsymbol{B}|\), accompanied by one-sided anti-sunward proton jets that exhibit near-perfect Alfvénic correlation [10]. They also display field-aligned electron strahls characterized by one-sided pitch angles both inside and outside the field reversal regions, suggesting topologically open magnetic field lines [11]. Consequently, some of these fluctuations are termed magnetic “switchbacks” (see also the recent review by [12]).

Due to the constancy of \(|\boldsymbol{B}|\), these Alfvénic fluctuations are termed spherically polarized Alfvén waves (SPAWs)—exact nonlinear solutions of the ideal MHD equations [13], [14]. Conventionally, the background magnetic field \(\boldsymbol{B}_0\) in SPAWs is defined as the ensemble average of \(\boldsymbol{B}\) (and hence by construction \(|\boldsymbol{B}_0| < |\boldsymbol{B}|\)). This definition, however, is problematic: (1) In in situ observations, the Alfvénic correlation \(\delta \boldsymbol{B}/\sqrt{\mu_0 \rho}=\pm \delta \boldsymbol{u}\) uniquely determines the perturbative part of \(\boldsymbol{B}\). Hence the local \(\boldsymbol{B}_0\) should be obtained as \(\boldsymbol{B}-\delta \boldsymbol{B}\). (2) The ensemble average of \(\boldsymbol{B}\) is inherently arbitrary, as it can vary substantially within a single solar wind stream, and consequently so does the Alfvén velocity \(\boldsymbol{V}_A = \boldsymbol{B}_0/\sqrt{\mu_0 \rho}\). [15] demonstrated that Alfvén waves in the solar wind are predominantly one-sided. This concept was subsequently extended by [16] to explain local radial velocity enhancements in Ulysses observations, where \(\boldsymbol{B}_0\) is defined on the constant sphere of \(|\boldsymbol{B}|\). The groundbreaking PSP observations in the upper corona have since reinforced this perspective: Alfvénic fluctuations in the most pristine solar wind appear as solitary perturbations on top of an otherwise unperturbed coronal magnetic field. By removing the perturbative component of \(\boldsymbol{B}\) constructed from the proton flow \(\boldsymbol{u}\) via the Alfvénic correlation, the background \(\boldsymbol{B}_0\) emerges naturally as the unperturbed coronal field. Consequently, Alfvén waves should be derivable from the ideal MHD equations as solitary solutions assuming only incompressibility (quasi-constant \(|\boldsymbol{B}|\), \(\rho\), and \(p\)) without a priori assumptions about \(\boldsymbol{B}_0\).

However, modeling SPAWs presents substantial challenges. Constructing a SPAW model involves three steps: (1) finding a magnetic field \(\boldsymbol{B}\) satisfying both \(\nabla \cdot \boldsymbol{B}=0\) and \(|\boldsymbol{B}|\simeq \mathrm{const}\); (2) decomposing \(\boldsymbol{B}\) into a background component \(\boldsymbol{B}_0\) and a fluctuating component \(\boldsymbol{B}_1\); and (3) constructing the velocity perturbation \(\boldsymbol{u}_1\) from \(\boldsymbol{B}_1\) via the Alfvénic correlation. Step (1) alone poses a highly non-trivial mathematical problem. [17] demonstrated the non-existence of 2D solutions (with \(\boldsymbol{B}\) restricted to a plane). Furthermore, assuming a constant far field (i.e., solitary solutions), 2.5D configurations (two spatial coordinates, three vector components) exist only when the vector field contains topologically closed regions (see Appendix 6). Therefore, maintaining open field-line topology requires \(\boldsymbol{B}\) to be genuinely 3D. Steps (2–3) are also non-trivial because \(\boldsymbol{u}_1\) depends on \(\boldsymbol{B}_1\), which in turn depends on the choice of \(\boldsymbol{B}_0\). This decomposition is well-defined only for solitary solutions; otherwise, the determination of \(\boldsymbol{B}_0\) remains ambiguous. Several studies have attempted to construct SPAWs. [18] and [19] developed analytic solitary 2.5D switchback models, which necessarily contain topologically closed regions. [20] and [21] constructed 3D turbulent switchbacks that lack spatial isolation, precluding a well-defined separation between \(\boldsymbol{B}_0\) and \(\boldsymbol{B}_1\).

In this study, we derive solitary Alfvén waves from the ideal MHD equations assuming only incompressibility (quasi-constant \(|\boldsymbol{B}|\), \(\rho\), and \(p\)), with \(\boldsymbol{B}_0\) naturally emerging as the constant far field. Based on this solution, we construct a solitary Alfvén wave model—an Alfvénon—characterized by quasi-constant \(|B|\) and open field-line topology. The remainder of this paper is organized as follows: Section 2 presents the solitary Alfvén wave solution. Section 3 constructs the Alfvénon. Section 4 presents results from MHD simulations of the Alfvénon. Finally, Section 5 discusses and concludes the study.

2 Solitary Alfvén Waves↩︎

We begin from the ideal MHD equations assuming adiabatic processes: \[\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot (\rho \boldsymbol{u}) &= 0,\tag{1}\\ \rho\left[\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \cdot \nabla) \boldsymbol{u} \right] &= -\nabla \left(p+\frac{B^2}{2\mu_0}\right) + \frac{1}{\mu_0} (\boldsymbol{B}\cdot \nabla) \boldsymbol{B},\tag{2}\\ \frac{\partial \boldsymbol{B}}{\partial t} &= \nabla \times (\boldsymbol{u} \times \boldsymbol{B}),\tag{3}\\ p \rho^{-\gamma} &= \mathrm{const},\tag{4}\\ \nabla \cdot \boldsymbol{B} &= 0,\tag{5} \end{align}\] where \(\rho\) is the plasma density, \(\boldsymbol{u}\) is the flow velocity, \(p\) is the pressure, \(\boldsymbol{B}\) is the magnetic field, \(\gamma\) is the adiabatic index, and \(\mu_0\) is the vacuum permeability. To proceed, we assume incompressibility: constant \(|\boldsymbol{B}|\), \(\rho\) and \(p\). Under these assumptions, we look for solitary Alfvénic solution of \(\boldsymbol{B}(\boldsymbol{r}, t)\) and \(\boldsymbol{u}(\boldsymbol{r}, t)\).

First, we make the conversion: \(\boldsymbol{b} = \boldsymbol{B}/\sqrt{\mu_0 \rho}\). Incompressibility guarantees \(\nabla\cdot \boldsymbol{u} = 0\) and \(\nabla \cdot \boldsymbol{b} = 0\). Consequently, Eqs. 2 –3 reduce to: \[\begin{align} \frac{\partial \boldsymbol{u}}{\partial t} &= (\boldsymbol{b}\cdot \nabla) \boldsymbol{b} - (\boldsymbol{u} \cdot \nabla) \boldsymbol{u},\tag{6}\\ \frac{\partial \boldsymbol{b}}{\partial t} &= (\boldsymbol{b} \cdot \nabla) \boldsymbol{u} - (\boldsymbol{u} \cdot \nabla) \boldsymbol{b}.\tag{7} \end{align}\] Under the constant-\(|\boldsymbol{b}|\) constraint, separating \(\boldsymbol{b}\) into a DC component \(\boldsymbol{b}_0\) and an AC perturbation \(\boldsymbol{b}_1\) is nontrivial. Previous derivations of shear, circularly polarized, and large-amplitude/spherically polarized Alfvén waves [1], [13], [14], [22] either leave the unperturbed state unspecified or set \(\boldsymbol{B}_0=\langle\boldsymbol{B}\rangle\); such an averaged background depends on the averaging window and generally does not satisfy \(|\boldsymbol{B}_0|=|\boldsymbol{B}|\), so it provides neither a deterministic AC/DC separation nor a unique Alfvén velocity. For a solitary wave packet, the separation is instead fixed by the far-field condition \[\boldsymbol{b}_1\to 0, \qquad \boldsymbol{b}_0 \equiv \lim_{r\to\infty}\boldsymbol{b}(\boldsymbol{r},t).\] Thus \(\boldsymbol{b}_0\) is the unique unperturbed (ground) state. Since \(|\boldsymbol{b}|\) is constant in the ideal formulation, \(|\boldsymbol{b}_0|=|\boldsymbol{b}|\), and \(\boldsymbol{b}_0\) uniquely defines the Alfvén velocity; \(\boldsymbol{b}_1=\boldsymbol{b}-\boldsymbol{b}_0\) is defined with respect to this state. Geometrically, the tip of \(\boldsymbol{b}\) remains on the constant-\(|\boldsymbol{b}|\) sphere, so the perturbation \(\boldsymbol{b}_1\) is constrained to the corresponding sphere defined by \(\boldsymbol{b}_0\), as illustrated in Fig. 1. Similarly, we decompose \(\boldsymbol{u}\) into \(\boldsymbol{u}_0 + \boldsymbol{u}_1\). Alfvénic solution dictates: \(\boldsymbol{u}_1 = \pm \boldsymbol{b}_1\). Because \(\boldsymbol{b}\) is Galilean invariant, without loss of generality, we can transform into a frame where \(\boldsymbol{u}_0 = 0\). In this frame, Eqs. 6 –7 reduce to: \[\begin{align} \frac{\partial \boldsymbol{u}_1}{\partial t} &= \boldsymbol{b}_0 \cdot \nabla \boldsymbol{b}_1, \tag{8}\\ \frac{\partial \boldsymbol{b}_1}{\partial t} &= \boldsymbol{b}_0 \cdot \nabla\boldsymbol{u}_1, \tag{9} \end{align}\] where \(\boldsymbol{b}_0 = \boldsymbol{B}_0/\sqrt{\mu_0 \rho}\) represents the Alfvén velocity. These yield the wave equations of solitary Alfvén waves: \[\begin{align} \frac{\partial^2 \boldsymbol{u}_1}{\partial t^2} &= (\boldsymbol{b}_0 \cdot \nabla)^2 \boldsymbol{u}_1, \tag{10}\\ \frac{\partial^2 \boldsymbol{b}_1}{\partial t^2} &= (\boldsymbol{b}_0 \cdot \nabla)^2 \boldsymbol{b}_1. \tag{11} \end{align}\]

This derivation should be read as the ideal constant-\(|\boldsymbol{B}|\) limit. In a finite domain, a finite-amplitude localized structure that reduces or reverses \(B_x\) cannot remain both strictly solitary and perfectly constant in \(|\boldsymbol{B}|\) while satisfying exact solenoidality and magnetic-flux conservation. The missing axial flux inside the perturbed region must be compensated by a slight compression or enhancement of neighboring field lines, so the numerical Alfvénon constructed below is a controlled quasi-constant-\(|\boldsymbol{B}|\) realization of the ideal solitary Alfvén-wave solution. The detailed flux-balance argument is given in Section 5.2.

An important consequence follows from Eqs. 8 –9 . Assuming \(\boldsymbol{b}_0 = b_0 \hat{x}\) and \(\boldsymbol{b}_1 = -\boldsymbol{u}_1\), Eq. 8 becomes \[\begin{align} \left(\frac{\partial}{\partial t}+b_0 \frac{\partial }{\partial x}\right)\boldsymbol{b}_1 = 0, \label{eq:propagation} \end{align}\tag{12}\] describing a forward-propagating (\(+\hat{x}\)) wave \(\boldsymbol{b}_1(x - b_0 t)\). Similarly, when \(\boldsymbol{b}_0 = -b_0 \hat{x}\), forward propagation requires \(\boldsymbol{b}_1 = \boldsymbol{u}_1\). Thus for forward-propagating waves, \(u_{1x}\) is always positive irrespective of the sign of \(\boldsymbol{b}_0\) (Fig. 1). This explains the one-sided anti-sunward proton jets associated with the SPAWs/switchbacks in the solar wind [10]–[12], [15], [16].

3 The Alfvénon Model↩︎

Despite the simplicity of the solitary Alfvén waves, constructing an Alfvénon is non-trivial. Without loss of generality, we seek a magnetic field \(\boldsymbol{B}(x,y,z)\) satisfying four constraints: (1) \(|\boldsymbol{B}|\simeq 1\), (2) \(\nabla\cdot \boldsymbol{B}=0\), (3) open field-line topology, and (4) constant far field (solitary). Constraints (1) and (4) are soft (approximate), while (2) and (3) are hard (exact).

3.1 The Iterative Algorithm↩︎

Given an arbitrary three-dimensional vector field \(\boldsymbol{F}(x,y,z)\), we apply the Helmholtz-Hodge decomposition: \[\begin{align} \boldsymbol{F} = \nabla \varphi + \nabla \times \boldsymbol{A}, \end{align}\] where \(\varphi\) is solved from the Poisson equation \(\nabla^2 \varphi = \nabla \cdot \boldsymbol{F}\). The divergence can be removed via: \[\begin{align} \boldsymbol{G} = \boldsymbol{F} - \nabla \varphi. \label{eq:remove95divergence} \end{align}\tag{13}\] The resulting field \(\boldsymbol{G}\) is then normalized to a unit vector field: \[\begin{align} \boldsymbol{F}' = \frac{\boldsymbol{G}}{|\boldsymbol{G}|}. \label{eq:normalization} \end{align}\tag{14}\] Such normalization reintroduces nonzero divergence. The iteration can therefore be viewed geometrically as an alternating-projection-type feasibility problem between two constraint sets, related to projection algorithms for convex feasibility and to alternating-projection methods for both convex and non-convex settings [23]–[25]. The Helmholtz-Hodge step is the \(L^2\)-orthogonal projection onto the solenoidal subspace within the chosen periodic Fourier representation and fixed mean field, whereas the normalization step is a pointwise projection onto the unit-\(|\boldsymbol{B}|\) manifold, i.e., a product of spheres wherever \(|\boldsymbol{G}|\) is nonzero. Because the unit-magnitude constraint set is non-convex, standard convergence theorems for alternating projections onto convex sets do not provide a global guarantee for the present algorithm. We therefore do not claim convergence for arbitrary initial fields; instead, convergence is treated as an empirical property of the smooth localized seeds considered here. The basin of attraction, possible failure modes, and rigorous convergence criteria remain open mathematical questions. However, with appropriate initial conditions \(\boldsymbol{F}\), alternating application of Eqs. 13 and 14 empirically converges. Denoting the fields at the \(n\)-th iteration by \(\boldsymbol{F}_n\) and \(\boldsymbol{G}_n\), we have:

\[\begin{align} \boldsymbol{G}_{n} &= \boldsymbol{F}_{n} - \nabla \varphi_n, \tag{15}\\ \boldsymbol{F}_{n+1} &= \frac{\boldsymbol{G}_n}{|\boldsymbol{G}_n|} \tag{16}, \end{align}\] where \(\nabla^2 \varphi_n = \nabla \cdot \boldsymbol{F}_n\). Equations 15 and 16 together constitute the iterative algorithm, hereafter referred to simply as the algorithm. In practice, the Helmholtz-Hodge decomposition is implemented in Fourier space (and thus assuming periodic boundary conditions) by solving the Poisson equation: \[\begin{align} \widehat{\nabla\cdot \boldsymbol{F}}=i \boldsymbol{k} \cdot \widehat{\boldsymbol{F}}=\widehat{\Delta\varphi}=-k^2\widehat{\varphi}. \end{align}\] and thus the potential field is obtained to be: \[\begin{align} \widehat{\varphi} = -\frac{i \boldsymbol{k} \cdot \widehat{\boldsymbol{F}}}{k^2}, \end{align}\] where the zero-\(k\) singular point \(\widehat{\varphi}(0)\) is set to be an arbitrary value. The solenoidal field can then be obtained: \[\begin{align} \begin{aligned} \widehat{\boldsymbol{F_\perp}}=\widehat{\boldsymbol{F}}-\widehat{\nabla \varphi}&=\widehat{\boldsymbol{F}}-\left[i\boldsymbol{k} \left(-\frac{i \boldsymbol{k} \cdot \widehat{\boldsymbol{F}}}{k^2}\right)\right]\\ &= \widehat{\boldsymbol{F}}-\boldsymbol{k} \left(\frac{\boldsymbol{k} \cdot \widehat{\boldsymbol{F}}}{k^2}\right). \end{aligned} \end{align}\] \(\widehat{\boldsymbol{F_\perp}}\) can then be inverse-transformed into real space, which is denoted as \(\boldsymbol{G}\). The values of \(\widehat{\boldsymbol{F}}\) and \(\widehat{\boldsymbol{F_\perp}}\) on the Nyquist planes are set to zero to avoid aliasing. Otherwise, significant divergence persists in \(\boldsymbol{G}\) for \(|k| \geq k_{\text{Nyquist}}\).

3.2 Magnetic Field↩︎

We construct the Alfvénon model on a \(128^3\) grid, with all spatial coordinates \(x\), \(y\) and \(z\) ranging from 0 to 1. We start from an initial field: \[\begin{align} \boldsymbol{F}_0 = \boldsymbol{B}_0 + A \cdot \left[ \cos(\phi) \hat{y} + \sin(\phi) \hat{z} \right] \cdot \exp\left[- \frac{\Delta r^2}{2\sigma^2} \right] \end{align}\] where \(\boldsymbol{B}_0 = (1,0,0)\), \(A=10.0\), \(\phi = 2\pi k_x x\), \(k_x = 4\), \(\Delta r = \sqrt{(x-0.5)^2+(y-0.5)^2+(z-0.5)^2}\), and \(\sigma = 1/30\). Because of the Gaussian envelope, \(\boldsymbol{F}_0\) is neither divergence-free nor of constant magnitude, making it a suitable input for the algorithm. Convergence is tracked with the vector field difference:

\[\|\Delta \boldsymbol{G}_n\|_2 = \sqrt{ \sum_{i,j,k} \big|\boldsymbol{G}_n - \boldsymbol{G}_{n-1}\big|^2 \,\Delta x\,\Delta y\,\Delta z },\] where \(\Delta x = \Delta y = \Delta z = 1/128\). The results are shown in Figure 2 (a). \(\|\Delta \boldsymbol{G}_n\|_2\) drops rapidly within the first 25 iterations and then asymptotically converges. Similarly, the standard deviation \(\sigma_{|\boldsymbol{B}|}\) of \(\boldsymbol{G}_n\) decreases rapidly and becomes negligible as the iterations proceed. Due to this convergence, the iteration is stopped at \(n=200\), and \(\boldsymbol{G}_{200}\) is taken as a candidate magnetic field for the Alfvénon model, hereafter denoted as \(\boldsymbol{G}_A\). We choose \(\boldsymbol{G}_A\) as the candidate magnetic field because it is strictly solenoidal, satisfying the hard constraint \(\nabla\cdot\boldsymbol{B}=0\) up to machine precision (\(10^{-10}\)), while retaining a quasi-constant magnitude that satisfies the soft constraint \(|\boldsymbol{B}|\simeq 1\) to high accuracy.

The resulting field \(\boldsymbol{G}_A\) exhibits a spatially localized, nontrivial twisting of otherwise unperturbed open magnetic field lines oriented along \(+\hat{x}\). Outside the perturbed region, the field is approximately uniform, \(\boldsymbol{B} \simeq \boldsymbol{B}_0 = (1,0,0)\). The distribution of \(|\boldsymbol{B}|\) in \(\boldsymbol{G}_A\) is shown in Fig. 2 (b), demonstrating that \(|\boldsymbol{B}|\) remains nearly constant on the unit sphere. To characterize the structure of the solution, we define the local deflection angle from \(\boldsymbol{B}_0\) as \[\theta(x,y,z) = \cos^{-1}\!\left(\frac{B_x}{|\boldsymbol{B}|}\right).\] Contours of \(\theta\) are displayed in Fig. 3, with its distribution shown in Fig. 2 (c). The maximum deflection angle, \(\theta_{\max} \simeq 45.0^\circ\), occurs at the pair of reflection-symmetric grid points \(P_1(ix = 59,\, iy = 70,\, iz = 64)\) and \(P_2(ix = 69,\, iy = 58,\, iz = 64)\), where \(ix\), \(iy\), and \(iz\) denote the integer indices on the \(128^3\) grid (ranging from 0 to 127). The red field line shown in Fig. 3 passes through \(P_1\), with neighboring field lines shown in gray. All field lines enter the computational domain through the \(x = 0\) plane and exit through the \(x = 1\) plane, confirming an open field-line topology throughout (see Appendix 7 for details of the field-line tracer). Outside the perturbed region, field lines remain essentially unperturbed. Field lines traversing the perturbed region undergo localized twisting while preserving field-line density (\(|\boldsymbol{B}| \simeq 1\)), and subsequently relax back to the unperturbed state \(\boldsymbol{B}_0\) downstream. The \(\theta_{\max} \simeq 45^\circ\) configuration studied here should be regarded as a moderate-amplitude prototype; the near-reversal, large-amplitude branch of solitary Alfvén waves will be treated separately in future work.

3.3 The Alfvénon Model↩︎

The quasi-constant far field \(\boldsymbol{B}_0\) enables a clear separation of the perturbation \(\boldsymbol{B}_1\) from \(\boldsymbol{B}\). Based on the Alfvénic correlation, we construct \(\boldsymbol{u}_1 = -\boldsymbol{b}_1 = -\boldsymbol{B}_1\) to ensure forward propagation along \(\boldsymbol{B}_0\), where we have adopted normalized units with \(\rho = 1\) and \(\mu_0 = 1\), and set \(\boldsymbol{u}_0 = 0\). These \(\boldsymbol{B}\) and \(\boldsymbol{u}\) fields serve as initial conditions for the MHD simulations.

4 MHD Simulations↩︎

To validate that the Alfvénon model constitutes a numerical representation of solitary Alfvén waves, we perform MHD simulations using the LAPS code [26], a pseudo-spectral solver for the ideal MHD equations with periodic boundary conditions. The plasma pressure is set to \(p = 0.05\), corresponding to a plasma beta \(\beta = 2p/B^2 = 0.1\) consistent with pristine solar wind conditions [27]. A polytropic index \(\gamma = 1.2\) is adopted to approximate realistic fast solar wind conditions [28]. Both viscosity and resistivity are set to zero; dissipation arises solely from numerical dealiasing. To minimize interactions between periodic images, the Alfvénon model (defined on a \(128^3\) grid) is embedded within a larger domain by appending uniform \(\boldsymbol{B}_0\) regions of size \(128^3\) on either side along the \(x\)-direction. This produces an initial grid of \(640\times128\times128\) spanning a domain of size \(5\times 1\times 1\) (\(L_x \times L_y \times L_z\)), with the Alfvénon centered in the middle third.

Simulation results are presented in Figure 4 at times \(t = 0.0\), \(5.0\), \(10.0\), and \(100.0\), where \(t\) denotes dimensionless simulation time. Column (a) shows the one-dimensional wave-packet profile along the \(x\)-direction at \(iy = 64\) and \(iz = 64\). Profiles are shifted to a common reference position assuming a propagation speed of 1, with only the central 128 grid points displayed. All quantities are nondimensionalized such that \(t = 1\) corresponds to the Alfvénon traversing one unit-length box, while crossing the full domain of length \(L_x = 5\) requires \(t = 5\). The near-perfect alignment of profiles demonstrates that the Alfvénon propagates at \(V_A =|\boldsymbol{B}_0|/\sqrt{\rho} \simeq 1\) while maintaining spatial coherence with only mild relaxation. The perturbations preserve Alfvénic correlations in all three components, including \(B_x\) and \(u_x\) (along \(\boldsymbol{B}_0\)). Column (b) displays the spatial distribution of \(\theta\), revealing negligible nonlinear evolution up to \(t = 10.0\) and gradual relaxation by \(t = 100.0\).

The deformation of the Alfvénon arises primarily from phase mixing. Column (c) displays the local \(V_A\). Inside the Alfvénon, \(V_A\) varies by approximately 1%, producing differential phase speeds that gradually deform the wave packet. Furthermore, for slabs perpendicular to the \(x\)-axis, \(V_A\) near the Alfvénon slightly exceeds that of the unperturbed region (by \(\lesssim 0.5\%\)), accounting for the gradual forward drift visible in column (a). This \(V_A\) distribution becomes progressively disrupted as the simulation proceeds, growing increasingly irregular by \(t=100.0\) due to nonlinear evolution.

Nonlinear evolution also generates density fluctuations, violating perfect Alfvénicity. Column (d) shows the density profiles. By construction, the initial Alfvénon contains no density perturbations. However, the model inevitably includes minor defects, such as high-frequency modes near the Nyquist frequency (suppressed by dealiasing procedures) and magnetic pressure imbalances arising from the slight non-constancy of \(|\boldsymbol{B}|\). These defects induce density perturbations as the simulation progresses, which in turn drive nonlinear evolution. Notably, the spatial pattern of \(\rho\) coincides with that of \(V_A\), indicating that density variations dominate over \(|\boldsymbol{B}|\) variations in causing the phase mixing.

To further investigate the sources of the observed relaxation, we performed two additional simulations beyond the initial run (Run 1: \(L_x=5\), \(\alpha=0.495\), where \(\alpha\) is the controlling parameter of dealiasing option 2 in LAPS; dealiasing becomes less effective as \(\alpha \to 0.5\); see Appendix 8 for details). These additional runs either suppress the dealiasing effect (Run 2: \(L_x=5\), \(\alpha=0.499\)) or enhance the effects of periodic boundary conditions (Run 3: \(L_x=1\), \(\alpha=0.495\)). Figure 5 presents the temporal evolution of energy diagnostics across all three runs: magnetic fluctuation energy \(E_B = \frac{1}{2}|\boldsymbol{B}-\boldsymbol{B}_0|^2\); kinetic energy \(E_K = \tfrac{1}{2} \rho |\boldsymbol{u}|^2\); incremental internal energy \(\Delta E_{\mathrm{int}} = (p - p_0)/(\gamma - 1)\), where \(p_0\) denotes the initial pressure; total fluctuation energy \(E_B + E_K\); residual energy \(E_r = E_K - E_B\), which vanishes for perfectly Alfvénic fluctuations; and total energy \(E_B + E_K + \Delta E_{\mathrm{int}}\). All diagnostics are integrated over the simulation domain. For comparison, the total fluctuation energy from Runs 2 and 3 is also displayed in panel (a).

Across all three runs, \(E_B\) and \(E_K\) gradually decrease while \(\Delta E_{\mathrm{int}}\) increases, with total energy \(E_B + E_K + \Delta E_{\mathrm{int}}\) conserved to high precision (\(\sim 10^{-9}\)). This confirms that fluctuation energy is converted entirely to internal energy via compressive work, with negligible numerical dissipation. Comparison of the runs reveals distinct influences from dealiasing and boundary effects. In Run 2, reduced dealiasing lowers the heating rate (dash-dotted line in Fig. 5 (a)), indicating that a substantial fraction of the compression is numerical. Nevertheless, perturbations in both Runs 1 and 2 remain perfectly Alfvénic (purple lines in panels a and b), with small density fluctuations (panel d). In Run 3, the shorter domain (\(L_x = 1\)) enhances nonlinear interactions from periodic boundaries (dashed line in panel a). While Runs 1 and 3 evolve nearly identically for the first 10 time units, Run 3 subsequently develops significant density fluctuations (panel d), generating non-zero residual energy (purple line in panel c) and accelerating heating (green line in panel c).

This domain-size dependence also clarifies the role of parametric decay instability (PDI) in the present simulations. PDI is naturally associated with extended, periodic, large-amplitude Alfvénic wave trains, whereas the object considered here is a localized solitary packet embedded in an otherwise unperturbed far field. In a periodic numerical domain, however, the packet is repeated by the boundary conditions; if the box is too short, these periodic images effectively increase the filling factor of the perturbation and make the system more wave-train-like, allowing PDI-like compressive sidebands to grow more efficiently. The comparison between the longer-domain runs and Run 3 shows that, in the finite-time periodic simulations presented here, the PDI-like signatures are strongly enhanced by finite-domain and periodic-image effects and are suppressed when the Alfvénon is more isolated from its periodic images. We therefore do not identify PDI as the dominant process controlling the simulated Alfvénon evolution over the present time interval. A full stability theory of isolated or open-boundary solitary Alfvén waves with respect to PDI is beyond the scope of this paper and will be addressed in future work.

Overall, these results demonstrate that the Alfvénon behaves as a nonlinear solitary Alfvénic solution under ideal MHD evolution.

5 Discussion and Conclusions↩︎

5.1 Properties of Solitary Alfvén Waves↩︎

In Section 2, we derived solitary Alfvén waves from the ideal MHD equations under the assumption of incompressibility. Although Eqs. 10 –11 formally resemble those governing classical shear [1], circularly polarized [22], or spherically polarized [13] Alfvén waves, two critical distinctions arise:

• \(\boldsymbol{b}_0\) is the unperturbed field, uniquely determined by the constant far field.

• \(\boldsymbol{b}_1\) is solitary, representing a localized perturbation of an otherwise unperturbed background field.

These distinctions lead to three non-trivial consequences: First, \(\boldsymbol{b}_1\) is intrinsically nonlinear, since if both \(\boldsymbol{b}_1\) and \(\boldsymbol{b}_1'\) individually satisfy the constant-\(|\boldsymbol{b}|\) constraint, their superposition \(\boldsymbol{b}_1 + \boldsymbol{b}_1'\) generally violates that constraint, so linear superposition fails and challenges the conventional Fourier decomposition employed in Alfvénic turbulence studies. Second, the Alfvén speed is uniquely defined, because \(\boldsymbol{b}_0\) is fixed by the constant far field and, expressed in Alfvén units, directly represents the unperturbed Alfvén velocity, eliminating the ambiguity inherent in classical ensemble-averaged definitions of the background field. Third, field-line twisting is non-trivial, because Alfvén’s theorem dictates that magnetic field topology is preserved in ideal MHD and simultaneously maintaining \(\delta|\boldsymbol{B}|\ll|\delta\boldsymbol{B}|\) with a constant unperturbed far field \(\boldsymbol{B}_0\) necessitates non-trivial twisting to preserve quasi-constant field-line density, i.e. \(|\boldsymbol{B}|\).

5.2 Subtleties of the Alfvénon Model↩︎

The construction of the Alfvénon model in Section 3 involves operations in Fourier space, and thus inherits limitations associated with spectral methods, most notably the Gibbs phenomenon. The algorithm naturally generates sharp gradients at the edges of the perturbed region, which cannot be accurately represented due to spectral ringing. These overshoots violate the constant-\(|\boldsymbol{B}|\) constraint and are subsequently suppressed by the iterative procedure, thereby shifting spectral energy toward high-\(k\) modes. This effect also introduces complications in MHD simulations: the LAPS code advances in Fourier space and therefore cannot properly resolve strong discontinuities, while the dealiasing procedure suppresses high-\(k\) components, inevitably affecting the evolution of the Alfvénon. For the parameters adopted in this study, such effects remain modest; however, for larger-amplitude Alfvénons, the accumulation of spectral energy at high wavenumbers becomes more pronounced, potentially leading to numerical instabilities.

As noted already in Section 2, a more fundamental limitation concerns the constant-\(|\boldsymbol{B}|\) constraint itself. We enforce constant \(|\boldsymbol{B}|\) only as a soft constraint because a strictly localized solution with exactly constant \(|\boldsymbol{B}|\) is incompatible with magnetic-flux conservation in a finite domain. As argued in [19] [see their Eq. (10)], a cylindrical tube aligned with the \(x\)-axis and enclosing the perturbation must have the same total magnetic flux through every cross-section; however, any finite-amplitude reduction or reversal of \(B_x\) inside the perturbed region removes axial flux that must be compensated elsewhere. The compensation appears as a slight compression or enhancement of neighboring field lines, and hence as a small departure from perfectly constant \(|\boldsymbol{B}|\), unless strict localization is relaxed. The Alfvénon model is therefore most accurately regarded as a controlled quasi-constant-\(|\boldsymbol{B}|\) and quasi-solitary realization of the ideal solitary Alfvén-wave formulation, consistent with the two soft constraints introduced in Section 3.

Figure 4 (c1) is a manifestation of this effect in the \(x\)-\(y\) plane. Since \(V_A = |\boldsymbol{B}|/\sqrt{\rho}\) and \(\rho=1\) in the initial condition, the panel effectively visualizes \(|\boldsymbol{B}|\), which is slightly compressed (by \(\sim 0.3\%\)) in the slab containing the solution. Figure 6 illustrates this effect in greater detail through 2D slices of \(B_x\) and \(|\boldsymbol{B}|\) in both the unperturbed (\(i_x=0\)) and perturbed (\(i_x=64\)) regions. In the unperturbed region (panel a), \(\boldsymbol{B}\) closely approximates \(\boldsymbol{B}_0\) with only minor deviations (\(\sim\) 0.1%). In the perturbed region (panel b), the local reduction in \(B_x\) is compensated by a slight enhancement of \(B_x\) in the surrounding area, ensuring that the total magnetic flux is strictly identical across the two panels—a direct consequence of the strictly enforced solenoidal condition \(\nabla\cdot\boldsymbol{B} = 0\). Panels (c) and (d) show the corresponding \(|\boldsymbol{B}|\) distributions, revealing that the perturbed region slightly compresses the neighboring field lines and thereby increases the local field-line density, i.e.\(|\boldsymbol{B}|\), by \(\sim\) 0.3 %. This compression is fully consistent with Fig. 4 (c1).

This effect is expected to diminish if the solution is constructed in a larger domain while preserving the size of the perturbed region. Solitary Alfvén wave packets with \(B_x\) reversals are therefore inherently space-filling, which may have interesting consequences for the solar wind. Within the bounded coronal hole outflow, the presence of switchbacks should compress the surrounding field lines while preserving the total magnetic flux through each cross-section. Observations have consistently shown that \(|\boldsymbol{B}|\) decreases more slowly than \(R^{-2}\) in near-Sun coronal hole outflows [10], [29], deviating from the spherical expansion trend expected from remote sensing. The space-filling nature of switchbacks may therefore directly account for this deviation. Future work will investigate this in detail.

5.3 Conclusions↩︎

To the best of our knowledge, the Alfvénon represents the first numerical realization of a solitary Alfvén wave packet. The Alfvénon exhibits nontrivial three-dimensional twisting of magnetic field lines while preserving quasi-constant \(|\boldsymbol{B}|\). Its complex structure and coherent finite-time propagation suggest that fundamental aspects of Alfvén wave physics remain incompletely understood, more than eight decades after their original conceptualization. Moreover, the ubiquity of solitary Alfvén waves in the highly magnetized solar corona suggests that such localized structures may be an important form of Alfvénic fluctuations in astrophysical environments.

Despite its simplicity, the iterative algorithm yields nontrivial results. Two natural follow-up studies arise: (1) investigating the dependence of the Alfvénon structure on the initial amplitude \(A\). The nontrivial field-line twisting revealed here suggests novel behaviors of magnetic fields in highly magnetized plasmas where \(|\boldsymbol{B}|\) is constrained to remain constant. Our iterative algorithm thus provides a foundation for future investigations of solitary Alfvén wave physics across a variety of contexts. (2) Conducting direct MHD simulations of Alfvénon collisions. Nearly all MHD turbulence phenomenologies depend on the interaction of counter-propagating Alfvén wave packets [30]–[37]. Our model, for the first time, enables direct simulation of such collisions in isolated setups. Both topics will be addressed in forthcoming companion papers.

Finally, we note that spherically polarized Alfvén waves—and hence Alfvénons—are also exact solutions of relativistic MHD [38]. In strongly magnetized environments where the Alfvén speed approaches the speed of light, relativistic Alfvénons could transport energy far more efficiently than classical shear Alfvén waves via their ultra-relativistic jets. Recent studies have proposed Alfvén waves as viable drivers of fast radio bursts (FRBs) from magnetars [39]–[43]. Exploring the behavior, stability, energy transport, and collision properties of relativistic Alfvénons represents an exciting direction for future research.

6 3D Dependence of Simply Connected Magnetic Solitary Solution↩︎

Let \(\mathbf{V}:\mathbb{R}^3 \to \mathbb{R}^3\) be a vector field that approaches a constant value at spatial infinity: \[\mathbf{V}(\mathbf{x}) \to \mathbf{V}_\infty \quad\text{as }\|\mathbf{x}\|\to\infty.\] No assumptions are imposed on the divergence, curl, or magnitude of \(\mathbf{V}\). Suppose that in some coordinate system \((q_1,q_2,q_3)\) the field depends on only two variables: \[\mathbf{V}(q_1,q_2,q_3)=\mathbf{V}(q_1,q_2).\]

Claim↩︎

A nontrivial isolated configuration of this type is only possible if the suppressed coordinate \(q_3\) parametrizes a compact loop in physical space. If \(q_3\) is unbounded, the only vector field satisfying the asymptotic condition is the trivial uniform field \(\mathbf{V}_\infty\).

Argument↩︎

For each fixed \((q_1,q_2)\), the field remains constant as \(q_3\) varies. Define the set \[\Gamma_{(q_1,q_2)} = \{ (q_1,q_2,q_3) : q_3 \in I \},\] where \(I\) is the range of \(q_3\). All points in \(\Gamma_{(q_1,q_2)}\) share the same field value \(\mathbf{V}(q_1,q_2)\).

Case 1: \(q_3\) unbounded. If \(q_3\) ranges over an unbounded interval, then \(\Gamma_{(q_1,q_2)}\) contains points with arbitrarily large Euclidean norm. Along this direction, \[\mathbf{V}(q_1,q_2,q_3) = \mathbf{V}(q_1,q_2) \quad\text{for all } q_3.\] If \(\mathbf{V}(q_1,q_2) \neq \mathbf{V}_\infty\), then the field fails to approach \(\mathbf{V}_\infty\) along these unbounded lines, contradicting the assumed asymptotic behavior. Thus the condition at infinity forces \[\mathbf{V}(q_1,q_2)=\mathbf{V}_\infty \quad\forall (q_1,q_2),\] so \(\mathbf{V}\) is constant everywhere.

Case 2: \(q_3\) compact. If \(q_3\) parametrizes a topologically closed region (e.g.an angular coordinate), then each set \(\Gamma_{(q_1,q_2)}\) is bounded. Nontrivial dependence on \((q_1,q_2)\) does not produce directions along which the field propagates to infinity while remaining fixed. Therefore the field may coincide with \(\mathbf{V}_\infty\) outside a sufficiently large ball while retaining nontrivial structure in a compact region.

Conclusion↩︎

A vector field on \(\mathbb{R}^3\) that becomes uniform at infinity and depends on only two coordinates can exhibit nontrivial localized structure only if the remaining coordinate is spatially bounded, forming topologically closed regions. If the unused coordinate is unbounded, the asymptotic condition forces the field to be identically equal to the far-field constant \(\mathbf{V}_\infty\).

7 Field Line Tracer↩︎

A Runge-Kutta method is employed to trace the field lines. Given the magnetic grid \(\{B_x,B_y,B_z\}\) on a periodic box of size \(\mathbf{L}\), each component is trilinearly interpolated to an arbitrary point \(\mathbf{x}\): \[\mathbf{B}(\mathbf{x}) = \bigl(B_x(\mathbf{x}),\,B_y(\mathbf{x}),\,B_z(\mathbf{x})\bigr),\] with indices wrapped periodically before interpolation. The local direction field is then \[\hat{\mathbf{b}}(\mathbf{x}) = \begin{cases} \dfrac{\mathbf{B}(\mathbf{x})}{\lVert \mathbf{B}(\mathbf{x}) \rVert}, & \lVert \mathbf{B}(\mathbf{x}) \rVert > 10^{-12},\\[0.75em] \mathbf{0}, & \text{otherwise}. \end{cases}\]

Starting from a seed \(\mathbf{x}_0\), Runge–Kutta 4 (step size \(h\)) advances along this interpolated, unit-magnitude field: \[\begin{align} \mathbf{k}_1 &= \hat{\mathbf{b}}(\mathbf{x}_n),\\ \mathbf{k}_2 &= \hat{\mathbf{b}}\!\left(\mathbf{x}_n + \tfrac{1}{2} h \mathbf{k}_1\right),\\ \mathbf{k}_3 &= \hat{\mathbf{b}}\!\left(\mathbf{x}_n + \tfrac{1}{2} h \mathbf{k}_2\right),\\ \mathbf{k}_4 &= \hat{\mathbf{b}}(\mathbf{x}_n + h \mathbf{k}_3),\\ \mathbf{x}_{n+1} &= \operatorname{mod}\!\left(\mathbf{x}_n + \tfrac{h}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4),\, \mathbf{L}\right). \end{align}\]

Each step samples the trilinearly interpolated \(\hat{\mathbf{b}}\) at the intermediate RK4 positions and wraps the result back into the periodic domain, producing the traced field line forward and/or backward from the seed.

8 Dealiasing in LAPS↩︎

LAPS provides two dealiasing options, summarized below.

Option 1 (circular-padding / \(2/3\) rule) applies a sharp spectral cutoff, retaining only modes with \(k \le \tfrac{2}{3}k_{\max}\): \[G_1(k)= \begin{cases} 1, & k \le \tfrac{2}{3}k_{\max},\\ 0, & k > \tfrac{2}{3}k_{\max}. \end{cases}\]

Option 2 (smoothing filter) applies a smooth rational filter with \(\theta=\pi k/k_{\max}\): \[G_2(\theta;\alpha)=\frac{a_j+b_j\cos\theta+c_j\cos(2\theta)}{1+2\alpha\cos\theta}, \quad a_j=\frac{5+6\alpha}{8},\; b_j=\frac{1+2\alpha}{2},\; c_j=-\frac{1-2\alpha}{8}.\] As \(\alpha \to 0.5\), the filter becomes less dissipative (closer to unity), so high-\(k\) suppression weakens.

References↩︎