Acceleration of polytropic solar wind: Parker Solar Probe observation and one-dimensional model

As an introduction, we briefly review the 1D isothermal solar wind model with adiabatic index $\gamma=1$, which was first analyzed by Ref. [Parker, 1958]. Considering a flux tube with cross section area $A(r)$, the mass flux conservation gives

where $\rho,V$ are the density and radial solar wind speed. The momentum equation is

with $C=\sqrt{k_{B}T/m_{p}}$ being the sound speed where $k_{B}$ is the Boltzmann constant, $m_{p}$ is the proton mass and $T$ is the sum of proton and electron temperatures (we have assumed the two temperatures are identical). In a spherically expanding solar wind ($A(r)=r^{2}$), the momentum equation can be re-arranged in the following form

with the mass-conservation relation. The l.h.s. term indicates that if we require a solution of $V(r)$ that starts from a very small value at the inner boundary and increases to larger than the sound speed, the equation has a singular point

at which $V(r)=C$. We can integrate equation (3) from this critical (sound) point to any radial location $r$ and get the Bernoulli’s equation

In FIG. 2, we show the solution of the isothermal wind model (equation (5)) for different coronal temperatures by solid curves. We note that there are actually two branches of solutions that pass through the critical point, but the other solution has $dV/dr<0$, i.e. wind speed decreasing with distance and passing through the sound speed from above, thus it is not the solar wind solution we are interested in. However, it is worth noting that since the momentum equation is invariant under the transformation $V\rightarrow-V$, the decelerating solution represents an accretion flowBondi (1952); Velli (1994) which is likely related to the star formationMcKee and Ostriker (2007). But as accretion is not important in the solar system, we do not analyze the decelerating solution in detail.

In this section, we show how to self-consistently find the transonic solutions, i.e. solutions that pass through the sound point smoothly, with $\gamma>1$. The polytropic case is more complicated than the isothermal case in that the critical point cannot be explicitly determined. A self-consistent analysis of the polytropic wind solution can be found in the textbook Ref. [Meyer-Vernet, 2007], while the approach presented in this section is different from that in Ref. [Meyer-Vernet, 2007].

The mass-conservation relation gives

where $\rho_{0}$, $V_{0}$, and $A_{0}$ are the density, velocity, and cross section area at the inner boundary $r_{0}$. The polytropic relation gives $P\rho^{-\gamma}=Const$ or equivalently $T\rho^{-(\gamma-1)}=Const$. The momentum equation is

where $f(r)$ is the external force per unit mass. In this section, we set $f(r)=0$. Using the polytropic relations, the pressure gradient term can be written as

Plug it into equation (7) and integrate from $r_{0}$ to $r$, we get

The subscript “0” indicates quantities at the inner boundary $r_{0}$. The above equation includes two integral constants: $V_{0}$ and $T_{0}$. The temperature $T_{0}$ should be a given inner boundary condition but $V_{0}$ is not a free parameter and should be determined by the constraint that the solution $V(r)$ passes through the critical point smoothly. So, let’s revisit the momentum equation (equation (7)). We can write the pressure gradient term as

where

is the square of local sound speed and $C_{0}^{2}=\gamma k_{B}T_{0}/m_{p}$ is the square of sound speed at $r_{0}$. Plug the above relation into equation (7), we get

One can easily show that for $\gamma=1$ and $A=r^{2}$, the above equation reduces to equation (3). Equation (12) implies that the location of the critical point $r_{c}$ and the velocity at the critical point $V_{c}$ satisfy

or

One more equation is needed to relate $V_{0}$ that appears in $C_{s}$ (equation (11)) with $r_{c}$ and $V_{c}$. Writing the Bernoulli’s equation (equation (9)) at $r=r_{c}$ and using the relation $T_{c}=T_{0}\left(A_{0}V_{0}/A_{c}V_{c}\right)^{\gamma-1}$ to substitute $T_{c}$, we get

Equations (13a), (14), and (15) form a closed equation set:

from which we can solve $V_{0}$, $V_{c}$, and $r_{c}$ simultaneously. Here we note that this three-equation model can be easily extended to a multi-temperature fluid, e.g. a wind where the ion and electron have different temperatures and polytropic indicesDakeyo et al. (2022). The only modifications are the definition of the sound speed and the pressure gradient term in equation (16c) as we need to sum up the contributions from all species. For simplicity, we assume a single-fluid solar wind throughout this study.

