Boundary Treatment for the Subsonic/Alfvénic Inner Boundary at 2.5 R ⊙ in a Time-dependent 3D Magnetohydrodynamics Solar Wind Simulation Model

A new magnetohydrodynamics (MHD) simulation model of the global solar corona and solar wind is presented. The model covers the range of heliocentric distance from 2.5 solar radii, so that coronal mass ejections at the earliest phase near the Sun can be treated in the future. This model is constructed by introducing a characteristics-based boundary treatment to an existing heliosphere 3D MHD model. In tailoring a set of characteristic equations for this new model, we assume that the coronal magnetic field is open to interplanetary space and that the solar coronal plasma is flowing outward everywhere at 2.5 solar radii. The characteristic equations for the subsonic/Alfvénic inner boundary surface are satisfied by altering the plasma density and/or temperature to maintain a polytropic relationship. In this article, the details of the characteristics-based boundary treatment for the middle of the corona (named CharM) are provided. The quasi-steady states of the solar wind derived from simulations with various choices of a parameter in the boundary treatments are compared and examined. Although further improvements are needed, we apply the new boundary treatment to simulations for three Carrington rotation periods from the minimum to maximum phase of the solar activity cycle, and show that an optimal choice yields a reasonable quasi-steady state of the transonic/Alfvénic solar wind matching the specified subsonic/Alfvénic plasma speed at 2.5 R ⊙.


