Constraints on Solar Wind Density and Velocity Based on Coronal Tomography and Parker Solar Probe Measurements

CRT is a method used to estimate the distribution of coronal electron density. White light coronagraph observations from the CoR2 coronagraph from the Sun–Earth Connection Coronal and Heliospheric Investigation (SECCHI; Howard et al. 2008) suite on board the Solar Terrestrial Relations Observatory (STEREO) spacecraft, are calibrated as outlined in Morgan (2015). A step is then taken to reduce the signals from coronal mass ejections (CMEs) which cause rapid large-scale changes in the density state of the solar corona (Morgan et al. 2012). Finally, the resultant data are passed through an inversion technique based on spherical harmonics outlined in Morgan (2019) to estimate the coronal electron density. This method uses a regularized least-squares regression approach to optimize the tomographic inversion parameters and to best fit the tomographic and observed polarized brightness. The final step narrows the width of structures present in the high-density streamer belt (Morgan & Cook 2020). See Bunting & Morgan (2022) for a more in-depth overview of the CRT method. An example of the tomography density map can be seen in the top panel of Figure 2 for Carrington rotation (CR) 2238 (2020 November). A large archive of tomography maps is available. ³

Zoom In Zoom Out Reset image size

The tomographic densities given by the CRT method can be used to estimate the solar wind velocity, which are then used as an inner boundary condition for solar wind models (Bunting & Morgan 2022, 2023). Bunting & Morgan (2022) found that a simple inverse linear density–velocity conversion model generated an inner boundary condition that, when coupled with a heliospheric model, gave model predictions that agreed well with in situ data at 1 au. However, it was found that this method gave inconsistent results. Bunting & Morgan (2023) found that an inverse exponential relationship generated inner boundary conditions that gave more consistent and accurate predictions of solar wind conditions.

In this work, we use the same exponential model outlined in Bunting & Morgan (2023). The tomographic densities are extracted from the equatorial region of the tomography maps at heliocentric distance 8 R_⊙. These electron densities give an estimate of the solar wind speeds by

$$\begin{eqnarray}&&V=\left({V}_{\max }-{V}_{\min }\right)\ast {e}^{-\alpha * {\rho }_{N}}+{V}_{\min },\end{eqnarray} \tag{ 1 }$$

where V_max and V_min are scaling parameters for the exponential model, α is a unitless exponential factor, and ρ_N is the value of the coronal electron density to be converted into velocity, normalized between the minimum and maximum values of density of all available tomographic time periods (2008–2021), namely 4.40 × 10³ cm⁻³ and 5.64 × 10⁴ cm⁻³.

The time-dependent HUXt model is a lightweight, highly efficient open source heliocentric solar wind model written in Python (Owens et al. 2020; Barnard & Owens 2022). This model ingests solar wind velocities at the inner boundary and propagates solar wind conditions into the heliopshere. The efficiency of this model is increased by assuming an incompressible flow and that pressure gradient, gravitational, and magnetic forces in the solar wind can be negated. This model also solves the solar wind velocities along a single line of latitude, negating non-radial aspects of the ambient solar wind. Despite these assumptions, the model predictions of HUXt over a 40 yr period agree to 6% with full 3D magnetohydrodynamic models (Owens et al. 2020).

Recent advancements have allowed HUXt to perform in three dimensions, as a collection of two-dimensional slices at a range of latitudes (Barnard & Owens 2022). For input, the 3D HUXt requires a "velocity map" which is formed of the velocity state at the inner boundary. The 3D HUXt model becomes extremely useful when extracting predictions at a location that is moving across a range of latitude, such as PSP for example.

HUXt adopts a residual acceleration model, in which a velocity V(r) can be calculated at heliocentric distance (r) by

$$\begin{eqnarray}&&V(r)={V}_{0}\left(1+{\alpha }_{\mathrm{IP}}\left[1-\exp \left(\displaystyle \frac{-\left(r-{r}_{0}\right)}{{r}_{{\rm{H}}}}\right)\right]\right),\end{eqnarray} \tag{ 2 }$$