We consider a radially-expanding solar wind ($A(r)=r^{2}$). Then from equation (16b), we get

where we have defined the normalized radius $s=r/r_{0}$ and

which is one half of the escape velocity. We note that $r_{0}$ does not necessarily equal to the solar radius $r_{s}$. As will be discussed in Section III.3, $r_{0}$ could be the outer radius of the isothermal layer formed due to large heating or thermal conduction in the lower corona Dakeyo et al. (2022). Plug equation (18) in equation (16a), we get

With equations (17) and (19), we can eliminate $V_{0}$ and $V_{c}$ in equation (16c) and get a single equation for $s_{c}$:

This is an algebraic equation, so we can numerically find all the roots $s_{c}(\gamma,C_{0}/C_{g})$. An alternative form of the equation can be found in Ref. [Rèville, 2016]. In Appendix A, we describe in detail how to determine the number of roots of equation (20) with the knowledge of $C_{0}/C_{g}$ and $\gamma$. In FIG. 3 we show the phase diagram of the number of roots of equation (20) on the $(C_{0}/C_{g})^{2}$-$\gamma$ plane. Dark red region has two roots, orange region has one root, and the white region has no root. The yellow line is $(C_{0}/C_{g})^{2}=2(\gamma-1)$. We have numerically verified, by comparing $V_{c}$ and $V_{0}$, that if there is one root, only the decelerating solution exists, i.e. solar wind solution can only be found in the dark red region, bounded by the two lines $C_{0}/C_{g}=1$ and $(C_{0}/C_{g})^{2}>2(\gamma-1)$. Naturally, $\gamma<3/2$ must also be satisfied as the two lines cross at $\gamma=3/2$. Early study Parker (1960) using asymptotic analysis points out that physically meaningful solar wind solution exists only with $\gamma<3/2$ and $1\gg(C_{0}/C_{g})^{2}>(\gamma-1)/2\gamma$. Our result gives a more precise constraint on $(C_{0}/C_{g})^{2}$. The cyan lines in FIG. 3 mark the realistic parameters of the Sun with $r_{0}=r_{s}$ and varying temperature at the inner boundary $T_{0}$. One can read that for $T_{0}=2$MK solar wind does not exist for $\gamma\gtrsim 1.1$, and even for a very hot corona with $T_{0}=5$MK, the solar wind does not exist for $\gamma\gtrsim 1.3$.

Figure 4 shows $s_{c}$ as a function of $T_{0}$ with fixed $\gamma$ (panel (a)) and $s_{c}$ as a function of $\gamma$ with fixed $T_{0}$ (panel (b)). In this plot, we have set $r_{0}=r_{s}$. Blue curves show the accelerating (solar wind) solution and orange curves show the decelerating solution. The black curve in panel (a) and the red dots in panel (b) correspond to an isothermal plasma ($\gamma=1$). Different from the isothermal case where the accelerating and decelerating solutions pass through the same critical point, the two branches of solutions with $\gamma>1$ have different critical radii because the temperature depends on the wind speed profile. As $T_{0}$ decreases and $\gamma$ increases, the critical radii for both the two branches of solutions increase, and $r_{c}$ of the accelerating solution increases much faster than that of the decelerating solution. Consistent with FIG. 3, one can read from panel (b) of FIG. 4 that $r_{c}$ for the accelerating solution diverges to extremely large values (actually infinity) as $\gamma$ approaches certain critical values, corresponding to the cross points between the cyan and yellow lines in FIG. 3.

After $s_{c}$ is solved, we can easily calculate $V_{c}$ from equation (16b) and then integrate the momentum equation from $r_{c}$:

to acquire the profile of $V(r)$. We note that, the above equation gives two branches of $V(r)$ starting from one critical point, but only one branch satisfies the given inner boundary condition. This is similar to the early work on polytropic accretion flowBondi (1952), which shows that there are two branches of solutions that cross the critical point but only one of them satisfies the boundary condition at infinity. In Figure 2, we plot $V(r)$ for $\gamma=1.09$ and varying $T_{0}$ in dashed-dotted curves. Compared with the isothermal case, the solar wind speed drops significantly even though $\gamma$ is only slightly larger than 1. In Figure 5, we show the profiles of the solved solar wind speed (top) and the corresponding temperature (bottom) for a fixed inner boundary temperature $T_{0}=3$MK and different values of $\gamma$. We see that the wind speed decreases very rapidly with the increasing $\gamma$.

