Stream-Aligned Magnetohydrodynamics for Solar Wind Simulations

We present a reduced magnetohydrodynamic (MHD) mathematical model describing the dynamical behavior of highly conducting plasmas with frozen-in magnetic fields, constrained by the assumption that, there exists a frame of reference, where the magnetic field vector, $\mathbf{B}$, is aligned with the plasma velocity vector, $\mathbf{u}$, at each point. We call this solution"stream-aligned MHD"(SA-MHD). Within the framework of this model, the electric field, $\mathbf{E} = -\mathbf{u} \times \mathbf{B} \equiv 0$, in the induction equation vanishes identically and so does the electromagnetic energy flux (Poynting flux), $\mathbf{E}\times\mathbf{B}\equiv0$, in the energy equation. At the same time, the force effect from the magnetic field on the plasma motion (the Ampere force) is fully taken into account in the momentum equation. Any steady-state solution of the proposed model is a legitimate solution of the full MHD system of equations. However, the converse statement is not true: in an arbitrary steady-state magnetic field the electric field does not have to vanish identically (its curl has to, though). Specifically, realistic tree-dimensional solutions for the steady-state ("ambient") solar atmosphere in the form of so-called Parker spirals, can be efficiently generated within the stream-aligned MHD (SA-MHD) with no loss in generality