Introduction
3D time-dependent magnetohydrodynamics (MHD) simulations are a powerful tool in the field of solar physics and space weather research. The simulation approach can trace numerically the propagation of coronal mass ejections (CMEs) from the innermost heliosphere to interplanetary space to the position of the Earth (e.g., Han et al. 1988;Smith & Dryer 1990;Dryer et al. 1991;Odstrcil et al. 2004;Hayashi et al. 2006;Shen et al. 2011Shen et al. , 2014Wold et al. 2018;An et al. 2019;Singh et al. 2022;Wu et al. 2022). Before simulating interplanetary disturbances such as CME propagation, the ambient quiet solar wind must be prepared. To obtain the quasisteady state of the solar wind as a proxy of the ambient solar wind, time-relaxation MHD simulation approaches have been widely used (e.g., Steinolfson et al. 1982;Han et al. 1988;Linker et al. 1990;Usmanov 1993;Wang et al. 1993;Washimi & Sakurai 1993;Usmanov & Dryer 1995;Hayashi 2005;Feng et al. 2010Feng et al. , 2021Hayashi et al. 2022aHayashi et al. , 2022b. The solar wind originates from the solar surface base of the coronal hole as subsonic/Alfvénic flow. The plasma flows are heated and accelerated, and the flow becomes supersonic/Alfvénic at some distance from the solar surface. A time-relaxation approach with a time-dependent MHD simulation model can solve the nonlinear interactions between the coronal magnetic field and plasma and numerically reproduce the transonic/Alfvénic solar wind in a physics-based manner. The same MHD model for the time-relaxation simulation can be used to simulate the temporal evolutions of CME events as responses to numerical perturbation mimicking the initial phase of the CME.
MHD simulation studies for the solar corona, solar wind, and disturbance propagations, in general, have used one of two choices in setting the inner boundary surface. The first choice is to place the inner boundary surface at 1 R e or the base of the corona on which the MHD variables, such as the radial components of the magnetic field (B r ), can be straightforwardly specified as realistic input boundary values from observations (e.g., Usmanov 1993;Usmanov & Dryer 1995). The other is to place the inner boundary surface in a supersonic/Alfvénic region. By setting the inner boundary sphere well beyond the critical points (i.e., at r = 15 ∼ 20R e ), simple boundary conditions such as the fixed boundary condition are allowed. In this case, it is necessary to introduce assumptions to infer the MHD boundary values on the inner boundary surface from observations and/or models (e.g., . It is often desired to simulate the solar corona and solar wind in a range of the heliocentric distance from the middle of the corona. For example, one can use information about CME events at an initial phase from coronagraph observations and/ or magnetic flux-rope models (e.g., Chen 1996;Thernisien et al. 2006;Wood & Howard 2009) to straightforwardly simulate the evolution of the CME event. Because the information is often about CMEs after their launches from the lower corona, the heliocentric distance of the simulated region is not necessarily set from 1R e in this example. In addition, we want to exclude the lowermost part of the solar corona, as the Courant-Friedrich-Levy condition is severe for regions with large Alfvén speeds and steep radial gradients of plasma quantities. By excluding the lowermost corona, we can allocate more computational resources to solving the temporal evolutions in the solar wind in distant regions. Therefore, an MHD model covering the region from the upper part of the solar corona to interplanetary space can bring benefits to space weather studies.
The plasma flow in the solar corona is overall subsonic/ Alfvénic. In the subsonic/Alfvénic regions, the temporal variations of MHD variables are determined by the states of the neighboring regions in the lateral (latitudinal and longitudinal) directions and the radial direction. Because some of the necessary information is missing, MHD simulations with the subsonic/Alfvénic boundary surface often suffer computational difficulties, such as unphysical numerical vibrations and instabilities near the boundary surface that grow exponentially in time and eventually cause computational failures. Characteristics-based boundary treatments (e.g., Nakagawa et al. 1987;Wu & Wang 1987) are a powerful mathematics-based tool to compensate for the lack of information about the subsonic/ Alfvénic boundary surface. In this boundary treatment, the relationship among the temporal variations of the MHD variables is treated in a manner such that the characteristic equations of the MHD equations and the specified boundary constraints will be satisfied simultaneously (e.g., Hayashi 2005).
Recently, we constructed a new subsonic/Alfvénic boundary treatment for MHD simulations with the inner boundary sphere set at the height of the middle of the corona (henceforth named CharM). We implement the CharM boundary treatment in an existing MHD simulation model, the heliosphere 3D MHD (H3DMHD) model (Han et al. 1988;Wu et al. 2007Wu et al. , 2015. In this new model, the heliocentric distance of the inner boundary sphere is chosen to be 2.5 solar radii (R e ). It is reasonable to assume that the radial component of the plasma flow is always positive (V r > 0) on the inner boundary surface at 2.5R e , because the solar coronal plasma appears to be flowing outward nearly radially at this distance. Under this assumption, the CharM boundary treatment can be simple, because it does not have to consider the case of stagnant plasma (V r ∼ 0). The heliocentric distance of the inner boundary sphere coincides with the radius of the outer boundary sphere or the source surface of the potential field source surface (PFSS) model (Schatten et al. 1969;Altschuler & Newkirk 1969). In this article, the two variables, V r and B r , on the inner boundary surface will be fixed.
The characteristic method plays another role than providing computational stability: it seeks an optimal set of unspecified plasma temperatures and/or densities such that will be compatible with the specified plasma speed (V r ) and the magnetic field B r . Because it is difficult to determine the plasma temperature and density in the middle of the corona, this feature can be a useful numerical tool. Test simulations are conducted for demonstrating these capabilities and the limitations of this model. This article is organized as follows: details of the timedependent 3D MHD model and the characteristics-based boundary treatment (CharM) are given in Sections 2 and the Appendix. The simulation results are shown in Section 3. A summary and discussion are given in Section 4.

MHD Model
The H3DMHD (Wu et al. 2007) model is a 3D timedependent MHD simulation model for the global solar wind.
The model solves time-dependent 3D MHD equations in the Sun-centered spherical coordinate system: where ñ, V, B, P g ,  , r, t, and g are the mass density, the plasma flow velocity, the magnetic field vector, the gas pressure, the energy density (=ñv 2 /2 + P g /(γ − 1) + B 2 /8π), the position vector originating at the center of the Sun, time, and the solar gravitational force (−GM ☉ r/r 3 ), respectively. The colons (:) denote the dyadic products of vectors. The plasmas are assumed to consist of fully ionized hydrogen, and the plasma pressure and plasma density are related with P g = 2ñk B T p /m p , where k B , m p , and T p are the Boltzmann constant, the mass of protons, and the plasma temperature, respectively. For the test simulations, we here set the specific heat ratio γ = 5/3. The current model does not include the effects of coronal heating; hence the solar wind plasma quantities in distant regions may not be realistic. Nonetheless, we here choose this value because γ will have to be equal to 5/3 when the heat source term is included in the coronal simulations in the future. This value of the specific heat ratio yields reasonably good agreements in the CME propagation simulations (e.g., Liou et al. 2014).
The H3DMHD model solves the MHD equations in the rest frame. In the simulations shown in this article, 110 grids are placed along the radial direction to cover the heliocentric distance from 2.5 R e to 19.0 R e , with a constant grid size of Δr = 0.15R e . We note here that the heliocentric distance of the upper boundary surface can be set farther away, at 1 au or beyond. We choose a rather small number to reduce the computation times. The angular sizes of the grids in the latitudinal and longitudinal directions are both uniform and set to be 5°. In order that the effect of the solar rotation will be included, the boundary maps of specified MHD variables (B r , V r and ñ or T p ) will be longitudinally shifted at the angular velocity of the solar sidereal rotation, Ω, of 360°.0 per 24.3 days.

Characteristic Boundary Treatments
The set of the practical forms of the equations is given in the Appendix. In brief, under the concept of the characteristics for the hyperbolic eight-variable MHD equation system, we are allowed to specify five constraints at the 2.5 R e inner boundary surface, because there are at least five MHD wave modes directing outward from the Sun (or inward to the simulation domain). Among several possible choices, we choose the one as follows.
Two sets of conditions specify four of the five constraints: the first set is that the radial components of the magnetic field (B r ) and plasma flow (V r ) are unchanged in the frame rotating with the Sun at the solar sidereal angular velocity (Ω). Here we choose to fix V r to test the present boundary treatment. The second set is that the plasma and magnetic field are parallel when seen in the rotating frame, in order to satisfy the requirement for magnetic solenoidality (Yeh & Dryer 1985). The first set is expressed with two equations for the temporal variations of the vector components, ∂ t V r = 0 and ∂ t B r = 0. The second condition is expressed as In this work, the remaining one constraint is given as a polytropic relationship, ∂ t (P g /ñ α ) = 0. The exponent (α) is the only parameter to control the CharM boundary treatment in this study, and it is not intended to represent any physics in the lower part of the solar corona. Hence, the exponent (α) is not necessarily equal to the specific heat ratio in the MHD equations (γ = 5/3). For example, with α = 1, the CharM boundary treatment alters the boundary density and simultaneously keeps the temperature constant in time, ∂ t T p ∝ ∂ t (P g / ñ) = 0. Setting α equal to a large number (such as 100) is equivalent to setting a condition of the (near-)constant density, Setting α equal to zero is equivalent to setting a condition of the (near-)constant plasma pressure, although we do not test this case in this study. In the following sections, we conduct the simulations with various values of α -1, 1.2, 1.4, 5/3, 3.0, and 100.0. The value α = 5/3 is chosen because the value is the same as the specific heat ratio (γ) in the governing MHD equations. The values 1.2 and 1.4 are between the two cases α = 1 and α = 5/3. The value 100 is chosen as a proxy of an infinitely large value of α, and the value 3.0 is set for an example case between α = 100 and α = 5/3. Hereafter, α is called the boundary-polytrope index.
In the Appendix, the characteristic equations to be solved are given. It must be mentioned here that the characteristic method is flexible and that equations other than those shown in the Appendix are possible.

Boundary Values at R 2.5 ☉ and Initial Values
In this study, we use the data of the source surface   which generally yields a weaker source surface field than anticipated from the 1 au in situ measurements. To compensate for the magnetic flux, typically a factor of 5 is sufficient, if we assume that no flux cancellation takes place in interplanetary space. In this study, we choose a factor value of 10, to run the codes with a slightly stronger magnetic field strength for testing.
We apply the speed prediction formula, (kilometers per second), to infer the radial component of the plasma speed (V r ) at r = 18R e from f s . Figure 1(c) shows the inferred V r at 18 R e . The initial and boundary values of V r in this simulation study are estimated from the inferred 18 R e values through a linear relationship with the height from 1 R e , V r , , where r¢ is the heliocentric distance in units of solar radius.
The initial plasma number density n (=ñ/m p ) is determined as n n V V r 215 ¢. The bottom-boundary plasma temperature (T p ) is determined from an assumption that the Notice that in each column, the overall shapes well resemble each other, but the values are substantially different. A smaller (greater) value of α yields higher (lower) boundary plasma density and lower (higher) boundary plasma temperature. scaled sum of the kinetic energy, gravitational potential, and thermal energies, r 2 {ñV 2 /2 − ñGM ☉ /r + 2k b ñT p /m p /(γ − 1)}, is constant and equal to the average at 1 au. For this estimation, the average mass density, plasma speed, and temperature at 1 au (n 1au , V 1au , and T 1au ) are set to 8.0 (count cm −3 ), 420.0 (kilometers per second), and 10 4 K, respectively. The initial value of the plasma density and the temperature above the inner boundary surface at r > 2.5R e are reduced to 50% and 10%, respectively, so that the initialized solar wind will start flowing outward smoothly. Above the inner boundary surface at r > 2.5R e , the initial values of the latitudinal and longitudinal components of the magnetic field and plasma velocity (B θ , B f , V θ , and V f ) are set to zero at t = 0.
We assume that the inner boundary sphere (at r = 2.5R e ) is rotating rigidly at the rotation rate Ω, as the coronal magnetic field below this height is sufficiently strong to control the plasma flow. The longitudinal component of the plasma flow at r = 2.5R e is given as V r sin q = W f , where θ is the colatitude counted from the solar north pole. The parallel condition between the magnetic field and the plasma flow must be satisfied in the frame rotating with the Sun; hence the longitudinal component of the magnetic field is given as . The latitudinal component of the magnetic field is given as B θ = B r V θ /V r . The CharM boundary treatment uses the MHD variables as they are at the boundary grids, except that the longitudinal component of plasma flow in the rotating frame (V r sin q -W f ) is used instead of V f in the nonrotating simulation frame.
The boundary variables that will not be altered through the CharM boundary treatment (such as B r and V r ) are calculated with the linear longitudinal interpolation of the original map data at each simulation time step. The variables altered through the CharM boundary treatment (X) are further altered as ∂ t X = − Ω∂ f X with the upwind differencing method, to take into account the effect of the solar rotation.
It is worth noting that the relationships among the components of B and V can be set rather arbitrarily. Any choice that satisfies the induction equation, where V is defined in the frame rotating with the boundary map of the fixed B r , can be allowed here. This is a condition for maintaining the solenoidality of the simulated magnetic field with the boundary B r (Yeh & Dryer 1985).

Results from Time-relaxation Simulations
The simulated solar wind system evolves from the initial state in accordance with the MHD equations and the CharM boundary treatment. Figure 2 shows the profiles of the radial component of plasma velocity (V r ) in the radial direction at a selected location of the latitude of 0°and the Carrington longitude of 180°at several physical times simulated with α = 1 (fixing T p at 2.5R e ), as an example. As seen in Figure 2, the initial plasma flow immediately starts evolving, to eventually reach a (quasi-)steady state. The typical relaxation time under the present simulation settings is found to be about 15 to 20 hr on the physical timescale. We regard the simulated state at t = + 20 hr as a steady time-relaxed state.

Time-relaxed States with Various Values of Boundarypolytrope Index, α
In the present model, the boundary-polytrope index, α, controls the behavior of the CharM inner boundary treatment. Figure 3 shows boundary maps of the plasma density (in the left column) and temperature (in the right column) obtained with α = 1, 1.2, 5/3, and 100, at t = 20 hr.
As designed, the boundary temperature with α = 1 ( Figure  3(e)) is not altered from the initial setting of the boundary value. Similarly, the boundary density with α = 100 ( Figure  3(d)) has little changed from the initial setting. In Figures 3(a)-(c) and (f)-(h), the distributions are rather blurred, because of the numerical error in transporting in the longitudinal direction as ∂ t ñ = − Ω∂ f ñ or ∂ t T p = − Ω∂ f T p . Overall, the shapes of the simulated boundary ñ and T p are similar with different α. However, the values are substantially different: the difference in the boundary density between the cases with α = 1 (Figure  3(a)) and α = 100 (Figure 3(d)) is about a factor of 4. The difference between the boundary temperature with α = 1 (Figure 3(e)) and α = 100 (Figure 3(h)) is about a factor of 5. During the time relaxation, the boundary treatment with α = 1 (α = 100), in general, increases the boundary density (temperature) to achieve a sufficiently large plasma pressure on the boundary surface so that the boundary plasma can flow outward at the specified speed under the presence of the nonuniform boundary magnetic field. The lateral gradient of the boundary magnetic field pressure can result in the direction of the simulated magnetic field being slightly diverted from the radial direction. The nonradial magnetic field above the inner boundary surface can obstruct the smooth plasma outward flows, and the influence of the obstruction can propagate backward to the inner boundary surface, because the plasma flows near the inner boundary surface are yet subsonic/ Alfvénic. Without the characteristics-based boundary treatment adjusting the plasma pressure (P g ∝ ñT p ), the simulated solar wind may not reach the steady state, as shown in Section 3.2.
The higher plasma density at 2.5 R e results in higher total mass flux with the fixed V r . Because the initial values of the mass density were calculated so that the total mass will well match those measured at the Earth (or 1 au), such an increase is not indeed desired. Similarly, the higher boundary temperature results in unrealistic high solar wind speeds in distant regions. Figure 4(a) shows the speed profile in the radial direction at the selected location of the latitude of 0°(at the solar equator) and the Carrington longitude of 180°, derived with various different values of α. Substantial differences in the speed near the outer boundary (19 R e ) are seen.
In Figures 4(b) and (c), the profiles of V r in time-relaxed states are shown for two other periods, CR 2059 and 2095. These two periods, CR 2059 and 2095, correspond to the minimum phase between the solar activity cycles 23 and 24 and the ascending phase of solar cycle 24, respectively. These two periods are chosen in order to examine whether the present model can handle the transonic/Alfvénic solar wind at periods other than CR 2126 (near the maximum phase of solar cycle 24 or at the first peak of the sunspot number in solar cycle 24). As seen in the plots of Figure 4, the present MHD model can indeed yield the steady state of the transonic/Alfvénic flow through the time relaxation (for 20 hr on a physical timescale) for 2.5R e r 19R e for these two periods as well. The shapes of the obtained profiles of V r shown in Figures 4(b) and (c) are similar to those in 4(a) (for CR 2126), except that the values of V r are dependent on the boundary values specified for each period and their gradients in the longitudinal and latitudinal directions.
Although it is not perfectly satisfactory, among the tested cases, the choice of α = 1.2 appears to yield most reasonable results. Figure 5 shows the latitude-longitude maps of the simulated variables on the inner boundary surface at r = 2.5R e and the outer boundary surface at r = 19.0R e . As seen in Figure 5(e), the flow speed (V r ) at r = 19.0R e has moderate contrast, ranging from about 350 to 480 km s −1 . The velocity contrast is a reasonable one for the quiet solar wind, although further examinations and test simulations are needed. Figure 6 shows the profile of the plasma flow speed (V r ) in the radial direction from three simulation cases, sampled at a location (S22.5°, 180°): (a) the (quasi-)steady state derived with α = 1 (thick line) as a reference; (b) the case with fixed boundary conditions for all eight MHD variables (without the characteristics-based boundary treatment); and (c) the same as the second case, except that the temperature on the inner boundary surface is multiplied by 5.

Results with and without the CharM Boundary Treatment
As seen in Figure 6, the case without the characteristicsbased boundary treatment, but with the same initial boundary values as in case (a), is unstable, yielding even negative values of V r . The selected location for the radial profile was near the source of the heliospheric current sheet and the lateral (latitudinal and longitudinal) gradient of B r is relatively large. We think that this magnetic configuration is one of the factors causing the falling plasma in the second case (b). The plasma flow with insufficient thermal energy (temperature) or density cannot pass through the oblique narrow paths near a laterally expanding flow region due to the gradients of the magnetic pressure. The falling plasma cannot be steady. The current setting (enforcing V r > 0 on the inner boundary surface) allows the numerical vibrations move outward, rather than stay at a certain region. Hence the vibration may not grow exponentially in time. This allows us to continue simulation (b), but the results may not be usable for other simulations, such as CME propagation simulations.   1.01, 1.1, 1.2, 1.3, 1.46, and 5/3) in the MHD equations. The case with γ = 1.01 corresponds to the near-isothermal case. The value of 1.46 is the one inferred from the Helios data analysis (Totten et al. 1995). The boundary values of V r at 2.5 R e are fixed and identical among these cases. The radial profiles are sampled at the latitude of 0°(the solar equator) and the Carrington longitude of 180°.
In the reference case (a), the CharM boundary treatment alters (usually increases) the mass density, to maintain the steady flow that matches the specified V r on the inner boundary surface. Without such a (numerical) mechanism, we cannot obtain the steady-state matching, given that V r has been rather arbitrarily determined. Profile (c) shows that the transonic/Alfvénic simulation with the fixed boundary conditions but a higher bottom-boundary temperature can yield a smooth V r profile. However, the speeds in the distant regions are about 900 km s −1 , which are unreasonably high as a quiet solar wind speed. From these results, we claim that the present CharM boundary treatment can reasonably reduce the chance of obtaining unreasonable or unstable solutions. This is an important advantage, because it is difficult to know in advance whether a certain combination of the plasma density, temperature, and velocity will result in a stable smooth flow or not. The Parker solution (Parker 1958(Parker , 1965) with a larger specific heat ratio will infer such combinations of the boundary values that are usable as the subsonic/Alfvénic fixed boundary values, but only for the case without the presence of the magnetic field.

Summary and Discussion
A characteristics-based boundary treatment, named as CharM, is introduced. This boundary treatment model is designed to provide numerical stability and robustness to MHD simulations for the transonic/Alfvénic solar wind, starting from a middle height in the corona. Results from time-relaxation simulations with various values of the boundary-polytrope index (α), a free parameter in the boundary treatment, are compared and examined to find that α = 1.2 is a tentatively optimal number in the present model.
The characteristics-based boundary treatment is designed to assist transonic/Alfvénic MHD simulations by adjusting the boundary temperature and/or density to obtain a reasonable steady-state solution matching the specified boundary velocity and magnetic field. Insufficient thermal energy given/specified on the inner boundary at 2.5 R e can lead to the downward flow seen in the simulation with the fixed boundary condition. It is found that the CharM boundary treatment can alter the plasma values within the specified polytropic constraint, to achieve a steady state of the solar wind flow at r > 2.5R e . Because the plasma quantities at the middle of the corona are not well determined, this functionality is an important model feature for assisting Sun-to-Earth simulation models, such as those for CME propagations in the upper corona and interplanetary space.
There are several items to be solved and improved in the present CharM boundary treatment. The present model with α = 1 tends to yield a steady-state solar wind with higher temperatures at 2.5R e . Similarly, the runs with a large value of α (=100) tend to yield a steady-state solar wind with a higher boundary plasma density.
The optimal boundary temperature and/or density sought by CharM must be different if we include coronal heating and/or acceleration in the governing MHD equations. One of the next steps is to include additional energy and momentum sources or to apply a smaller value of γ, such as those inferred from in situ measurements (γ = 1.46; Totten et al. 1995), to mimic the thermal energy supply.
In Figure 7, the velocity profiles of the steady state derived with various values of the specific heat ratio (γ) are shown. In this comparison, the boundary-polytrope index (α) is set to 5/3, and the boundary values of V r are the same. The smaller value of the specific heat ratio (closer to 1) yields a higher speed near the outer boundary. The values of V r of the nearisothermal case (γ = 1.01) at the outer boundary surface are more than twice as large as those with γ = 5/3. The profiles in the region close to the inner boundary surface are rather similar to each other, though the profile with γ = 5/3 becomes rather flat for r > 5R e , while the profiles with smaller γ keep increasing with respect to the heliocentric distance. The plasma density and temperature on the inner boundary surface in the time-relaxed steady states differ with different values of γ as a consequence of the nonlinear interaction between the inner boundary and the simulated solar wind above the inner boundary surface during the time-relaxation simulation. The interactions will be more complicated if we include the source terms of the energy and/or momentum in the governing MHD equations. It is an advantage of the present model that it is capable of studying such nonlinear interactions between the subsonic/Alfvénic boundary surface and the transonic/ Alfvénic solar wind above the boundary surface, although careful simulation setups are necessary.
As the present model adjusts the boundary plasma temperature and density in a manner to maintain the polytropic relationship, P g /ñ α , the final boundary temperature and density of the steady state are dependent on the initial guesses. To prepare initial values of the boundary plasma temperature, we simply assumed that the temperature and density are functions of the plasma velocity that is inferred from the PFSS model and the flux tube expansion factor ( f s ). New methods of inferring the plasma quantities at the middle of the corona are desired. The present model can be improved with such new input information, which can in turn help to check the appropriateness of the model.
In this article, we fix the boundary V r then construct the characteristic equations shown in the Appendix. The values of the specified 2.5 R e V r are rather small, in part because we want to examine how well CharM can handle the boundary MHD variables in a situation far away from being supersonic/ Alfvénic. It is possible to construct sets of characteristic equations where the boundary V r is allowed to change. For example, we can assume that the constant mass flux and the temporal variations of the boundary V r and ñ can be related as ∂ t (ñV r ) = ñ∂ t V r + V r ∂ t ñ = 0, instead of Equation (A2). Starting with this, we can construct a new set of characteristic equations, replacing the equation set (A7). We will test other possible choices as well, as the concept of the characteristics of the hyperbolic MHD system allows us to find suitable choices for retrieving desired solar wind features, such as good agreements of the mass flux at the Earth or reasonable velocity profiles with respect to the heliocentric distance near the Sun.
where L l,m is the element of the left eigen matrix of the MHD equation system, for lth wave modes (numbered in nonincreasing order of eigenvalues, and mth variables (numbered in a order: ñ, V r ñ, V θ ñ, V f ñ, B θ , B f , and  ). The characteristic equation for the radial component of magnetic field, ∂ t B r = − V r ∂ r B r + L, can be isolated from the seven other characteristic equations. The components of the left eigen matrix (L l,m ) and the eigenvalues (λ l ) used in this study are from Cargo & Gallice (1997) and also given in Hayashi (2005). The right-hand side of the mth MHD variable is calculated in the same manner as the grids above the inner solar surface boundary sphere. In this present simulation study, we want to fix the radial component of plasma flow (V r ) and magnetic field (B r ). The specified V r is always set to be positive, V r > 0. In this case, four of seven eigenvalues, λ 1 = V r + V F , λ 2 = V r + V A , λ 3 = V r + V S , and λ 4 = V r , are always positive. The number of outgoing wave modes with negative eigenvalues, λ l < 0, representing the modes propagating from the domain of computation toward the Sun, is equal to or less than three. Therefore, we are allowed to specify at least five constraints to complete the equation system for determining the temporal variations of all eight MHD variables.
We first consider a case where all three remaining eigenmodes (with λ 5 = V r − V S , λ 6 = V r − V A , and λ 7 = V r − V F ) are negative. In this study, we choose the five constraints as follows. The first condition is to fix B r in time, ∂ t B r = 0. This condition must be accompanied by two constraints, V r B θ − V θ B r = 0 and V r B f − V f B r = 0. We want to keep the boundary V r distributions, hence the fourth condition is expressed as ∂ t V r = 0. The fourth condition is expressed as With this, the temporal variations of B θ and B f can be written as As the last condition, we choose the polytropic relationship between the plasma density and pressure (or temperature), which is controlled with the boundary-polytrope index α. Setting α = 1 is equivalent to setting a fixed boundary condition for the temperature (∂ t T p ∝ ∂ t (P g /ñ) = 0). Setting α = ∞ is equivalent to setting a fixed boundary condition for the plasma density (∂ t ñ = 0). By setting α = 0, the boundary surface gas pressure is unchanged in time (∂ t P g ). In this work, we only test the cases with α 1. It is worth noting that the relationship is not necessarily an exponent polytropic one: we introduce the polytropic relationship for simplicity. Indeed, any relationship between the plasma density and pressure (or temperature), ∂ t P g ≔ F∂ t ñ, where F is any function, can be used for constructing the characteristics-based boundary treatment.
With the equations above, the temporal evolution of the energy density (the left-hand side of Equation (4)) is rewritten as By coupling these constraints among the temporal variations and characteristic equations (Equation (A1)), we finally obtain a set of three characteristic equations with three unknown variables, ∂ t ñ, ∂ t (ñV θ ), and ∂ t (ñV f ): for the three modes with eigenvalues, λ 5 = V r − V S , λ 6 = V r − V A , and λ 7 = V r − V F . This equation set is always solvable. When λ 5 > 0, λ 6 > 0, or λ 7 > 0, we are allowed to specify more than five constraints. It is possible to design the boundary treatment to switch a set of characteristic boundary treatments for the five constraints (with five positive eigenvalues) to/from those for six, seven, or eight positive eigenvalues; however, it is rather complicated, and such switching is often very frequent and causes numerical instability. In this present model, we simply set the total contribution of an incoming wave mode to zero when its eigenvalue is positive (the wave is incoming and directed outward from the Sun): for l = 5, 6, 7. It is evident that when the boundary flow is super-Alfvénic, all RHS l will be zero; hence, all of the three temporal variations, ∂ t ñ, ∂ t (ñV θ ), and ∂ t (ñV f ), will be calculated to be zero (equivalent to the fixed boundary condition). The denominators in the equations shown in this Appendix are never equal to zero. The radial component of plasma flow, V r , is the only variable that can be zero in general cases; however, the present model specifies always positive V r (outward plasma flow) at r = 2.5R e . The plasma density, ñ, is always positive unless the vacuum is considered. A parameter, (γ − 1), is not equal to zero unless the isothermal case is considered.