In conclusion, results from this section show that a simple polytropic wind model cannot explain the observed solar wind speed with the polytropic index deduced from in-situ data as shown in Section II. In the following two sections, we will discuss two mechanisms that possibly contribute to solve the problem.

In the lower corona, heating mechanisms such as magnetic reconnection and high thermal conduction may prevent the plasma temperature from decaying much. Hence, one way to overcome the difficulty that solar wind solution does not exist with large polytropic index is to assume there is an isothermal zone Dakeyo et al. (2022) with radius $r_{iso}>r_{s}$ such that $T(r\leq r_{iso})\equiv T_{0}$. In the recent study by Ref. [Dakeyo et al., 2022], this “iso-poly” solar wind model is solved numerically with quite thick isothermal layers such that the critical point always falls inside the isothermal layer. In this section, we discuss the model in more details and analyze the case when the critical point is outside the isothermal layer. Especially, we point out that, due to the discontinuity in $\gamma$ at the interface between the isothermal and polytropic layers, the values taken by $r_{iso}$ present a radial interval (around the critical point), for which iso-poly solar wind solutions do not exist.

For the polytropic layer, if the critical point is still inside the layer ($r_{c}>r_{iso}$), the procedure described in the prior section to solve equation (16) remains exactly the same. We only need to set the inner boundary at $r_{0}=r_{iso}$ and re-define $C_{g}$ using equation (18). In FIG. 6, we show the phase diagram of the solutions to equation (16) on the $r_{0}/r_{s}-\gamma$ plane, with a base temperature $T_{0}=3$MK. The plot is similar to FIG. 3 but here the number of solutions is determined by numerically searching roots of equation (20) in the range $r\in[1,20000]r_{0}$ instead of using the method described in Appendix A. Thus it also serves as a verification of the method described in Appendix. In this plot, the yellow curve marks $C_{0}/C_{g}=1$ and the cyan curve marks $(C_{0}/C_{g})^{2}=2(\gamma-1)$. The horizontal dashed line marks the critical radius for the isothermal wind $r_{c,iso}$. Note that the yellow curve intersects with the horizontal line at $\gamma=1$ (equation (4)).

If $(r_{iso},\gamma)$ falls in the red region, i.e. the iso-poly critical point is in the polytropic layer. We can first calculate the solar wind solution in the polytropic layer following the procedure in Section III.2. We will also acquire the solar wind speed at the inner boundary $r_{iso}$. As the wind speed must be continuous, $V(r_{iso})$ serves as the upper boundary condition for the isothermal layer. Thus, we can then integrate the isothermal momentum equation back from $r_{iso}$ to $r_{s}$. In panel (a) of FIG. 7, we show profiles of solar wind speed for $T_{0}=3$MK, $\gamma=1.1$ and varying radius of the isothermal layer $r_{0}$ ($r_{iso}$). Black curve is $r_{0}=r_{s}$, i.e. no isothermal layer, and green curve is $r_{0}=2r_{s}$, a case that the critical point falls in the polytropic layer. The circles mark the locations of $r_{0}$. One can see that raising the height of the isothermal layer indeed increases the solar wind speed. Besides, FIG. 6 shows that changing $r_{iso}$ also changes the range of $\gamma$ in which a solar wind solution can be found, though the maximum $\gamma$ value allowed is always $3/2$ with $r_{c}$ in the polytropic layer.