where V₀ is the initial velocity at the inner boundary distance (r₀), r_H is a scale height and α_IP is the interplanetary residual acceleration parameter. Riley & Lionello (2011) proposed values of r_H and α_IP of 50 R_⊙ and 0.15 that best matched the HelioMAS model between the distances of 30 R_⊙ to 1 au. Within HUXt this acceleration is applied to all inner boundary velocity values independent of magnitude, that is, there is no consideration that slow and fast streams have different accelerations.

Effectively, r_H determines the distance at which the solar wind velocity reaches its near-final velocity. The α_IP parameter determines the residual percentage increase of the velocity relative to the inner boundary. However, it is known that the fast solar wind accelerates more rapidly and reaches its final velocity at distances closer to the Sun, compared to the slower solar wind as seen in Figure 1 and suggested by previous studies (e.g., Breen et al. 2002; Sanchez-Diaz et al. 2016; Morgan & Cook 2020; Halekas et al. 2022, etc.). In the context of the acceleration model, this suggests that the magnitude of r_H and α_IP must be dependent on the magnitude of the solar wind velocity at the inner boundary.

Figure 1 demonstrates the effect of varying the r_H and α_IP parameters on the evolution of the solar wind. The blue, green, and red plots show HUXt radial solutions with different initial velocities at the inner boundary. With increasing initial velocity, r_H is increased and α_IP is decreased in order to approximate the measured distributions. Due to the intermittent and uneven sampling of the solar wind by PSP at various distances, these profiles have been changed manually to fit, approximately, the slowest and fastest winds, and the most probable intermediate velocities. Parameters for slow, intermediate, and fast streams are listed in Table 1. The blue plot approximates the evolution of the slowest solar wind, green represents the intermediate, and red represents the fastest. The slowest solar wind velocity is still accelerating, albeit gradually, at 1 au (∼215 R_⊙) and the fastest solar wind reaches its near-final velocity at ∼150 R_⊙. These plots also approximately follow the distribution of the OMNI data at 215 R_⊙. In the following work, we alter the HUXt acceleration parameters seen in Equation (2) by linearly interpolating between the r_H and α_IP depending on the initial solar wind velocities at the inner boundary as listed in Table 1.

Examples of HUXt results with a varying inner boundary velocity are shown in Figure 3. In this plot, we show the effect of varying the velocity at the inner boundary on the subsequent evolution of the velocity when the r_H and α_IP values are interpolated based on values shown in Table 1. As the initial velocity increases, the acceleration between 8 R_⊙ and 1 au also increases due to the increasing magnitude of the α_IP parameter. Furthermore, the distance over which the velocity increases considerably moves to lower distances due to the decreasing magnitude of the r_H parameter.

Zoom In Zoom Out Reset image size

Schwenn (1990) analyzed in situ velocity measurements from the Helios mission spacecraft during a period when the Helios-A and Helios-B spacecraft were in radial alignment, showing that the fastest solar wind showed little increase in solar wind velocity between 0.3 and 1 au. This work was later revisited by Maksimovic et al. (2020) who applied a linear acceleration to the solar wind velocity between these heights. The slowest solar wind accelerated at a rate of 90 km s⁻¹ au⁻¹. When a similar approach of a linear solar wind acceleration between 0.3 and 1 au is adopted for the model seen in Figure 3, the slowest solar wind accelerates at a rate of ∼75 km s⁻¹ au⁻¹ over this heliocentric distance range. This provides an improvement on the original HUXt inbuilt acceleration model, which for the same initial velocity at 0.3 au, gave an acceleration of ∼14 km s⁻¹ au⁻¹: a significant underestimation compared to both the findings by Maksimovic et al. (2020) and our PSP-based acceleration model.