INTRODUCTION
Space weather describes the dynamic state of the Earth's magnetosphere-ionosphere system, which is driven by the solar wind and solar ionizing radiation. The greatest disturbances in space weather are geomagnetic storms, the most severe of which are caused by coronal mass ejections (CMEs) (see (Gosling 1993)). To simulate, forecast or just nowcast space weather events, which, from a mathematical standpoint, are essentially time-dependent processes, global models are needed to simulate, as the initial stage, the steady state of the solar-terrestrial environments. The resulting solution for the "pre-event" solar wind and terrestrial magnetosphere serves as the background through which the space weather disturbance propagates and shocks and discontinuities are generated. The solar energetic particle (SEP) transport as well as the formation of their seed population is mostly controlled by this background solution, as well.
In realistic solar atmosphere simulations, the -essentially three-dimensional (3-D) -interplanetary magnetic field can be steady-state only in coordinate systems corotating with the Sun (such as HGR). In corotating frames solar active regions and coronal holes, that are mostly responsible for shaping the structure of the solar atmosphere, can be treated as steady-state sources. Despite being steady-state globally, the solution for the solar wind is highly variable at Earth's location, since in any corotating systems the Earth orbits the Sun with the Carrington rotation period. This relative motion allows us to compare time series of solar wind parameters observed by Earth or L1 orbiting spacecraft, SW observed (t), with the simulated steady-state parameters for the 3-D background solar wind, SW background (R), by extracting the time series of simulated values, SW background (R Earth HGR (t)), at the time-dependent Earth location in the HGR system, R Earth HGR (t). The time series of observations by the two STEREO spacecraft provide additional validation possibilities.
In the last two decades several 3-D solar wind (Usmanov et al. 2000;Suzuki & Inutsuka 2005;Verdini et al. 2009;Osman et al. 2011;Lionello et al. 2014a,b;Usmanov et al. 2018), and coronal heating models (Tu & Marsch 1997;Hu et al. 2000;Dmitruk et al. 2002;Li & Habbal 2003;Cranmer 2010) that included, or were exclusively driven by, Alfvén wave turbulence became increasingly popular. Our group also includes Alfvén wave turbulence in our corona and inner heliosphere model (cf., Sokolov et al. 2013Sokolov et al. , 2021aOran et al. 2013;van der Holst et al. 2014;Gombosi et al. 2017Gombosi et al. , 2021. Although popular, this physics-based approach to modeling the background solar wind is not the only way to model the solar corona and the solar wind. Empirical descriptions of the solar wind, like the widely used Wang-Sheeley-Arge (WSA) model (Arge & Pizzo 2000) are also attractive because of their simplicity and ability to predict the solar wind speed in the inner heliosphere. In addition, the WSA formalism can be readily incorporated into global 3-D models for the solar corona and inner heliosphere, e.g., via a varying polytropic index distribution (see Cohen et al. 2006), as proposed by Roussev et al. (2003), or via the boundary condition set at the solar surface (Shen et al. 2018) or at 0.1 AU (see, e.g. Pomoell & Poedts 2018).
There is, however, a conceptual difference between global solar atmosphere models that are based on numerical solutions of the ideal (or resistive) MHD equations and simple (but efficient) semi-analytic solar wind models that imply a classical Parker-spiral structure for the interplanetary magnetic field lines. Specifically, in semi-analytic models it is assumed that in the corotating frame of reference under steady-state conditions the solar wind velocity vector u(R) is always aligned (parallel or anti-parallel) with the interplanetary magnetic field vector, B(R). This means that stream-lines and magnetic field lines are always identical and they form Archimedean or Parker (1958) spirals. The assumed coincidence of solar wind stream-lines and interplanetary magnetic field lines (in corotating frames) is at the heart of almost all analytic, semi-analytic, semi-empirical and empirical models of interplanetary space. This assumption results in simple methods for tracing interplanetary magnetic field lines from any point of interest in the solar system back to the Sun. One can just trace back a Parker spiral to find out where the local solar wind and magnetic field line originates from on the Sun. For instance, in the WSA semi-empirical model (Arge & Pizzo 2000) this assumption is used to relate the Earth's location to the point on the source surface at which to extract the parameters, predicted by the model.
The problem with the concept of aligned interplanetary stream-lines and magnetic field lines is that numerical solutions of the full set of MHD equations have not been able to reproduce this feature so far. Full MHD numerical solutions of the solar atmosphere obtained within the framework of regular MHD are not stream-aligned. Even if the initial conditions describe parallel B and u vectors, the smallest numerical error (such as truncation, finite representation, etc) will result in an uncontrollable growth of misalignment. One of the reasons for this discrepancy is the numerical reconnection across the heliospheric current sheet near the top of helmet streamers: the reconnected field is directed across the current sheet, while the global solar wind streams along the current sheet, thus resulting in significant misalignment. Within regular MHD there is no mechanism to re-establish the stream-line-field line alignment.
Tracing magnetic field lines from a point of interest in the heliosphere to the solar surface (magnetic connectivity) is needed to identify the point of origin of the magnetic field line. This same magnetic connectivity is also necessary to solve the field-aligned transport of energetic particles that create space radiation hazards. With standard computational MHD calculating the magnetic connectivity can become so challenging that sometimes it is preferable to trace stream-lines instead, or even Lagrangian trajectories of fluid elements.
In this paper we present a new global model of the solar atmosphere/solar wind system using a reduced MHD model which ensures stream-aligned magnetic field under steady-state conditions. We call this solution "stream-aligned MHD" (SA-MHD). In this reduced set of MHD equations the electric field, E = −u × B ≡ 0, identically vanishes in the induction equation and so does the electromagnetic energy flux (Poynting flux) in the energy equation, E × B ≡ 0. At the same time the impact of the magnetic force on the plasma motion (the Ampère force) is explicitly enforced in the momentum equation. It must be emphasized that any steady-state solution of the reduced equation set is a solution of the full MHD system of equations. However, the converse statement is not true, since in an arbitrary steady-state magnetic field the electric field does not have to vanish identically (though its curl must). The main point is that the practically interesting Parker (1958) spiral solutions for the steadystate background solar atmosphere and solar wind can be efficiently generated within the SA-MHD with no loss of generality, in 3-D geometry with realistic solar magnetic field.
where B = |B|, u = |u|, andÎ is the unit tensor. We use the equation of state, P = 2n i k B T , ρ = m p n i for the coronal plasma with the polytropic index, γ = 5/3. We omitted many important effects (effects of gravity and inertial forces are discussed in section 3.3, see Sokolov et al. 2021a, for detailed description of a model used in simulations here) in the governing equations, while we explicitly included terms proportional to the divergence of magnetic field. If the numerical solution for the magnetic field is not divergence-free (see Godunov 1972;Powell et al. 1999), this approach allows us to handle maintain basic physical properties, such as pressure positivity. Indeed, if one takes Eq. (1d) divided by temperature, T , subtracts the dot product of u/T with Eq. (1c) and the dot product of B/(µ 0 T ) and Eq. (1b) and finally adds Eq. (1a) times s + 1 T u 2 2 − h , an extra conservation law is obtained: (see also Godunov 1961Godunov , 1972Harten et al. 1983;Powell et al. 1999), where s and h are the entropy and enthalpy per unit mass and the thermodynamic equations were used. Eq. (2) is only valid for smooth solutions and it states that the entropy satisfies the conservation law for a passively advected scalar. At shocks and discontinuities the entropy conservation breaks down, but, entropy can only increase, not decrease. Therefore, the total entropy (its volume integral) does not decrease: ∂ t ρsdV ≥ 0, since the integral of ∂ t (ρs) over the smooth solution region vanishes due to Gauss' theorem applied to Eq. (2), while the discontinuities result in entropy production. A consequence of this entropy property is that for such systems finitevolume numerical schemes ensure that in the numerical solution the physical entropy does not decrease, and consequently, the pressure can never become negative (should the pressure go to zero, the entropy, s ∝ log(P ), would drop to minus infinity).

Separating the Potential Magnetic Field in the Solar Corona
Realistic models of the 3-D coronal magnetic field use boundary conditions taken from full disc photospheric magnetograms incorporating current and past observation results. The potential magnetic field solution provides the minimum free magnetic energy for a given boundary condition, therefore, in the "ambient" (background) solution for the solar wind the magnetic field is approximately equal to the potential configuration in the proximity of the Sun. Following Ogino & Walker (1984) and Tanaka (1994) (see also Powell et al. 1999;Gombosi et al. 2002) we split the total magnetic field to a potential (B 0 ) and non-potential (B 1 ) component B = B 0 + B 1 . We note that at R = R the potential B 0 field dominates. If we use photospheric maps of the radial magnetic field (B r ), then the potential B 0 field can be recovered from the observed magnetogram using the Potential Field Source Surface Method (PFSSM) that has been originally described by Altschuler et al. (1977). The Laplace Eq. for scalar magnetic potential is solved in terms of spherical harmonics between R ≤ R and R ≤ R SS = 2.5R with the given radial gradient of the potential (the observed radial field) at R = R and with vanishing magnetic potential (i.e. purely radial magnetic field) at R = R SS . Accordingly, the non-potential (B 1 ) field (used by Sokolov et al. (2013); Oran et al. (2013); van der Holst et al. (2014)) satisfies zero boundary condition for the radial field component, Using the split magnetic field approach the governing equations become the following (the continuity equation does not change, therefore we do not repeat it here): An advantage of splitting the magnetic field is that there are no quadratic terms in B 0 in equation (4b). Near the Sun (especially near active regions) such terms can be so large, that the error in their approximation might exceed the whole physical effect from the B 1 field. For the same reason we cannot, in general, rely on the assumption that numerical approximations for ∇ × B 0 and ∇ · B 0 are exactly zero (as it is theoretically assumed), so that to maintain entropy conservation we need to introduce the corresponding source terms in equations (4a-4c). Particularly, all contributions from the field to the momentum Eq. (4b) can be combined with the and it is the last term in the RHS of (4b), which ensures that the potential, B 0 , contributes to the total field, but not to the current, in the Ampère force.
In a solar corona model different approximations, with split and non-split magnetic fields, can be used depending on the heliocentric distance, R = |R|. For instance, the "threaded field line model" (Sokolov et al. 2021a) that includes the transition region and covers heliocentric distances between the photosphere (R ) and the low boundary of the solar corona, R b ≈ 1.1R , neglects the non-potential B 1 component, assuming that in this region B 1 ≡ 0, and the plasma flow occurs only along the potential field lines (threads). This approach allows us to move the boundary conditions for all other physical quantities from the top of the transition region (where they are poorly known) to the photosphere. Accordingly, the boundary condition, ( Finally, above the magnetogram source surface, at R ≥ R SS , the B 0 field, which is purely radial at R SS , can be naturally continued beyond the source surface as a steady-state, radial and divergence-free field decaying as 1/R 2 : This magnetic field, however, is no longer a potential field, therefore, its contribution to the current cannot be neglected. To fully account for the total force effect when using split fields and to ensure their continuity at R SS , one needs to add the extra source term to the RHS of the momentum Eq. (4b) in the R ≥ R SS region: which is easy to compute, since only the last term (which is to cancel the last term in the RHS of 4b) depends on time via B 1 , while all other terms are steady-state terms and can be calculated only once.
With the newly added source the momentum equation, in effect, reduces to Eq. (1c) expressed in terms of the total magnetic field, B = B 0 +B 1 and the induction equation (Eq. 4a) reduces to (1b). Nevertheless, it is the unknown field B 1 which is solved from these equations throughout computational domain. In order not to break the entropy conservation, for R ≥ R SS the corresponding source term, S E = u · S M , should be added to the RHS of the energy Eq. (4c). For derivations, but not for computations, the RHS of resulting energy equation can be written as follows: Although the energy Eq. (4c) with the split field and the transformed RHS is formally reducible to Eq. (1d) for the full field, such reduction would require us to redefine the magnetic energy density in terms of the total field too, as ( , balances the last term in the RHS of Eq. (7). This redefinition, however, may compromise the whole idea of only solving equations for B 1 beyond R SS . To use the energy equation, (4c), with an extra source, u · S M , seems to be a better option.

Intuitive solution
Now, we show how to proceed from standard MHD equations (1a-1d) to equations for SA-MHD, in which where α(t, R) is a scalar function of time and coordinates. To derive a governing equation for α, in a steady-state stream-aligned solution, one can integrate the MHD equations over a magnetic flux tube element of a length of d bounded by two nearby cross sections of the flux tube, dS 1 and dS 2 , and a bundle of magnetic field lines about a chosen magnetic field line, which all pass through the contours of these cross sections. The equation, ∇ · B = 0, gives BdS = const (9) along the flux tube, which allows us to relate the change in the cross-section area along the thread to the magnetic field magnitude. Now, since the stream-lines are aligned with the magnetic field lines, one can also integrate the continuity Eq. (1a) for a steady-state flow along the flux tube, which gives Therefore, the ratio of constant mass flux and the constant magnetic flux, is constant along the field line and along the stream-line. Both of these relations can be expressed as Now, if we assume that this ratio is constant along the fluid particle trajectory (which in steady state coincides with the stream-line but in the dynamic system this statement is a matter of a model assumption): we can combine equations (1a) and (13) to obtain the conservation law for α: 3.2. Rigorous solution A less intuitive, but mathematically more rigorous derivation of the same set of model equations can be obtained if we consider an MHD state that is close to the force equilibrium, ∂ t u → 0, and add the aligning source (that tends to zero) to the induction equation (Eq. 1b): where the definition of α from Eq. (8) is generalized for the case of non-aligned vectors u and B; for the aligned vectors both definitions are identical. To keep the entropy conservation, a similar source should be added to the RHS of the energy (Eq. 1d): Now, if in an arbitrary MHD initial state the magnetic field is aligned according to Eq. (8), the aligning term ensures, that and keeps the alignment forever. By substituting B = αu into Eq. (15), we arrive at equation (14) (multiplied by u).
In the momentum equation the stream-aligned magnetic field should be applied: which can be also written in a non-conservative form: Eq. (19) contains the explicit expression for the Ampère force, which happens to be always perpendicular to the velocity vector (similarly to the Coriolis force). In the equation of energy (1d) the Poynting flux, E × B vanishes identically. In addition, the sources in the RHS of Eq. (16) for the stream-aligned field reduce to ∂ t B 2 /2µ 0 and entirely balance the change in the magnetic energy density in the LHS. Therefore in effect, the stream-aligned magnetic field does not contribute to the energy density, so that the energy equation does not include the magnetic field at all: This is because the Ampère force is perpendicular to the velocity vector and it does not affect the kinetic energy, hence, does not exchange energy between the field and plasma. For the same reason, under steady-state conditions, equation (20) conserves the Bernoulli integral along the stream-line in its classical form: Equations (1a,14,18,20) present a full set of governing equations, which, as the solution approaches steady state, becomes fully conservative (the nonconservative source term that is proportional to ∂ t α in the momentum equation tends to zero), satisfies the ∇ · B = 0, constrain (since, ∇ · B = −∂ t α → 0), and their solution in this limit presents a legitimate solution of the MHD equations, since the aligning sources (∝ ∂ t u → 0) tend to zero too. Moreover, this is a desired stream-aligned solution. On the other hand the full MHD equations with no aligning sources, once applied to any non-equilibrium MHD state, with ∂ t u = 0, will immediately produce nonstream-aligned solution even from stream-aligned initial state. It seems that in any discretized implementation (actual numerical solution) it is extremely difficult (if not impossible) to reach stream-aligned steady-state 3-D solutions within the framework of regular computational MHD, which justifies the approach proposed here.

Stream Aligned MHD in a Rotating Frame of Reference
For applications to solar wind the momentum and energy equations of SA-MHD are used the frame of reference rotating with angular velocity, Ω , with respect to the inertial frame. In addition to the Coriolis force density, −2ρΩ × u, the centrifugal and gravitational forces are accounted for via the gradient of the potential, Φ(R): In the rotating frame the Bernoulli integral (Eq. 21) is modified by the inertial force potential (Eq. 22): The Bernoulli integral can be used to relate the model predictions at the source surface close to the Sun, such as the WSA predictions at R = 2.5 R , to the observed solar wind parameters at 1 AU (see e.g., Cohen et al. 2007). However, Eq. (25) is not convenient for use, since the velocity in the inertial frame, u inert , relates to that in the rotating frame, u, with the well-known equation: so that at 1 AU the velocity in the rotating frame has too large contribution from rotational velocity. In terms of the observable parameters in the inertial frame, the Bernoulli integral implies the conservation of the following quantity: Within the framework of isothermal approximation for solar atmosphere, one can use generic expression enthalpy per a unit of mass, h, instead of a particular model for ideal gas with constant γ assumed in Eq. 27. For isothermal two-component plasma with T i = T e = T , the expression for enthalpy is h = 2k B T mp log ρ. Note, that in neglecting the effect from the gas-kinetic pressure and from the magnetic field, for a test particle of a unit mass co-moving with the solar stream velocity, u, both the energy integral, u 2 inert 2 − GM R = const, and the vector of particle angular momentum, [R × u inert ] = const, conserve separately and so does their linear combination in Eq. 27. However, the force effect of stream-aligned magnetic field increases the angular momentum of the solar wind while moving outward the Sun. Therefore, according to Eq. 27 at large heliocentric distances there is a finite gain in the solar wind energy per unit of mass equal to Ω · [R × u inert ] 1 AU The closed system of equations describing self-consistent variation of the angular momentum and azimuthal magnetic field component while the solar wind propagating outward the Sun (Parker 1958) can be obtained by solving the equations of SA-MHD in spherical coordinates, R, θ, ϕ, which are the heliocentric distance, co-latitude and longitude correspondingly. The analytical solution may be found assuming no dependence on ϕ (axially symmetric field) and no motion over θ (radially divergent flow). The Coriolis force forces the particles to rotate in the negative ϕ direction, while it has no component in the θ direction. Consequently, the fluid motion occurs only in the (R, ϕ) direction and all stream-lines will be on conical surfaces of constant θ. Steady-state stream-aligned solutions of Eqs. (1a, 23, 24) may be integrated along stream-lines yielding following equation: In the absence the magnetic field, i.e. at α ≡ 0, Eq. (28) would express the conservation of angular momentum (R sin θu ϕ,inert ≡ |R × u inert |) in the inertial (not rotating!) system. However, due to the effect from realistic interplanetary magnetic field the heliocentric distance in the right hand side should be taken at the Alfvén point where the local stream speed, u, equals to the Alfvén wave speed, Typically, the factor M −2 A is very large near the base of the solar corona, because in this region the solar magnetic field is large and the solar wind not accelerated yet. In the opposite limiting case at large heliocentric distances, where the solar wind is fast and the magnetic field is weak, the plasma stream is hyperalfvénic, M −2 A is negligibly small. In between the accelerating and rarefying solar wind must go through the Alfvén point and this is the point at which the integration constant is determined (this point to avoid discontinuous solutions). With these considerations the azimuthal velocity is (see Weber & Davis 1967):   Parker (1958) spiral by about 9 • at Earth's orbit. Taking into account the solar rotation (with respect to Earth) this means that a stream-aligned field line originating from the same point of the solar surface will cross Earth about 15 hours later than the corresponding Parker (1958) spiral aligned field line. In space weather applications, tracing the stream-line in the co-rotating frame from the observation point toward its origin site can be used to associate in situ solar wind properties to remotely observed structures of the solar corona (see, e.g. Biondo, Ruggero et al. 2021). The large difference demonstrated in Fig. 1 means that field line associated phenomena (like SEP events) will arrive to Earth some 15 hours later than predicted by simple Parker (1958) spiral connectivity models. It is, therefore, important that our computational model fully incorporates all principles and capabilities of earlier models.

Numerical Solution with Full MHD
For comparison with the SA-MHD numerical result presented later in this paper, we provide numerical solution for the steady state solar corona prior to the SEP event occurred on April 11, 2013 (Lario et al. 2014). GONG magnetogram at 2013 Apr 11, 06:04 UT serves as the inner boundary condition and the main driver of the solar wind model (see the left panel of Fig. 2). In the simulation, the initial condition of the 3-D magnetic fields in the solar corona is reconstructed by a potential field source surface (PFSS) solution using a spherical harmonic expansion with the inner boundary provided by the magnetogram. The radial magnetic fields on the surface of the Sun calculated using PFSS is shown in the right panel of Fig. 2. The PFSS reconstruction has been used to obtain the potential (B 0 ) field below the source surface located at 2.5R . Below the source surface we solve the full set of MHD equations (1a and 4a-4c), above it we add the source term S M given by Eq. (6) to the momentum equation, (equation 4b), as well as the extra energy source, u · S M , to the energy equation, (equation 4c). The governing equations for the solar corona and inner heliosphere models described by Sokolov et al. (2021a) are solved using the SWMF (Tóth et al. 2005;Toth et al. 2012) developed at the University of Michigan. The magnetic field in the equatorial plane ( Fig. 3 top panels) and in the y = 0 meridional plane of the rotating HGR coordinate system (Fig. 3 bottom panels) exhibits visible imperfections in the form of disconnected ("V-shaped") field lines below 20R and too long closed loops in (the middle bottom panel of Fig. 3). The magnetic field in the equatorial plane of the inner heliopshere (right panels in Fig. 3) is totally disordered and is not applicable to solve magnetic connectivity. We emphasize that this feature is independent of the numerical solution algorithm: while steadystate field lines and stream-lines are identical at the PDE level, this is not true in any discretization within the framework of computational MHD. The stream-aligned model is designed to address all these issues.

The Leontovich Boundary Condition as the Aligning Source
The computational SA-MHD cannot be applied to the entire solar corona that starts with a slowly expanding atmosphere with speeds as low as a few km/s on top of the transition region and in closed field regions. The slow-speed region extends to heliocentric distances of 2.5 − 3.5 R . In this region the SA-MHD would fail to describe closed magnetic field lines. In addition, the problem to align the strong magnetic field with the slow and somewhat random motion is not only physically inappropriate, but also mathematically ill-posed, since the characteristic perturbation speeds tend to infinity if the plasma speed tends to zero. Therefore, we successfully apply the SA-MHD above some spherical heliocentric boundary with a radius of R SA .
Below this boundary we apply source terms in the MHD equation, which are proportional to the electric field describing the non-alignment between the field and flow and turning to zero, when the alignment is achieved. These terms reduce the electric field below R SA and entirely eliminate it at R = R SA . However, they do not prevent slow motions and magnetic field line reconnection and closure. The sources are heuristically derived from the well-known Leontovich (1948) boundary condition which make them physics-based.
The Leontovich boundary condition (Leontovich 1948;Landau et al. 1985) is applicable to electromagnetic fields when a high-frequency wave interacts with a metal and deep inside the metal the field vanishes. The currents shielding the field are assumed to be concentrated in a thin layer near the boundary, having the impedance, Z. Under these assumptions the transverse components of electric and magnetic fields outside the metal are related to boundary condition: wheren is a unit vector normal to the metal surface directed inward the metal. A straightforward (however, insufficient) application to the boundary between the usual and SA-MHD gives us an opportunity to consider the electric-field-free SA-MHD region as if this is the metal and to assign the impedance Z to the boundary surface. In this case at the MHD side of this boundary the finite electric field will cause a change in the B 0 field given by Eq. (31): It is important to emphasize that this field is induced by the current which is not conducted by the moving solar wind plasma, that is why Eq. (32) describes B 0 field, rather than B 1 . Another comment is that the field undergoes a jump across the surface at which the finite current is concentrated, that is why the field given by Eq. (32) exists only "outside the metal," that is in the SA-MHD region that starts at a heliocentric distance of R = R SA . The modified field contributes to the flux function resulting in the evolution of the MHD state toward higher alignment (lower electric field). The effect is easy to evaluate if the impedance surface is orthogonal to the magnetic field. In this case the electric field, E = B nn × u t is orthogonal tô n. The normal components, B n = B ·n, u n = u ·n of the field and velocity are not affected by the extra field given by Eq. (32), while the transverse velocity and field in the layer of a thickness of ∆x, orthogonal to the magnetic field evolve since the extra field Eq. (32) is present at the upper boundary of the slab and absent at the lower one. Using the flux functions in equations (4a) and (4b), one can find the source terms for the MHD velocity, (∂ t u) L , and for the induction field, (∂ t B) L , caused by the shielding currents implied in the Leontovich (1948) boundary condition (that is why we use subscript L): in the last equation one can also put B n u n = B · u, since the normal vector is aligned with the field, hence, B = B nn . Now, one can derive aligning sources assuming the distributing conducting layers instead of concentrated impedance by substituting d(1/Z)/dx for 1/(Z∆x) ratio, which may be conveniently parameterized as d(1/Z)/dx = 1/(V 2 A + u 2 n )τ , so that: These sources at some choice of τ can be added to the MHD equations in a chosen region, not necessarily near the boundary. The advantage of the chosen parameterization is that τ directly characterizes the electric field relaxation due to the aligning sources, since: Although greatly improving the MHD solution at R < R SA , the source terms do not allow us to align the solution at R → R SA . The latter goal is achieved if while integrating numerically the MHD equations over time after advancing the solution through the time step, ∆t, we add the source terms and assume that the relaxation rate, 1/τ , is equal to 1/∆t at R → R SA and modulated with some geometric factor, ξ, which equals one at R → R SA and decays at smaller heights. In other words, the linear aligning operator is applied to the velocity and magnetic field vectors in the following manner: The total energy density decreases by the following quantity: At ξ = 1, that is in the proximity of the SA-MHD computational domain, operator 38, 39 provides a perfect alignment about some weighted average direction such that the weak field is aligned with the strong stream without noticeable modification in the stream direction and, in the opposite limiting case, slow flow is aligned with the strong field.

Numerical Solution with SA-MHD
The simulation results with SA-MHD are provided for the same solar configuration as described in Section 4.1. Below the PFSS source surface at R = 2.5 R the MHD model is not modified and the numerical solution is not much different from that provided in Fig.2. In the intermediate region 2.5 R < R < 3.5 R we apply the aligning operator 38, 39 based on the Leontovich boundary condition as described in Section 4.2. The geometric factor ξ = R/R − 2.5 is chosen equal to zero at the lower boundary and one at the upper boundary, of this region.
At R = R SA = 3.5 R the MHD solution is enforced to be aligned, so that above this boundary the SA-MHD equations are solved numerically. The characteristic perturbation speeds needed to construct the numerical scheme are discussed in Appendix. Fig. 4 shows the magnetic field structure in the corona and inner heliosphere obtained with SA-MHD and the Leontovich (1948) boundary condition. One can see that the unphysical features seen in Fig. 3 are gone: There are no long closed loops or V-shaped field lines and the magnetic field lines in the inner heliosphere clearly form a Parker (1958) spiral.

Magnetic Connectivity
Magnetic connectivity (the magnetic field line connecting a point of interest in the heliosphere to the solar photosphere) is a very important component of space weather simulations. Among other things magnetic connectivity controls the simulated intensity and time profiles of SEP events (cf., Borovikov et al. 2018;Li et al. 2021;Young et al. 2021). For the synoptic map shown in Fig. 2 the stream-lines are not as simple as in the case of uniform solar wind speed (see Fig. 1). As one see in Fig. 5 the delay between the Earth crossing of MHD and SA-MHD stream-lines can vary between a few hours and a few days, depending on actual angular radial speed profile.
Most existing SEP models use MHD stream-lines as proxies for IMF lines and solve the SEP transport along these stream-lines. Our results indicate that this approach should be used with great care: the Earth crossing of an interplanetary stream-line can be off by a few hours or a few days depending on actual solar wind conditions. In addition, the sign of the radial magnetic field component can flip back and forth depending on the heliospheric current sheet crossings, so the µ = cos(pitch− angle) quantity (appearing in the focused transport equation of Skilling (1971)) might abruptly change sign along the stream-line. An equation describing energetic particle transport in SA-MHD and a method to solve such equation numerically is described in a follow-up paper (Sokolov et al. 2021b).

CONCLUSION
We showed how the MHD equations can be modified to describe stream-aligned steady-state flows of plasma parallel to the magnetic fields with no electric field. We call this approach stream-aligned MHD (SA-MHD). We applied the new approach to modelling the solar corona and heliosphere in a coordinate frame co-rotating with the Sun. Domains using the SA-MHD approach can gradually be merged to regions in which the ordinary MHD equations are solved by distributed external currents similar to those implied by the Leontovich (1948) boundary condition. With this approach, the numerical solution in the heliosphere is perfectly aligned with the IMF lines shaped as classical Parker spirals. Our model addresses the magnetic connectivity problem as well as the simulation of the SEP acceleration and transport. This work was supported by the NASA Heliophysics DRIVE Science Center (SOL-STICE) at the University of Michigan under grant NASA 80NSSC20K0600. We would also like to acknowledge high-performance computing support from: (1) Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation, and (2) Pleiades operated by NASA's Advanced Supercomputing Division. In order to demonstrate an effect of the aligning terms on the characteristic waves, one can linearize Eqs. (15) and (16), together with Eqs. (1a) and (1c). For linear characteristic perturbations, (δρ, δu, δP, δB) T , propagating along the x axis, i.e. depending on the coordinates and time via the combination, (x − Dt). The background for the wave propagation is (ρ 0 , u 0 , P 0 , B 0 ) where B 0 = (B x0 , B y0 , 0). The z component of the magnetic field is made zero by choosing a coordinate system in which the magnetic field vector is in the (x, y) plane. In the equations for the perturbation, all time derivatives should be substituted with ∂ t = −D∂ x , and, after this substitution, one can omit the x derivatives by denoting ∂ x (δρ) → δρ etc. Note, that De Sterck et al. (1999) analyzed eigenvalues of two-dimensional steady-state SA-MHD, however, such eigenvalues only characterize the direction of wave propagation, rather than their speed needed to construct numerical schemes.
Equations (A1h) and (A1d) can be solved separately: This pair of waves can be considered in two limiting cases. In the first limiting case we consider "traditional" MHD with no alignment source. Now one can put α = 0,Ṽ A,x = 0 in Eq. (A4). Note, that for B x0 > 0 and u x0 > 0 the perturbations of velocity and magnetic field are parallel in the left propagating wave (Λ < 0), but they are anti-parallel in the wave propagating to the right, creating or increasing field misalignments with the stream. Both of these Alfvén waves can exist in the absence of an aligning source and propagate with the Alfvén speed, V A,x = B x0 / √ µ 0 ρ 0 .
In the MHD description without the aligning sources this would be the slow magnetosonic wave. We see that in MHD (two 3-component-vector equations + two scalar equations, 8 characteristic waves) with the aligning sources and stream aligned background two kind of waves (right Alfvén and right slow magnetosonic) degenerate. These are the only waves in which the perturbations are not stream aligned. The distinction of SA-MHD (one 3-component-vector equation + three scalar equations with 6 characteristic waves) is that it lacks these two degenerate waves since this approximation does not allow non-aligned perturbations.
In the remaining three branches of waves the magnetic field is stream-aligned and there is a relationship between δu y and δB y : In the particular case of magnetic field aligned propagation (B y0 = 0) this reduces to the expression for another Alfvén wave of different polarization: (W 3a) : Λ 3 = −Ṽ A,x , δ(ρ, u x , u y , u z , P, α) T = (0, 0, 1, 0, 0, 0) T δW 2 ,