If we set $r_{iso}>r_{c,iso}$ (in the green region on FIG. 6), no transonic solution exists in the polytropic layer and the critical point, if exists, falls inside the isothermal layer such that $r_{c}=r_{c,iso}$. In this case, we can first determine the solution in the isothermal layer, which is simply the classic Parker’s model (Section III.1). The solution then gives $V(r_{iso})$, which serves as the inner boundary condition for the polytropic layer. In panel (a) of FIG. 7, the orange dashed curve is $r_{0}=5.78r_{s}=1.5r_{c,iso}$. We get a higher wind speed than the cases $r_{0}=2r_{s}$ and $r_{0}=r_{s}$. However, one limitation appears: the difference between the isothermal sound speed and the polytropic one, does not ensure that the wind speed continues to increase beyond $r_{0}$. For instance, let us consider the following scenario: $r_{0}$ is only slightly larger than (or equal to) $r_{c,iso}$ so that the isothermal layer gives $V(r_{0})\gtrsim\sqrt{k_{B}T/m_{p}}=C$. However, for the polytropic solar wind, we have $C_{0}=\sqrt{\gamma k_{B}T/m_{p}}>\sqrt{k_{B}T/m_{p}}$, which may lead to $V(r_{0})<C_{0}$ in the polytropic layer, i.e. the supersonic solar wind suddenly becomes subsonic across $r=r_{0}$ for this case. In other words, since the wind speed $V(r)$ does not overcome $C_{0}$ before going out of the isothermal layer, there is no possible supersonic expansion in the polytropic layer (no overlap between the green and red regions in FIG. 6), and the flow can not remain supersonic above $r_{0}$. Then one will get either a blow-up solution with $dV/dr\rightarrow\infty$ at some point or a “breeze” with $V(r\rightarrow+\infty)\rightarrow 0$. The blue curve on the panel (a) of FIG. 7 illustrates the case of a "breeze" solution with $r_{0}=3.85r_{s}=r_{c,iso}$, for which beyond $r_{0}$, the accelerating solar wind suddenly starts decelerating and the wind speed tends to zero at large $r$. Another limitation of the iso-poly model appears when we consider the white parameter space (region between the yellow ($C_{0}/C_{g}=1$) curve and the $r_{c,iso}$ line) of FIG. 6. Within this parameter space, the isothermal layer is thinner than $r_{c,iso}$ so that the solar wind cannot be accelerated to a supersonic speed below $r_{iso}$. Meanwhile, the polytropic layer does not have an accelerating transonic solution. Thus, for this parameter region, we can not find a transonic iso-poly solution.

In panel (b) of FIG. 7, we plot iso-poly wind speed profiles for $T_{0}=3$MK, $r_{0}=5.78r_{s}=1.5r_{c,iso}$ and different values of $\gamma$. The black curve represents the isothermal case, and dark to light colors correspond to increasing values of $\gamma$. The red circle marks $r_{iso}$ and the red square marks $r_{c}$(=$r_{c,iso}$). Below $r_{iso}$, the solutions are exactly the same for all the $\gamma$. Indeed, within the isothermal layer, all the iso-poly models have the same initial temperature, and they passe through the same sonic point. Slightly above $r_{iso}$, a larger $\gamma$ value leads to a larger acceleration of the wind on a short radial distance, but it also gives a smaller asymptotic wind speed. We note that, in this case (solutions in the green region of FIG. 6), $\gamma$ can be larger than $3/2$.

In conclusion, the iso-poly model is able to produce transonic wind solutions under either of the following two conditions: $(i)$ $r_{iso}>r_{c,iso}$ and $V(r_{iso})>(C,C_{0})$, i.e. the wind is supersonic and faster than both sound speeds as it enters the polytropic layer. This case corresponds to the orange curve in panel (a) of FIG. 7. $(ii)$ ($\gamma$, $r_{iso}$) falls in the dark red region in FIG. 6, i.e. the wind is subsonic when it enters the polytropic layer and is accelerated to supersonic speeds in the polytropic layer. This case corresponds to the green curve in panel (a) of FIG. 7.

In this section, we discuss the influence of the external force on the polytropic wind model. In the solar wind, Alfvén wave pressure gradient may provide such force as the amplitude of the wave magnetic field can be comparable to the mean magnetic field Velli, Grappin, and Mangeney (1991); Kasper et al. (2019); Bale et al. (2019). Rigorously, we need to model the amplitudes of outward and inward propagating Alfvén waves with two additional equations Verdini et al. (2009); Lionello et al. (2014); Réville et al. (2020); Chandran (2021). But these new equations will greatly complicate the system as they introduce a new critical point, the Alfvén pointHeinemann and Olbert (1980); Velli (1993). Thus, in this section we assume the force $f(r)$ is a known function of $r$, so that the problem can be solved semi-analytically like we did in previous sections. The expression of $f(r)$ is