As the acceleration of the solar wind is now dependent on the solar wind conditions at the inner boundary, we must now optimize the inner boundary conditions to reflect this. We revisit the tomography-derived inner boundary condition initially proposed by Bunting & Morgan (2022) and further developed by Bunting & Morgan (2023). In the latter study, an empirical conversion function was derived to convert tomographic densities into solar wind velocities at 8 R_⊙. However, because of the changes in solar wind acceleration, this inner boundary condition must be re-optimized. As in Bunting & Morgan (2023), we use a quasi-exhaustive fitting method with the aim of optimizing the V_max, V_min, and α parameters in the density–velocity conversion seen in Equation (1). In this method, a statistical approach is used to fit the distribution of density-derived velocities to a reduced in situ velocity distribution.

The in situ measurements of velocities are backmapped from 1 au to 8 R_⊙; this is done by taking into account the acceleration. We use the updated acceleration model presented in Section 3.3 to find the specific acceleration parameter and, thus, percentage velocity increase between 8 R_⊙ and 1 au. We then reduce the in situ measured velocities by this percentage, giving a distribution of expected velocities at 8 R_⊙ which is compared to the density-derived velocity distribution. An example of an initial in situ velocity distribution can be seen in the left panel of Figure 4, with the backmapped distribution of the in situ distribution shown in the right panel.

Zoom In Zoom Out Reset image size

The fitting method adopts a quasi-exhaustive approach in that, from an initial range of 160–240 km s⁻¹ for V_max, 80–140 km s⁻¹ for V_min, and 1–12 for α, optimal values are found. These initial ranges are determined based on the initial acceleration parameters of the model seen in Table 1. Up to this point, the solar wind velocity at 8 R_⊙ has been constrained between 100 and 180 km s⁻¹ due to rough extrapolation from the PSP data shown in Figure 1. However, in the context of the fitting method, it is important to allow the fitting algorithm to generate model parameters that give the best possible fit without being strongly constrained on the initial extrapolated values. Therefore, the initial parameter ranges to be tested are chosen to be intentionally wide.

Each parameter range is split into 10 separate bins. We then exhaustively generate a density-derived velocity distribution with every possible combination of the parameter bins. Each modeled distribution is then tested and scored against the backmapped in situ data via the mean absolute error (MAE). The parameter bin which gives the best agreement (lowest MAE) between modeled and backmapped distribution is then chosen and the range is updated to the two adjacent bins for each parameter. This process is repeated until the α bin range becomes less than 0.01 and the V_max and V_min parameter ranges become less than 1 km s⁻¹, or until the process has been repeated 11 times (Bunting & Morgan 2023).

The quasi-exhaustive fitting method presented in Section 3.4 is applied via a sliding window approach, with a "window" of 1 yr's worth of tomography data generating velocities compared with one year's worth of in situ data and a "sliding" in three month increments. Once completed, the average of the fitted values between 2018 October and 2021 September are recorded, and are 5.16 ± 1.09, 183.04 ± 6.13 km s⁻¹, and 107.80 ± 5.45 km s⁻¹ for α, V_max, and V_min respectively. The errors are one standard deviation of the values over time.

This result of the fitting process are utilized to generate inner boundary conditions for single CRs. A slice of tomography maps of ±7° from the solar equator is taken. The tomographic densities are then converted into velocities via Equation (1), with the fitted values of α, V_max, and V_min stated above. This generates an inner boundary condition that drives a 3D HUXt model run. 3D is chosen instead of the faster 2D model run in order to accommodate the location of the PSP spacecraft. The PSP orbit varies in latitude by over 3° in a single CR. The efficiency of HUXt allows the solar wind velocity to be extracted at any given time at any location with a computational time of only a few tenths of a second. Therefore, to extract and accurately model the PSP spacecraft orbital path, we extract the solar wind velocity at the location of the PSP spacecraft in hourly time steps. The results are presented in Section 4.2. For model predictions at the location of Earth, a 2D HUXt run is used. The comparison of modeled velocity and in situ OMNI data is presented in Section 4.3.

