Role of Longitudinal Waves in Alfv\'en-wave-driven Solar Wind

We revisit the role of longitudinal waves in driving the solar wind. We study how the the $p$-mode-like vertical oscillation on the photosphere affects the properties of solar winds under the framework of Alfv\'en-wave-driven winds. We perform a series of one-dimensional magnetohydrodynamical numerical simulations from the photosphere to beyond several tens of solar radii. We find that the mass-loss rate drastically increases with the longitudinal wave amplitude at the photosphere up to $\sim 4$ times, in contrast to the classical understanding that the acoustic wave hardly affects the energetics of the solar wind. The addition of the longitudinal fluctuation induces the longitudinal-to-transverse wave mode conversion in the chromosphere, which results in the enhanced Alfv\'enic Poynting flux in the corona. Consequently, the coronal heating is promoted to give higher coronal density by the chromospheric evaporation, leading to the increased mass-loss rate. This study clearly shows the importance of the longnitudinal oscillation in the photosphere and the mode conversion in the chromosphere in determining the basic properties of the wind from solar-like stars.


INTRODUCTION
Inspired by the observed double-tail structure of comets, which indicates the presence of gas outflow (Biermann 1951), Parker (1958) predicted the outward expansion of the hot coronal plasma, which results in the formation of transonic outflow. Later on, the insitu measurement by the Mariner 2 Venus probe confirmed the existence of supersonic plasma streams from the Sun, which is now called the solar wind (Neugebauer & Snyder 1966). Hot coronae and stellar winds are also ubiquitously observed in low-mass main sequence stars that possess a surface convection zone (Wood et al. 2005(Wood et al. , 2021Güdel et al. 2014;Vidotto 2021) In the framework of the thermally-driven wind model, the energy source of the solar wind is the thermal energy of the solar corona. The thermally-driven wind model therefore predicts that faster solar wind emanates from hotter regions of the corona, and vice versa. In reality, however, several observations indicate that highspeed solar wind emanates from relatively cool parts of the corona. Fast solar wind is known to originate from coronal holes (Krieger et al. 1973;Kohl et al. 2006), which exhibit cooler temperature than the other regions of the corona (e.g., Withbroe & Noyes 1977;Narukage et al. 2011). The observed anti-correlation between the freezing-in temperature and the velocity of the solar wind (Geiss et al. 1995;von Steiger et al. 2000) also supports the fact that fast solar wind originates in cool portions of the corona. These observations indicate that magnetic field plays a substantial role in the solar wind acceleration.
It is believed that the convection beneath the photosphere is the source of the energy for the hot corona and the solar wind (Klimchuk 2006;McIntosh et al. 2011). Convective fluctuations excite various modes of waves that propagate upward (Lighthill 1952;Stein 1967;Stepien 1988;Bogdan & Knoelker 1991). Magnetic reconnection between open and closed field lines is another possible source of transverse waves (Nishizuka et al. 2008), in addition to the direct ejection of heated plasma (Fisk 2003). Among various types of waves, Alfvén(ic) waves have been highlighted as a reliable agent to effectively transfer the kinetic energy of the convection to the corona and the solar wind via the Poynting flux (e.g., Belcher 1971;Shoda et al. 2019;Sakaue & Shibata 2020;Matsumoto 2021). This is first because they are not so much affected by the shock dissipation owing to the incompressible nature, unlike compressible waves, which easily steepen to form shocks as a result of the amplification of the velocity amplitude in the stratified atmosphere, and second because they do not refract, unlike fast-mode magnetohydronamical (MHD hereafter) waves (e.g., Matsumoto & Suzuki 2014), but do propagate along magnetic field lines (Alazraki & Couturier 1971;Bogdan et al. 2003).
In recent years transverse waves have been detected in the chromosphere (Okamoto et al. 2007;De Pontieu et al. 2007;McIntosh et al. 2011;Okamoto & De Pontieu 2011;Srivastava et al. 2017), whereas it is still under debate whether the sufficient energy required for the formation of the corona and the solar wind propagates into the corona (Thurgood et al. 2014).
In contrast to Alfvénic waves, acoustic waves have not been considered to be a major player in the coronal heating because the acoustic waves that are excited by p-mode oscillations at the photosphere (e.g., Lighthill 1952;Felipe et al. 2010) rapidly steepen to form shocks before reaching the corona (Stein & Schwartz 1972;Priest 2014;Cranmer et al. 2007) . However, Morton et al. (2019) pointed out the contribution of p-mode oscillations to the generation of Alfvénic waves via the mode conversion from longitudinal waves to transverse waves (Cally & Hansen 2011). The aim of the present paper is to investigates roles of the acoustic waves that are excited by vertical oscillations at the photosphere in the Alfvén wave-drive wind. For this purpose, we perform MHD simulations that handle the propagation, dissipation, and mode-conversion of both transverse and longitudinal waves from the photosphere to several tens of solar radii with self-consistent heating and cooling.
In Section 2 we describe the setup of our simulations. We present main results in Section 3. We discuss related topics in Section 4 and summarize the paper in Section 5.

METHODS
We consider the magnetohydrodynamics of the solar wind in one-dimensional (1D hereafter) open magnetic flux tubes from the photosphere at r = R (solar surface) to several tens of solar radii. Figure 1 shows an overview of our model.

Basic Equations
We consider a one-dimensional (spherical symmetric, ∂/∂θ = ∂/∂φ = 0), super-radially expanding flux tube. The cross section of the flux tube is defined by the filling factor of the open flux tube f op (r), which is lower than unity on the photosphere and asymptotically approaches unity as r gets larger. The conservation of the open magnetic flux Φ op yields the following relation.
where X represents the value of X on the photosphere. We note that Φ op is constant in each simulation. We solve the one-dimensional MHD equations along the flux tube characterized by f op (r). For simplicity, we consider the polar wind, which is not affected by the solar rotation. In deriving the MHD equations in a super-radially expanding flux tube, the scale factors of the coordinate system are required. Here, we assume that the magnetic flux tube expands isotropically in θ and φ directions. In terms of scale factors, this assumption yields Using these scale factors, the MHD equations in an expanding flux tube is derived (see Shoda & Takasao (2021) for derivation). The basic equations are given as follows. 1 where v, B, ρ and p are velocity, magnetic field, density and gas pressure, respectively. v ⊥ and B ⊥ are the perpendicular (θ and φ) components of v and B, respectively, that is, where e θ and e φ are unit vectors in θ and φ direction, respectively. M is the solar mass. e denotes the total energy density per unit volume given by where e int is the internal energy density per unit volume. p T denotes the total pressure: represent the phenomenological turbulent dissipation of Alfvén waves (see Section 2.3 for detail).
Q C and Q R represent the conductive heating and radiative cooling per unit volume, respectively. In terms of conductive flux q cnd , Q C is given by For q cnd , we employ the Spitzer-Härm type conductive flux (Spitzer & Härm 1953) that strongly depends on temperature and transports energy preferentially along the magnetic field line. Besides, to speed up the simulation without loss of reality, we quench the conductivity in the low-density region. q cnd is then employed as follows.
The radiative cooling rate is given by a linear combination of optically thick and thin components as follows.
where Q thck R and Q thin R correspond to the optically thick and thin radiative cooling rates, respectively. The control parameter ξ rad should satisfy ξ rad ≈ 1 in the photosphere and ξ rad ≈ 0 from above the transition region. Although the profile of ξ rad is given as a solution of radiative transfer, here we simply model it as follows.
where we set p chr = 0.1p . Thus the radiation is assumed to be optically thick in p p chr and optically thin in p p chr .
Following Gudiksen & Nordlund (2005), we approximate the optically thick cooling by an exponential cooling that forces the local internal energy to approach the reference value e ref int : where we set the time scale as follows.
whereρ = 1.87 × 10 −7 g cm −3 , is the mean (timeaveraged) mass density in the photosphere. The reference internal energy is calculated once the corresponding reference temperature T ref (r) is given. Here, we set T ref (r) = T . The optically-thin cooling function is composed of two different contributions. In the chromospheric temperature range, we employ the radiative cooling function given by Goodman & Judge (2012) (Q GJ ), while in the coronal temperature range, the loss function Λ(T ) is given by the CHIANTI atomic database.

Equation of state
The hydrogen in the lower atmosphere (photosphere and chromosphere) of the Sun is partially ionized because the temperature there is not sufficiently high. In this work, the effect of the partial ionization is considered in the equation of state. The internal energy is composed of the random thermal motion of the particles and the latent heat of the hydrogen atoms, which is given by where n H is the number density of hydrogen atoms, χ is the ionization degree and I H = 13.6 eV is the ionization potential of hydrogen. For simplicity, the formation of H 2 molecules is not considered. The thermal equilibrium is assumed with respect to ionization, in which the ionization degree is given by the Saha-Boltzmann equation.
where λ e is the thermal de Broglie wavelength of an electron: Note that pressure and ionization degree are connected by

Phenomenology of Alfvén-wave turbulence
In heating the solar wind, energy cascading is required to convert the kinetic and magnetic energies to heat by viscosity and resistivity. Alfvén-wave turbulence, a type of MHD turbulence which is triggered by the collision of counter-propagating Alfvén waves (e.g., Goldreich & Sridhar 1995;Lazarian 2016), is a promising process for the energy cascading in the solar wind. Because Alfvén wave turbulence is a three dimensional process, and thus, to deal with the turbulent dissipation in the one-dimensional system, one needs to model the effect of turbulence. Here we adopt a phenomenological model of Alfvén wave turbulence (Hossain et al. 1995;Dmitruk et al. 2002;van Ballegooijen & Asgari-Targhi 2016), which yields the (averaged) turbulent heating rate in terms of mean-field quantities (Elsässer variables). Following Shoda et al. (2018a), the turbulent dissipation terms in Eq.s (5) and (7) are explicitly given by and where λ ⊥ is a perpendicular correlation length and z ± ⊥ are Elsässer variables (Elsasser 1950) defined by We assume that the correlation length increases with the radius of the flux tube.
Because the Alfvénic fluctuations are localized in intergranular lanes on the photosphere (Chitta et al. 2012), we set λ ⊥ as a typical width of inter-granular lanes.
The dimensionless coefficient c d is chosen following van Ballegooijen & Asgari-Targhi (2017) as

Geometory of Flux Tubes
In modeling the solar wind in a one-dimensional flux tube, we need to prescribe the filling factor of the open flux tube f op as a function of r. Since the open flux tube is localized on the photosphere and expands as the radial distance increases, f op (r) should be an increasing function of r that asymptotically approaches unity.
Following Shoda et al. (2020), we employ the two-step expansion of the flux tube, which is described in terms of f op (r) as where f exp 1 (r) and f exp 2 (r) represent the first and second expansions, respectively.
The first expansion occurs in the chromosphere until one flux tube merges with the adjacent flux tube. Although the direct observation of chromospheric magnetic field is still missing, because the expansion occurs in response to the exponential decrease in the ambient gas pressure, it would be straightforward to assume that the filling factor increases exponentially in height. For this reason, the following formulation is adopted.
where f op cor is the open-flux filling factor in the corona and H mag is the scale height of flux-tube expansion. We relate H mag to the pressure scale height on the photosphere H by where a = 6.9 km s −1 and g = 0.274 km s −2 are the sound speed and the gravitational acceleration on the photosphere, respectively.
The second expansion occurs in the extended corona with a typical length scale of ∼ R . Following Kopp & Holzer (1976), we adopt the coronal expansion as where We adopt the fixed values of r exp /R = 1.3, σ exp /R = 0.5. The values of f op and f op cor are summarized in Table  1.

Simulation Setup
The simulation domain extends from the photosphere (r = R ) to the outer boundary (r = r out ) located at nearly r = 100R in most cases. The radial distance of r out for each run is tabulated in Table 1. At r = r out , we set the free boundary conditions. A great advantage to set the inner boundary at the photosphere is that we can self-consistently calculate the density at the coronal base, which is one of the critical parameters to determine the mass loss rate,Ṁ w , of the solar wind (e.g., Lamers & Cassinelli 1999). The coronal base, where the density is nearly ten orders of magnitude smaller than the density at the photosphere, is frequently set as the inner boundary of simulations for solar and stellar winds (Verdini et al. 2010;Lionello et al. 2014;Shoda et al. 2019). However, the coronal-base density is determined by the chromospheric evaporation as a result of the energy balance between conductive heating and radiative cooling at the transition region (Rosner et al. 1978;Withbroe 1988). Specifically, when the heating in the corona increases, denser chromospheric material is heated up by the thermal conduction from the corona, resulting in an increase in the density at the coronal base. Since our numerical simulations solve these heating and cooling processes in a self-consistent manner, we can obtain re-liableṀ w independently from the treatment of the inner boundary at the photosphere.
The size of the spatial grid, which varies with r, is set as follows: (34) where we set ∆r m = 20 km, ∆r M = 2000 km, ε ge = 0.01, and r ge = 1.04R . Figure 2 shows ∆r as a function of r.
At the inner boundary, we fixed the temperature to the photospheric value, The initial temperature is set to T = T in the entire simulation region. We initially set the hydrostatic density distribution with T = T in the inner region that is extended with a power-law profile in the outer region: where we adopt ρ w,0 = 10 −19 g cm −3 unless otherwise stated. The inner hydrostatic profile switches to the outer power-law one at r/R − 1 ≈ 0.01. We note that although the outer density is larger than the hydrostatic value with T = T , it is still smaller than the observed density in the solar corona and wind by a factor of five. The transverse components of velocity and magnetic field correspond to the amplitudes of Alfvénic waves. The inner boundary condition of them are defined in terms of the Elsässer variables (Eq. (25)) in the photosphere. We set the free boundary condition to the incoming component at the inner boundary so that it is absorbed without being reflected there: To inject MHD waves from the photosphere, we impose time dependent boundary conditions for the density, velocity and perpendicular magnetic field. The transverse perturbation is injected via the outgoing component of the Elsässer variables with a broadband spectrum, where φ t N is a random phase and The longitudinal perturbation, which originates from the p-mode oscillation, is excited with the density, where ρ = 1.88 g cm −3 , and the radial velocity, where φ l N is a random phase and This corresponds to the period range between 100 seconds and 5 minutes, which is narrower than that of the transverse component. We tabulate the transverse and longitudinal components of the root-mean-squared velocity amplitudes, δv ⊥, and δv , , of the input fluctuations at the photosphere in Table 1. We take δv ⊥, = δv , = 0.6 km s −1 as standard values for the velocity perturbation on the photosphere. We here note that, because the only outward flux is selectively injected from the photosphere in both transverse and longitudinal components, the corresponding "random" velocity amplitude is about √ 2 times larger than these values, which are comparable to observed transverse (Matsumoto & Kitai 2010) and longitudinal (Oba et al. 2017a) amplitudes of ∼ 1 km s −1 .
Each model is labeled BxVyy, where "x" indicates the type of the magnetic flux tube and "yy" denotes the amplitude of δv , . We classify the cases into three groups by the effect of the magnetic field. The first group, which includes only one case, is labeled x = 0. In this case, we switch off Alfvénic waves by setting δv ⊥, = 0; we test whether the formation of the corona and wind is possible or not only by longitudinal waves. The result of this case is presented in Appendix A (see also Suzuki 2002).
The second and third groups are labeled x = s and w, which stand for "standard (or strong)" and "weak" magnetic fields, respectively. The aim of these groups is to investigate how the longitudinal-wave excitation on the photosphere affects the properties of the solar wind.
For this purpose, we compare cases with different amplitudes of δv , for a fixed transverse amplitude of δv ⊥, = 0.6 km s −1 . In the second group, we adopt the equipartition magnetic field, B = 1300 G, at the photosphere from 8πp /B 2 = 1 (Section 3) to model observed kilo-Gauss patches (Tsuneta et al. 2008;Ito et al. 2010). In the third group, in order to examine the effect of the geometry of magnetic flux tubes on the propagation and dissipation of waves in the chromosphere, we reduce B to 1/4 of that of the second group with keeping the field strength above the corona (Section 4.1).
In order to extensively investigate the effect of longitudinal waves on the solar wind, the second group in particular is investigating a wide range of 0 < δv , < 3.0 km s −1 . δv , 2 km s −1 , which is larger than the observed average value explained above, targets transient large-amplitude disturbances (e.g., Oba et al. 2017b).
We perform the simulations for a sufficiently long time in order to study the average behavior of the atmosphere and wind after they reach a quasi-steady state. To satisfy this requirement, the simulation time is set to 4500 minutes for the cases with δv , = 0 − 1.2 km s −1 and 6000 minutes for the cases with δv , = 1.5 − 3.0 km s −1 . Even after the quasi-steady state is achieved, the radial profile fluctuates in time. Therefore, when we compare average properties of different cases, we take the average of physical quantities for 1500 minutes before the end of the simulation. Time-averaged wind profiles of cases with δv , = 0 km s −1 (blue dashed; BsV00), 0.6 km s −1 (green solid; BsV06) and 1.8 km s −1 (red dotted; BsV18) in comparison with observations. (a): Temperature. The circles (Fludra et al. 1999) show the radial distribution of electron temperature observed by CDS/SOHO. (b): Radial velocity. The circles (Teriaca et al. 2003) and the stars (Zangrilli et al. 2002) represent proton outflow speeds in polar regions observed by SOHO. The location of the top of the chromosphere at T = 2 × 10 4 K for each case is shown by diamonds in both panels. Time-averaged density profiles in the chromosphere and the low corona in comparison to observations. The line types are the same as those in Figure 3. The squares represent electron density (right axis) obtained from observations of multiple total solar eclipses during solar minimum phases (Saito et al. 1970). The crosses and triangles (Wilhelm et al. 1998) are electron density observed by SOHO in interplume lanes and plume lanes, respectively. Open and filled diamonds show the location of the top of the chromosphere at T = 2 × 10 4 K and the location of the coronal base at T = 5 × 10 5 K, respectively. δv , cases, the terminal velocity is nearly invariant with δv , . This shows that the variety in the solar wind velocity is unlikely to come from the longitudinal-wave injection from the photosphere. Figure 4 shows the radial profiles of the mass density ρ (left axis) in the chromosphere and corona (0.005 ≤ r/R − 1 ≤ 1), in comparison to the observed electron densities n e (right axis) in the corona. In converting n e to ρ, we assume that the corona is composed of fully ionized hydrogen plasma, that is, ρ = m H n e . The line format is the same as that of Figure 3. In contrast to the temperature and velocity, the density depends significantly on δv , . Specifically, the coronal density is four times larger in δv , = 1.8 km s −1 than in δv , = 0.0 km s −1 . Given that the filling factor of the open flux tube f op is fixed and the wind velocity is nearly independent from δv , , the larger coronal density means the larger mass-loss rateṀ w , which is given byṀ Our simulation results show that the wind mass-loss rate is sensitive to the longitudinal-wave injection. The underlying physics is discussed in detail in the following sections.

Mass-loss rate and Wind Energetics
To see more quantitatively how the mass-loss rate depends on the photospheric longitudinal-wave amplitude δv , , we present in Figure 5 the relation between δv , and the mass-loss rate (blue stars) evaluated at r = r out ; shown on the right axis isṀ w and shown on the left axis is the enhancement of the mass loss rate ∆Ṁ w , which is defined by whereṀ 0 w (= 1.32 × 10 −14 M yr −1 ) denotes the mass-loss rate derived from the case with δv , = 0.0 km s −1 (BsV00). The time variability is also presented by vertical error bars taken from the maximum and minimum values during the time averages. Cases with larger δv , exhibit higher time variability, which is discussed later in Section 4.6.
The blue solid line in Figure 5 is the power-law fit to the time-averaged ∆Ṁ w in a range of δv , ≤ 2.7 km s −1 : The fitting formula indicates that ∆Ṁ w increases almost linearly with δv , until saturating above δv , 2.7 km s −1 . The linear dependence indicates that the increase ofṀ w is slower than the increase of the injected energy flux carried by the longitudinal waves ∝ δv , 2 . This implies that not all of the additional input energy of the longitudinal waves, but a portion of it, is used to enhance the mass loss. One possible reason is that, as δv , increases, a larger fraction of the input longitudinal waves dissipates in the chromosphere due to more efficient shock formation.
Although the mass-loss rate depends on the amplitude of the longitudinal wave in the photosphere, it does not mean that the longitudinal wave is the main driver of the solar wind. As shown in Appendix A, without transverse wave injection (B0V06), the atmosphere is heated only up to a few times 10 5 K and steady outflows do not occur by the acoustic waves from the photosphere. Thus, the interaction between longitudinal and transverse waves is possibly the key to understand the cause of the enhancement of the mass loss. To figure out what caused the increase ofṀ w , we investigate the global energetics of the wind, which is a key to understand the scaling law of mass-loss rate (Cranmer & Saar 2011;Shoda et al. 2020). In particular, we consider the radia-tive energy loss to discuss the energy conservation law from the photosphere to the solar wind (Suzuki et al. 2013).
In the quasi-steady state, the time averaged energy conservation Equation (8) is given by where are the surface-integrated kinetic energy flux, enthalpy flux, Poynting flux, conductive flux, and gravitational potential-energy flux, respectively. We note thatṀ w in Equation (52) can be assumed to be constant in the quasi-steady state. We define the radiation luminosity L R as follows: where r lch is the radial distance in the lower chromosphere. We set r lch − R = 0.7 Mm (r lch /R = 1.001).
Below r < r lch we assume L R = 0 because the exponential (Newtonian) cooling, which dominates the radiation in r r lch , should yield negligible net radiative loss. Equation (47) is then rewritten in terms of L R as follows.
where L tot is the total surface-integrated energy flux, which is expected to be constant in r in the quasi-steady state. By relating the values of L tot at different radial distances, several analytical relations are derived.
1. Photosphere: Because the kinetic, thermal, and conductive energy fluxes are negligibly small on the nearly static and low-temperature photosphere, the dominant terms in L tot are the Poynting flux and the energy flux of gravitational potential, that is, where v g, = 2GM /R = 617 km s −1 is the escape velocity. We note that L R, = 0 is assumed as described above.
2. Coronal base: Because the mean outflow velocity is small at the coronal base, L K and L E are negligible, and thus, L tot is approximated by where the subscript "cb" denotes the value at the coronal base, which we set r cb /R = 1.03. We have confirmed that the conductive luminosity is small at the coronal base because the temperature gradient is already shallow. Therefore, we can safely simplify Equation (56) to 3. Distant solar wind (outer boundary): Because the kinetic energy flux dominates the enthalpy, Poynting, and conductive fluxes in the super-Alfvénic region, L tot,out is approximated by where the subscript "out" denotes the value at the outer boundary (r = r out ). Because the radiative loss above the coronal base is generally negligible, we use L R,out ≈ L R,cb (but see discussion below).
The energy conservation between the photosphere and the coronal base (Eq.s. (55) and (57)) yields where we approximate L G, ≈ L G,cb . Cian stars and red pentagons in Figure 6 respectively denote L R,cb and L A,cb normalized by L A, . Equation (59) is satisfied if the sum of these two components is 100 % in Figure 6. As one can see, however, this conservation is not perfectly fulfilled, possibly because of the treatment of the radiative cooling in the low chromosphere. As described previously, L R excludes the contribution of the radiation cooling below r < r lch . By this assumption, we should underestimate L R , leading to L A, > L A,cb + L R,cb .
Although we have to bear in mind that L R,cb could be underestimated, an increasing trend of L R,cb for δv , is physically plausible; the density in the chromosphere and the low corona is higher for larger δv , (Figure 4), which yields larger radiative cooling. As a result, L A,cb /L A, does not monotonically increase with δv , but eventually saturates for δv , 2 km s −1 (Figure 6) because a large portion of the input Alfvénic Poynting flux is already lost via radiation below the coronal base.
Next, the energy conservation between the coronal base and the outer boundary (Eq.s (57) and (58)) yields The green open circles and black triangles in Figure 6 respectively represent L G,cb and L K,out normalized by L A, . L A,cb , L G,cb and L K,out in Figure 6 exhibit a similar trend on δv , ; they increase in δv , 2 km s −1 and saturate for δv , 2 km s −1 . Figure 6 also shows that Eq.(60) is reasonably satisfied, that is, The wind velocity can be well approximated by the escape velocity, v r,out ≈ v g, = 2GM /R , which is also a reasonable assumption in our numerical results (Table 1). Then, by using L K,out ≈ L G,cb ≈Ṁ w v 2 g, /2, we can rewrite Equation (60) as follows: as already found in Cranmer & Saar (2011). The comparison ofṀ w (blue stars) to L A,cb /v 2 g, (orange circles) in Figure 5 confirms that Equation (61) explains the obtained mass loss rate quite well particularly in the small δv , 1 km s −1 regime. In other words,Ṁ w is primarily determined by the Alfvénic Poynting flux at the coronal base. In contrast, Equation (61) slightly overestimates the obtainedṀ w in the larger δv , 1 km s −1 cases because the radiative cooling above r > r cb is not negligible in Equation (58) owing to the larger density in the corona (Figure 4). However, even in these cases with large δv , , Equation (61) still gives a reasonable estimate ofṀ w .
In summary, the increase and saturation ofṀ w on δv , directly reflect the trend of the Alfvénic Poynting flux at the coronal base. The saturation can be interpreted by the excess of the radiative cooling discussed previously. On the other hand, in order to understand the increase of L A,cb in δv , ≤ 2.7 km s −1 , we need to further examine detailed properties of waves below the coronal base, which is presented in the following subsection.

Dissipation and Mode Conversion of Waves
To understand the dependence of L A,cb on δv , in Figure 6, we examine the propagation and dissipation of transverse waves (≈ Alfvén waves) from the chromosphere to the low corona. Shoda et al. (2020) introduced an equation that describes the evolution of Alfvén waves: ∂ ∂t (62) where ↔⊥ and Q turb indicate the mode conversion from transverse waves to longitudinal waves and the energy loss by turbulent cascade, respectively. These are explicitly written as We note that the first term of Equation (63) denotes the nonlinear excitation of longitudinal perturbation from the magnetic fluctuation associated with transverse waves (Hollweg 1971;Kudoh & Shibata 1999;Suzuki & Inutsuka 2005;Matsumoto & Shibata 2010). Using ε r↔⊥ and Q turb , we define the the energy loss rates via the turbulent dissipation and the mode conversion as and ∆L A,mc (r) = r R dr4πr 2 f op ε ↔⊥ .
(66) Figure 7 shows properties of the damping of Alfvénic waves in the chromosphere; Panel (a) presents the radial profile of L A ; Panels (b) and (c) present ∆L A,turb (energy loss by turbulence) and ∆L A,mc (energy loss by mode conversion), respectively. We note that the net loss of Alfvénic waves, L A, − L A , is not always equal to the sum of these energy losses, possibly because of numerical dissipation.
As shown in Figure 7, the mode conversion rate and the turbulent loss rate behave differently. A general trend is that ∆L A,turb increases and ∆L A,mc decreases with increasing δv , . As δv , increases, the mode conversion from transverse (Alfvénic) waves to longitudinal (acoustic) waves is suppressed; instead, the "inverse conversion" (longitudinal-to-transverse wave energy transfer), ∆L A,mc < 0, takes place in the case with large δv , (red dotted; BsV18). This is because transverse waves are excited at the region with plasma β ≈ 1 in the chromosphere by the mode conversion from the large-amplitude longitudinal waves injected from the photosphere (Cally 2006;Schunker & Cally 2006;Cally & Goossens 2008). We here note that the conversion from the longitudinal mode to the transverse mode occurs even in the simple 1D geometry because the direction of magnetic field, (B rr + B ⊥ )/|B|, is not parallel with the direction of the wave propagation that is strictly along the r direction, wherer is the radial unit vector.
As a result of the excitation of transverse wave from longitudinal wave, L A increases near the surface (Figure  7(a)). The amplitude of the excited transverse waves increases with δv , , which raises the Alfvénic Poynting flux at the coronal base ( Figure 6). As a consequence of increased energy injection to the corona (increased L A,cb ), the coronal heating rate increases, which leads to larger coronal density (see discussion in Section 2.5). In fact, as shown in Figure 4, the density at the coronal base (where T = 5 × 10 5 K) is higher for larger δv , , even though the coronal base is located at a higher altitude.
Another interesting point is that the velocity of the wind is insensitive to the value of δv , (bottom panel of Figure 3); the increase ofṀ w (∝ ρv r ) is solely by the increase in the density. According to the standard model of the solar/stellar winds (Hansteen & Leer 1995;Lamers & Cassinelli 1999), the additional heating and momentum inputs in the subsonic region (v r < a) of a wind raise the mass-loss rate with negligible effects on the terminal velocity, while those in the supersonic region (v r > a) do not affect the mass-loss rate but result in the higher terminal velocity. Based on this background, to understand the behavior of wind speed with respect to δv , , we examine the radial distribution of the energy transfer rate from the Alfvénic wave to the gas from the corona to the wind. Specifically, we calculate the loss rate of the Alfvénic Poynting flux per unit mass, defined as Because of energy conservation, ζ A corresponds to the energy (= heating + work) transfer rate from the Alfvénic Poynting flux to the plasma. Figure 8 presents ζ A of the cases with δv , = 0.6 km s −1 (BsV06; green solid) and 1.8 km s −1 (BsV18; red dashed) normalized by ζ A of δv , = 0 (BsV00). The increase in δv , promotes the energy input in the subsonic region (r 2 − 3R ) but does not affect (or even reduces) the energy input in the supersonic region (r 2 − 3R ). In other words, the vertical oscillation on the photosphere affects only the subsonic region. For this reason, an addition of δv , does not affect the wind velocity but only increases the mass-loss rate.  Figure 9. Comparison of the radial profiles of time averaged βr (Equation (68)). Thin and thick lines correspond to the cases of Br, = 1300 G (BsVyy) and 325 G (BwVyy), respectively. Dashed and solid lines show the cases with δv , = 0 (BxV00) and 0.6 km s −1 (BxV06), respectively. Diamonds show the location of the top of the chromosphere at T = 2 × 10 4 K for each case.
We have shown that the nonlinear mode conversion between transverse waves and longitudinal waves in the chromosphere is the key to determine the wind properties when both transverse and longitudinal perturba-tions are input at the photosphere. The mode conversion rate sensitively depends on plasma β = 8πp/B 2 with peaked at β ≈ 1 (Hollweg 1982;Spruit & Bogdan 1992). Therefore, it is expected that the magnetic field strength in the chromosphere plays an essential role in determining the Alfvénic Poynting flux that enters the corona. Here, we perform simulations in a flux tube with weaker magnetic field from the photosphere to the chromosphere but with the same field strength above the corona ("BwVyy" in Table 1). Figure 9 compares the radial profiles of the time averaged plasma beta values of four cases, BsV00, BsV06, BwV00, and BwV06, that are evaluated from the only radial component of the magnetic field, We note that, although β r ≥ β, the difference is small because |δB ⊥ | < |B r |.
where v A = B r / √ 4πρ is the Alfvén velocity along the r direction. We note that L A = L + A + L − A . Let us begin with the comparison of the cases without δv , , BwV00 (light blue thick dashed lines) and BsV00 (deep blue thin dashed lines). Figure 10 (a) shows that L A of the weak field case (BwV00) is larger than L A of the standard field case (BsV00) near the photosphere. However, the former declines more rapidly in the chromosphere to be smaller than the latter above the upper chromosphere of r/R − 1 > 2 × 10 −3 . As a result, the mass loss rateṀ w of BwV00 is slightly smaller thanṀ w of BsV00 (Table 1), which is consistent with Equation (61). The rapid damping of the Alfvén waves is mainly because of more efficient turbulent dissipation (Figure 10(c)). Utilizing Equation (1), we can rewrite the correlation length (Equation (26)) as follows: where we are adopting the same λ ⊥, = 150 km in both cases (Equation (27)). Since the flux-tube expansion of the weak field case is smaller in our setup, f op /f op is smaller, which enhances the turbulent dissipation. The rapid turbulent damping in the chromoshere reduces the amplitude of the outgoing Alfvén waves in the upper region. The reflected Alfvén waves downward to the photosphere are also suppressed to give the smaller ratio of L − A /L + A near the photosphere of the weak field case (Figure 10(b)). Therefore, the net outgoing Poynting flux, L A (= L + A + L − A ), is larger there (Figure 10(a)). The comparison of BwV06 to BwV00 indicates thaṫ M w increases more than twice by the additional input of the longitudinal perturbation of δv , = 0.6 km s −1 in the weak field condition ( Table 1). The enhancement factor ofṀ w is considerably larger than the value (≈ 1.5 times) obtained in the standard field condition. This is because in the weak field case of BsV06 (orange thick lines) larger Alfvénic Poynting flux reaches the coronal base (Figure 10(a)) by the generation of transverse waves through the mode conversion ( Figure 10(d)) in spite of the higher turbulent loss (Figure 10(c)). Figure 9 shows that β r of this case decreases with height and crosses unity in the chromosphere, which induces the efficient longitudinal-to-transverse mode conversion (Figure 10(d)) as shown by Cally (2006) and Cally & Goossens (2008). In contrast, β r stays < 1 in the standard field case of BsV06 (Figure 9). As a result, ∆L A,mc remains positive (Figure 10 (d)), namely the excitation of transverse waves by the mode conversion is negligible. We can conclude that, in addition to the longitudinal fluctuation in the photosphere (Section 3), the magnetic field strength in the chromosphere is also an essential factor to determine the global wind properties through nonlinear processes of MHD waves.

Limitation of the 1D Geometry
We have simulated the propagation, dissipation, and mode conversion of MHD waves in 1D super-radially open flux tubes. While we took the phenomenological approach to the Alfvén-wave turbulence (Section 2.3) to consider the 3D effect, multi-dimensional effects are also important in other wave processes (e.g. Hasan & van Ballegooijen 2008;Matsumoto & Suzuki 2012;Iijima & Yokoyama 2017;Matsumoto 2021). The nonlinear mode conversion, which is a key process in the present paper, is probably one of those that have to take into account 3D effects because the conversion rate increases with the attacking angle between the direction of a magnetic field line and a wave-number vector (Schunker & Cally 2006). Since the attacking angle tends to be restricted to a small value in the 1D treatment, the amount of the generated transverse waves by the mode conversion may be underestimated in our simulations.
Although we have only considered shear Alfvén waves, torsional Alfvén waves are also expected to be excited (Kudoh & Shibata 1999). The nonlinear steepening of the torsional mode is slower than that of the shear mode (Vasheghani Farahani et al. 2012). Therefore, if we considered torsional waves in addition to shear waves, the dissipation rate of the transverse waves would be slower, which may affect the global wind properties.

Missing physics in the chromosphere
We described that the radiative cooling in the chromosphere governs the saturation of the Alfvénic Poynting flux that reaches the coronal base in the cases of the large δv , (Figure 6 and Sections 3.2 & 3.3). The local thermodynamical equilibrium is not strictly satisfied in the chromosphere, and the radiative cooling is governed by multiple bound-bound transitions (Carlsson & Leenaarts 2012). In addition, the radiative loss also affects the propagation of compressional waves such as acoustic waves (e.g., Bogdan et al. 1996). Ideally, detailed radiative transfer has to be solved to accurately handle these complicated processes, although we have taken the approximated prescription to consider them phenomenologically (Section 3.2). A more accurate treatment (e.g., Hansteen et al. 2015;Iijima & Yokoyama 2017) might modify the radiative loss rate, which we plan to tackle in our future works.
The gas in the chromosphere is partially ionized plasma. The relative motion between charged particles and neutral particles, which is called ambipolar diffusion, promotes damping of transverse waves and heating the gas (Khodachenko et al. 2004;Khomenko & Collados 2012). However, in the current condition of the solar chromosphere, the ambipolar diffusion does not give a large impact on the low-frequency Alfvén waves with < 10 −2 Hz considered in this paper (Equation (39); Arber et al. 2016), although it may affect higherfrequency Alfvén waves and rapid dynamical phenomena (e.g., Singh et al. 2011). Another interesting aspect is that magnetic tension, which is induced by ambipolar diffusion, can be an additional generation mechanism of transverse waves in the chromosphere (Martínez-Sykora et al. 2017). In future, the effect of partial ionization should be investigated also in the context of the solar/stellar wind studies. . Comparison of the radial profiles of the rootmean-squared density fluctuations normalized by the average density δρ / ρ = ρ 2 − ρ 2 ave / ρ of the three cases, BsV00 (dashed blue), BsV06 (solid green), and BsV18 (dotted red). Diamonds show the location of the top of the chromosphere at T = 2 × 10 4 K for each case.

Density fluctuation
from more efficient longitudinal-to-transverse mode conversion in the chromosphere (Figure 7) in addition to more rapid dissipation of longitudinal wave by shock formation.
Above that, all the presented three cases exhibit a first peak of δρ / ρ around r/R − 1 ∼ (2 − 5) × 10 −2 from the transition region to the low corona. This peak reflects time-variable spicule activities (see Section 4.6); density fluctuations are excited by the nonlinear mode conversion from transverse waves to longitudinal waves via the gradient of the magnetic pressure associated with the Alfvénic waves (Hollweg 1971;Kudoh & Shibata 1999;Matsumoto & Shibata 2010). As explained in Figure 7, the upward Alfvénic Poynting flux is larger for larger longitudinal-wave injection at the photosphere even though the same amplitude of transverse fluctuations is excited. Therefore, taller spicules are generated and the first peak is located at a higher altitude for larger δv , .
Although the location of the first peak depends on δv , , the radial profiles of δρ / ρ converge to a similar level above r/R 1 almost independently from δv , . A gentle second peak of δρ / ρ is formed around r/R ∼ 5 by the parametric decay instability of Alfvénic waves (Terasawa et al. 1986;Tenerani et al. 2017;Suzuki & Inutsuka 2006;Shoda et al. 2018b;Réville et al. 2018).
Since the density fluctuation in the corona and solar wind is observable by remote sensing, our model could  Figure 11. Red and blue circles respectively represent the observed amplitudes of z + ⊥ and z − ⊥ from PSP (Chen et al. 2020). Diamonds show the location of the top of the chromosphere at T = 2 × 10 4 K for each case.
be constrained by the detailed comparison with observation (Miyamoto et al. 2014;Hahn et al. 2018;Krupar et al. 2020). For example, it is reported that the relative density fluctuation in the coronal base is as large as 10 % or larger (Hahn et al. 2018;Krupar et al. 2020), which possibly indicates the non-negligible fraction of compressional waves present in the coronal base.

Alfvénicity of the solar wind
In most cases, the simulated solar wind is fast, in that the termination velocity mostly exceeds 500 km s −1 . The fast solar wind is known to be Alfvénic, that is, the outward Elsässer variable is much larger than the inward Elsässer variable in the fast streams. Here the Alfvénic nature of the simulated solar wind is discussed. Figure 12 compares the numerical results of the time averaged outgoing and incoming Elsässer variables z + and z − to observations at r ≈ 40R by the Parker Solar Probe (hereafter PSP) (Chen et al. 2020). The obtained z ± from these three cases are consistent with the observed z + and z − that show large scatters. Both z + and z − of our numerical results are larger for larger δv , in the coronal regions of 10 −2 < r/R − 1 < 1. However, larger z ± yields larger turbulent loss as shown in Figure 7 , which suppresses the increase of z ± . As a result, almost the same maximum amplitudes of z ± are obtained for the different δv , cases in a self-regulated manner, sim-ilarly to the density perturbations ( Figure 11). In the solar wind region, r/R 10, z ± is smaller for larger δv , because the density is higher (Figure 4) . Figure 13 shows time versus radial distance diagrams of the mass density in the low atmosphere. The left and right columns respectively present the cases with the standard (BsVyy) and weak (BwVyy) magnetic field in the chromosphere. The top, middle, and bottom rows correspond to δv , = 0, 0.6, and 1.8 km s −1 , respectively. The yellow lines represent the positions of T = 2 × 10 4 K, which correspond to the bottom of the transition region.

Time Variability
One can see that the transition region moves up and down in all the cases, which should be observed as spicules (e.g., Beckers 1972;De Pontieu et al. 2007;Shoji et al. 2010;Okamoto & De Pontieu 2011;Yoshida et al. 2019;Tei et al. 2020). The velocity of the upward motions can be derived to be an order of ∼ 10 km s −1 from the slope of yellow lines, which roughly coincides with the sound speed and accordingly the propagation speed of slow MHD waves. Namely the gas in the upper chromosphere is lifted up by longitudinal slow-mode waves that are generated from transverse waves through the mode conversion in the upper chromosphere (Hollweg 1982;Suematsu et al. 1982;Matsumoto & Shibata 2010;Sakaue & Shibata 2021).
In the cases without longitudinal perturbations (top panels of Figure 13), the stronger field in the chromosphere (top-left) gives the higher and more dynamical transition region because of the larger Alfvénic Poynting flux ( Figure 10 and Section 4.1). By adding the longitudinal perturbations at the photosphere, the transition region behaves more abruptly and the average height increases. The comparison of the middle right panel to the top right panel indicates that the activity of the transition region is drastically enhanced in the weak field case by the addition of δv , = 0.6 km s −1 . This is because in the weak field condition transverse waves are generated more effectively from acoustic waves around β r = 1 (blue line) in the chromosphere as discussed in Section 4.1.
The bottom panels of Figure 13 exhibit a dynamical behavior of transition regions with chromospheric gas being violently uplifted to higher altitudes by the injection of the large-amplitude vertical perturbation with δv , = 1.8 km s −1 . Multiple yellow lines are frequently plotted at single time slices. This indicates that cooler gas with T < 2 × 10 4 K is transiently distributed above hotter gas with T > 2 × 10 4 K. We investigated how the properties of solar winds are affected by the p-mode like longitudinal perturbation at the photosphere. We performed 1D simulations from the photosphere to beyond several tens of solar radii for Alfvén-wave-driven winds in the wide range of the amplitude of the vertical perturbation of 0.0 km s −1 ≤ δv , ≤ 3.0 km s −1 in super-radially open magnetic flux tubes.

SUMMARY
The coronal temperature and wind velocity are not significantly affected by the additional input of the lon-gitudinal perturbation (Figure 3). However, higher coronal density is obtained for larger δv , (Figure 4), and accordingly, the mass-loss rate increases with δv , by up to ≈ 4 times ( Figure 5) because larger Alfvénic Poynting flux enters the corona so as to drive denser outflows as a result of more efficient chromospheric evaporation. The p-mode like vertical oscillation excites acoustic waves, a part of which is converted to the transverse waves by the mode conversion in the chromosphere (Figure 7). These transverse waves contribute to the upgoing Alfvénic Poynting flux, in addition to the Alfvén waves that come from the photosphere. This result confirms the observationally inferred link between p-mode oscillations and Alfvénic waves in the solar corona (Morton et al. 2019).
Cases with larger δv , exhibit higher time variability and larger density perturbations in the low corona. The mass loss rate saturates when δv , 2.5 km s −1 , because an increase of δv , no longer leads to the excitation of transverse waves by the mode conversion but instead is compensated by the radiative loss by the direct shock dissipation of acoustic waves in the chromosphere.
Simulations with weaker field strength in the low atomosphere show that the magnetic field in the chromosphere controls the mode conversion between longitudinal and transverse modes. In the cases that include a region with plasma β ≈ 1 in the middle chromosphere, the mode conversion effectively generates transverse waves there even for a moderate amplitude of δv , = 0.6 km s −1 .
We conclude that p-mode oscillations at the photosphere play an important role in enhancing Alfvénic Poynting flux over the corona of the Sun and solar-type stars.
Numerical simulations in this work were partly carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. M.S. is supported by a Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows and by the NINS program for cross-disciplinary study (grant Nos. 01321802 and 01311904) on Turbulence, Transport, and Heating Dynamics in Laboratory and Solar/ Astrophysical Plasmas: "SoLaBo-X." T.K.S. is supported in part by Grants-in-Aid for Scientific Research from the MEXT/JSPS of Japan, 17H01105 and 21H00033 and by Program for Promoting Research on the Supercomputer Fugaku by the RIKEN Center for Computational Science (Toward a unified view of the universe: from large-scale structures to planets, grant 20351188-PI J. Makino) from the MEXT of Japan.

A. ACOUSTIC WAVE-DRIVEN WIND
We examine the properties of the atmosphere of the case with the only longitudinal fluctuation (B0V06) to see if the corona and solar wind are formed solely by acoustic waves. In order to avoid the initial infall of material from the upper layer, we set lower initial density (ρ w,0 = 10 −25 g cm −3 in Eq. (36)) than that of the other cases. Figure 14 presents the radial profile of the atmosphere averaged from t = 3500 min to t = 5000 min. The acoustic waves that travel upward from the photosphere rapidly dissipate at low altitudes, r/R − 1 0.1. Although the atmosphere is heated up by wave dissipation of longitudinal waves, the temperature of the "corona" remains low ( 2 × 10 5 K, Figure 14 (a)) As a result, the gas in the upper atmosphere does not stream out. Instead, it falls down to the surface, which is seen as negative mass loss rate,Ṁ w < 0, in r/R − 1 > 1 (Figure 14 (c)). The accretion reduces (raises) the density in the outer (inner) region of r/R − 1 <(>)0.5 from the initial value (Figure 14 (b)).
The accretion occurs partially because the initial density in the upper region is still higher than the hydrostatic density with T 10 5 K. We note, however, that the initial density is much lower by 6 -7 orders of magnitude than the observed density in the solar wind. This simulation demonstrates that even such low-density gas cannot be driven outward only by the acoustic waves. We thus conclude, through a direct numerical demonstration, that the solar coronal heating and the solar wind driving cannot be accomplished only by the sound waves from the photosphere.