where $\alpha$ and $\beta$ are two constants, $s=r/r_{0}$ and $f(r_{0})=f_{0}$. Asymptotically, there is

i.e., the external force decays as $r^{-2}$, consistent with a Wentzel-Kramers-Brillouin (WKB) decay of the Alfvén wave Alazraki and Couturier (1971); Belcher (1971); Hollweg (1974b), which predicts an asymptotic $r^{-3/2}$ decay of the magnetic field fluctuations. The integral of $f(r)$ can also be written analytically:

We use one dimensionless parameter $F$ to measure the strength of $f_{0}$ such that $f_{0}=F\times(GM/r_{0}^{2})$. In this study, we fix $\alpha=3$, $\beta=8$, such that $f(r)$ peaks at $r=2.22r_{0}$, and we set $r_{0}=r_{s}$. The radial profile $f(r)$ is plotted in panel (a) of FIG. 8, where the curves with different colors correspond to different values of $F$ and the black dashed line is the gravity force. We note that the choice of the analytic form of $f(r)$, including values of $\alpha$ and $\beta$, is not fully rigorous. The only physics-based requirements are that the wave pressure gradient increases in the lower corona and asymptotically decays as $r^{-2}$. One can adjust the profile of $f(r)$, which will result in different profiles of the solution $V(r)$. But a parametric study of the influence of the shape of $f(r)$ on the solution $V(r)$ is beyond the scope of this study and is unnecessary for drawing the main conclusion of this section.

We start from the isothermal case ($\gamma=1$). With the external force, the critical point at which $V=C=\sqrt{k_{B}T/m_{p}}$, is determined by

Obviously, adding an external force below the original sound point, i.e., the sound point without external force, does not change the location of the sound point nor the profiles of $V(r)$ beyond the sound point. After numerically solving the critical point from equation (25), we can integrate the momentum equation from the critical point, which gives

to recover the profile of $V(r)$. We assume a spherical expansion such that $A(r)=r^{2}$. In panel (b) of FIG. 8, solid curves show the solar wind speed $V(r)$ with $T=3$MK and different $F$. It is clear that a stronger external force lowers the critical point and increases the solar wind speed significantly.

In the polytropic case, the closed equation set with non-zero $f(r)$ is similar to equation (16):

From equation (27b), we get

where we have defined a “modified critical radius” $s_{c}^{\prime}$ for convenience such that

Plug it into equation (27a), we get

Then, plug $V_{c}$ and $V_{0}$ into equation (27c), we acquire the equation for $s_{c}$:

which is solved numerically. Finally, we integrate the montum equation from the critical point and get the profile $V(r)$. In panel (b) of FIG. 8, the dashed curves show $V(r)$ with $T_{0}=3$MK, $\gamma=1.10$ and varying $F$. Similar to the isothermal case, the external force increases the wind speed significantly. In addition, the external force weakens the constraint on the polytropic index $\gamma$. We have verified that a higher $F$ results in a larger range of $\gamma$ in which a solar wind solution can be found. For instance, with $T_{0}=3$MK, the solar wind solution is only possible for $\gamma\lesssim 1.15$ (FIG. 3) without the external force. However, with $F\gtrsim 0.04$, $\gamma$ can be larger than $3/2$, and with $F\gtrsim 0.05$, $\gamma$ can take any value between $1$ and $5/3$. In panel (c), of FIG. 8, we plot $V(r)$ with $T_{0}=3$MK, $F=0.05$ and $\gamma$ varying from $1.0$ to $5/3$. Even with $\gamma=5/3$, the solar wind is still accelerated to over 600 km/s. However, we emphasize that the external force used here is not a self-consistent physics model. The selected values of $F$ lack comparison with observations as the radial profile of the wave pressure in the corona cannot be easily obtained from remote-sensing data. Assuming the wave amplitude at solar surface is $\delta u=30$km/s, a typically value used in Alfvén wave modelsRéville et al. (2020), and the scale of variation of the wave pressure is roughly $L=0.1r_{s}$, one can estimate $F\approx(\delta u)^{2}/L/(GM_{s}/r_{s}^{2})=0.047$. Hence, a choice of $F=0.05$ may be reasonable, but future studies are necessary to quantify the Alfvén wave pressure in the solar corona and solar wind.