The time series comparison between PSP observational data (orange) and modeled velocity (blue) at the location of PSP is shown in Figure 5. There is a strong correlation between the observed and modeled velocities at PSP for all time periods. The model accurately predicts prolonged periods of relatively slow solar wind (2019 March–April) as well as periods of more variable solar wind conditions (2020 December–2021 February). The statistical analysis are promising, with an MAE of 64.3 km s⁻¹ and root mean square error (RMSE) of 87.8 km s⁻¹. This provides an improvement on the density–velocity-derived inner boundary condition optimized for the inbuilt acceleration of HUXt (Bunting & Morgan 2022), with the point-to-point metrics of these data being 64.8 km s⁻¹ for MAE and an RMSE of 92.7 km s⁻¹ over the same time period. These suggest that the tomography–HUXt model, with the updated acceleration model, provides a reasonable estimate of the solar wind conditions at heliocentric distances less than 1 au. Further statistics can be found in Table 2. The small-scale variations in the PSP data closely align with the general trend predicted by the model, with the general magnitude of the fast solar wind streams (e.g., 2018 November) showing a good agreement. Also, modeled periods of uncharacteristic slower solar wind (2021 August) demonstrate a good estimation of the in situ data. Due to the data gaps, it is difficult to suggest that the model predicts the solar wind conditions well for all times, but certainly something can be said for the strong agreement during periods where PSP data are recorded.

Zoom In Zoom Out Reset image size

Figure 6 shows a direct comparison between the model velocities at the location of PSP with heliocentric distance at the same times as shown in Figure 1 for models optimized using the inbuilt acceleration model as presented in Bunting & Morgan (2023) in (b) and optimized on the updated acceleration model shown in (c).

Zoom In Zoom Out Reset image size

The difference between the two acceleration models are demonstrated clearly in Figure 6. The adapted acceleration model histogram (Figure 6(c)) shows a strong noticeable increase in velocity with distance when compared to the inbuilt histogram (Figure 6(b)). A comparison between the two histograms shows that the velocity distributions at greater heliocentric distances agree. This is expected as both models are optimized on gaining velocity distributions at 1 au. However, the differences become more obvious at smaller distances. The inbuilt acceleration model gives a minimum velocity at the inner boundary of around 250 km s⁻¹ while the adapted acceleration model gives velocities that rarely exceed 220 km s⁻¹.

When the adapted acceleration model histogram (Figure 6(c)) is compared to the in situ data shown in Figure 6(a), the modeled data show a stronger concentration of slower solar wind velocities, specifically in the bins between 110 and 120 R_⊙. However, there are areas where the higher solar wind values seem to better match the range of velocities seen in Figure 6(a)), specifically ∼150 R_⊙ and greater. This is again due to the density–velocity conversion model being optimized for measurements at 1 au. Overall, the adapted acceleration model gives a distribution that gives far better agreement with PSP measurements in comparison to the inbuilt acceleration model.

The prime focus of space weather forecasting is accurately predicting the solar wind conditions at Earth. In this section, the model predictions at Earth are tested against in situ data provided by OMNI and processed as stated in Section 3.4. The densities are extracted at the latitude of Earth only. The tomographic densities are then converted into velocities at 8 R_⊙ and used as an inner boundary condition for HUXt. Figure 7 shows the model predictions (blue) compared to the in situ solar wind velocities (black) at Earth.

Zoom In Zoom Out Reset image size

Figure 7 shows that the model makes a good estimate of the ambient solar wind conditions at Earth. There are periods of excellent agreement both in terms of time of arrival and magnitude of the faster solar wind streams, such as the majority of 2021, and 2020 August–December. Statistically, the model provides good predictions of the solar wind conditions at 1 au with an MAE of 52.8 km s⁻¹ and RMSE of 72.4 km s⁻¹. Further statistics are provided in Table 2. There are periods where a fast solar wind stream is observed, but not predicted, by the model such as July 2020. Short periods of disagreement could be due to CMEs in the measurements that are not included in the model. There are occasions which the model predicts a fast solar wind stream that is not observed by OMNI. The main source of error is defects in the tomography density maps, including latitudinal uncertainty. As seen in Figure 2, there are steep gradients in density within the solar corona at the boundary between the equatorial streamers and coronal holes. A small latitudinal or longitudinal inaccuracy in the tomography model, or small deflections of streamers with distance, will cause a large inaccuracy in the boundary velocity, and a subsequent large error in the model solar wind velocity.

Given the modeled velocity profiles, and the tomographical densities at 8 R_⊙, we use mass flux conservation to propagate densities to 1 au by

$$\begin{eqnarray}&&{\rho }_{1}=\displaystyle \frac{{\rho }_{0}{V}_{0}{r}_{0}^{2}}{{V}_{1}{r}_{1}^{2}},\end{eqnarray} \tag{ 3 }$$

where ρ and V are density and velocity at heliocentric distance r. The 0 subscript indicates the inner boundary and the 1 subscript indicates 1 au.

An example of the evolution of the solar wind velocity and proton density is shown in Figure 8. This is shown for two cases showing an example of typical slow solar wind (blue) V₀ = 118 km s⁻¹, and typical fast solar wind (red) V₀ = 180 km s⁻¹. As well as the clear anti-correlation between density and velocity, Figure 8 also shows that, in the case of the fast solar wind, the density is close to a factor of 10 smaller than that of the slow solar wind. This agrees with other models and observations (e.g., Habbal et al. 1997; Schwenn 2006; Allen et al. 2020; Bunting & Morgan 2022).

Zoom In Zoom Out Reset image size

As this approach is applied to the modeled values at the inner boundary and 1 au, a back-mapping approach is used in order to directly relate the same parcel of plasma at the two distances of interest. This back-mapping approach entails ejecting and tracking "test particles" within the HUXt model at the inner boundary in 5° increments. The time taken for the test particles to reach 1 au is recorded and the differences in time of arrival, along with knowledge of the solar rotation rate set by HUXt, allow a longitudinal shift to be calculated and build a relationship between the parcels of plasma at 1 au and the inner boundary.

Once the back-mapping approach is undertaken, a post-processing step is applied to ensure there is no crossing of solar wind streams. This allows the electron density at 1 au to be derived from the tomographic densities and derived velocities via Equation (3). We convert the electron density to proton density by assuming a 90% composition of ionized hydrogen and 10% ionized helium is for a quasi-neutral solar wind. Therefore, the relative abundance between the two species is 9 H⁺ : 11 e⁻. This is used for the final comparison between modeled and observed proton densities which can be seen in Figure 9 with accompanying point-to-point metrics seen in Table 3.

Zoom In Zoom Out Reset image size

There is a strong correlation between the general trends of the in situ and modeled proton densities presented in Figure 9. Statistically both compare with an MAE of 2.9 cm⁻³ and a Pearson Correlation Coefficient of 0.38. These results also give confidence that the tomographic densities are accurate and that our velocity model is a good representation of the solar wind. However, the sharp peaks in proton density observed by OMNI are not all present within the modeled densities. One reason for this is that HUXt provides a relatively smoothed velocity output (see Figure 7), due to the resolution requirement of the model input. The inner boundary consists of a coronal resolution of ∼28. Thus, HUXt is restricted to resolving solar wind structures at Earth with an approximate 8 hr resolution. The reduced physical approach could also introduce discrepancies in the output here when accurately modeling the compression and rarefraction present in a stream interaction region.

Profiles of velocity and density allow for an estimate of proton temperatures. For this, the force equation that describes a steady-state flow of the solar wind through infinitesimal radial tubes is used (Munro & Jackson 1977):

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}+\rho V\displaystyle \frac{{dV}}{{dr}}+\rho \displaystyle \frac{{{GM}}_{\odot }}{{r}^{2}}-\rho D=0,\end{eqnarray} \tag{ 4 }$$

where $\tfrac{{dP}}{{dr}}$ is the pressure gradient, G is the gravitational constant, M_⊙ is solar mass, and ρ D represents the outward force per unit volume on the plasma possibly due to wave–particle interactions. In the following work, we only consider gas pressure and gravitational forces, thus D = 0. Substituting the ideal gas law into Equation (4) gives

$$\begin{eqnarray}&&{T}_{p}\displaystyle \frac{d\rho }{{dr}}+\rho \displaystyle \frac{{{dT}}_{p}}{{dr}}=\displaystyle \frac{-\mu {m}_{p}\rho }{{k}_{{\rm{B}}}}\left(V\displaystyle \frac{{dV}}{{dr}}+\displaystyle \frac{{{GM}}_{\odot }}{{r}^{2}}\right),\end{eqnarray} \tag{ 5 }$$

where μ is the mean atomic weight (∼0.62), m_p is the proton mass, k_B is Boltzmann's constant, and T_p is proton temperature.

We now consider two typical cases of the solar wind: fast and slow. The relationship of velocity and density of these cases can be seen in Figure 8. The fast solar wind case shows an initial velocity (V₀) of 180 km s⁻¹, which accelerates to ∼620 km s⁻¹ at 1 au. In the case of the slow solar wind, the velocity at the inner boundary is 115 km s⁻¹ which accelerates to ∼300 km s⁻¹ at 1 au.

For each of these two cases, the above parameters of ρ and V are used to derive a fitting parameter (rhs of Equation (5)). Using an iterative process, a function of T_p with distance is found that gives the best agreement of the rhs of Equation (5) to the fitting parameter. With the aim of generating a more stable fit, the log of the negative of the rhs of Equation (5) is fitted. The fitted parameters of each case can be seen in the dotted red (slow) and blue (fast) plots in the left panel of Figure 10.

Zoom In Zoom Out Reset image size

A scale-height function of T_p gives a stable agreement to the fitted parameter, given by

$$\begin{eqnarray}&&{T}_{p}={T}_{0}\exp \left(-\displaystyle \frac{h}{{H}_{0}}\right),\end{eqnarray} \tag{ 6 }$$

where T₀ and H₀ are free parameters and h = r − r₀.

In the case of the fast (slow) solar wind, T₀ was found to be 2.5 ×10⁶ K (8.2 ×10⁵ K) and H₀ had a magnitude of 79.0 R_⊙ (104.7 R_⊙). Proton temperature (T_p ) is shown as a function of distance in the right panel of Figure 10.

These results show that fast solar wind has a greater temperature, and aligns more closely with in situ measurements, than in previous studies (e.g., Feldman et al. 1974a, 1974b). Measurements of solar wind temperature have been performed by analyzing ionized helium temperatures from the Helios spacecraft at a range of heliocentric distances between 0.3 and 1 au. Marsch et al. (1982) found that the radial temperature of a parcel of faster solar wind (∼650 km s⁻¹ at 1 au) had an initial temperature of 1.3 MK at ∼65 R_⊙, decreasing to 0.9 MK at 1 au. Figure 10 finds a temperature of ∼1.5 MK at ∼65 R_⊙, decreasing to 0.2 MK at 1 au. Likewise for a parcel of slow solar wind (solar wind speed of ∼350 km s⁻¹ at 1 au), Marsch et al. observed a temperature of 0.35 MK at 65 R_⊙, cooling to 0.15 MK at 1 au. Figure 10 shows a temperature of 0.5 MK at 65 R_⊙, cooling to 0.15 MK. This comparison shows that the model is reasonable for the slow solar wind, but the fast solar wind temperature is greatly underestimated at greater heliocentric distances compared to the work of Marsch et al. (1982). There are two major reasons for this disagreement: (1) the work presented above neglects the wave pressure term which could be significant in heating the faster solar wind (van Ballegooijen & Asgari-Targhi 2016); (2) protons are expected to have a strong temperature anisotropy whereas the above work assumes an isotropic solar wind (Huang et al. 2020).

Elliott et al. (2012) used a statistical approach to derive a relationship between proton temperature and solar wind velocity within the ambient solar wind observed by OMNI. A quadratic relationship gave a good representation between proton temperature and velocity at 1 au. The results suggest a parcel of solar wind traveling at 620 km s⁻¹, such as the example given in Figures 9 and 10, yields a temperature of ∼1 ×10⁵ K, agreeing with our temperature model. Likewise for the slow solar wind case, Elliott et al. (2012) suggested a velocity of 300 km s⁻¹ at 1 au corresponds to a proton temperature of ∼2 × 10⁴ K. This value is significantly lower than our estimation (9 × 10⁴ K).