Measurement of the open magnetic flux in the inner heliosphere down to 0.13AU

The Parker Solar Probe (PSP; Fox et al., 2016) was launched on August 2018 into a highly elliptical heliocentric orbit with an inclination of $\sim 4^{o}$ to the solar equatorial plane. With a sequence of Venus gravity assists, PSP’s perihelion distance decreases over the course of the mission. As of November 2020, PSP has completed six orbits with its first three perihelia at 35.7$R_{\odot}$ (0.166 AU), and its fourth and fifth both at 27.8$R_{\odot}$ (0.129 AU). The most recent perihelion was even closer at 20.4 $R_{\odot}$ (0.09 AU). Thus PSP provides measurements of the heliospheric magnetic field more than twice as close to the Sun as the previous record holder, Helios, at 65$R_{\odot}$ (0.3AU). In addition to this unprecedented radial evolution, PSP also samples longitudinal structure in a unique way: as the PSP approaches perihelion, its orbital velocity increases faster than the co-rotation velocity of the Sun and, as a result, it crosses a threshold, where it first co-rotates and then super-rotates with respect to the solar photosphere. The upshot of this is that PSP samples longitudinal variation very slowly, with two intervals of co-rotation (inbound and outbound) where it measures the same solar wind for an extended period. The downside to this is that each PSP perihelion probes only a small range of solar longitude, although data from its cruise phase at larger heliocentric distances provides measurements all around the Sun.

In this work, we utilise DC magnetic field data from the FIELDS instrument (Bale et al., 2016) ¹¹1Data publicly available at https://fields.ssl.berkeley.edu/data/ from orbits 1-5 of PSP, spanning from August 2018 to July 2020 (data from orbit 6 will soon be available and allow for a further extension of the present study). We utilise the full time-series of the one minute (1 min) averaged data product (B_RTN_1min), in addition to producing histograms of the four-samples-per-cycle ($\sim$4.6 Hz) data product (B_RTN_4_per_cyc) over hour-long and day-long timescales. These data utilise the inertial Radial-Tangential-Normal (RTN) coordinate system, with the R component indicating the radial direction at the spacecraft position, the T component defined as the cross product of the radial direction with the solar rotation axis, and the N component completing the orthogonal triad. In terms of heliographic coordinates, T points along a line of constant solar latitude ($\phi$) and N points along lines of constant solar longitude ($\theta$). While $B_{R}$ is the main vector component of interest in this paper, as discussed at length in Sections 2.3-2.6, knowledge of the full 3D vector is key to understanding the evolution of individual components. Here, we also introduce the angular quantity $\alpha,$ which we refer to as the ‘clock angle’, which is the angle the HMF vector makes with the radial direction when projected into the R-T plane.

The data from PSP and the near-1AU spacecraft considered in this study are summarised in Figure 1. Here, each panel is a 2D histogram of data value on the y-axis and time on the x-axis. The time axis is binned into one day intervals and each column of histogram data shows the distribution of the full resolution data set for each of these intervals. A solid black line threads the modal value of the data and indicates which side of the heliospheric current sheet (HCS) PSP was located on each day. Dashed magenta lines indicate the times of perihelia and allow reference to the other panels. The colour scale of the histograms (light to dark) shows the density of measurement (low to high). The top panel shows PSP’s measurements of the radial component ($B_{R}$) of the Heliospheric Magnetic Field (HMF). A dotted envelope showing $B_{R}=\pm 3nT(1AU/R^{2})$ bounds the data and reveals that, as expected, the predominant variation in this component is proportional to 1/$R^{2}$. It also highlights the lower perihelia distances of orbits 4 and 5, where the envelope and data range gets larger than during orbits 1 through 3. We can also clearly see that no matter which side of the HCS the PSP is located on any given day, there is a non-negligible population of field measurements of the opposite polarity. Particularly at the closest approach, this corresponds to the prominent switchbacks, which were a key early discovery of the PSP (e.g. Bale et al., 2019; Kasper et al., 2019; Dudok de Wit et al., 2020; Horbury et al., 2020). The remaining rows show the quantity of principal interest in this paper, $B_{R}R^{2}$, displayed in the same histogram format over the investigated time interval as measured by PSP/FIELDS, the IMPACT/MAG (Luhmann et al., 2004) instrument on STEREO AHEAD (Kaiser, 2004) and the MFI investigation (Lepping et al., 1995) on the Wind (Harten & Clark, 1995) spacecraft. In these latter three panels, a solid red line shows the variation of each spacecraft’s heliographic latitude with time. We note PSP’s minimum latitude and minimum heliocentric distance are closely related.

Finally, to demonstrate the relation between the full cadence and one minute (1 min) average data products used in this work, 1D histograms of both data products are shown overlaid in a series of inset panels on the right. For 1AU, the distributions are virtually identical. For PSP, the one minute average distribution is slightly distorted compared to the four-samples-per-cycle data. In particular, the positive tail of values is foreshortened by the 1 min averaging data but the peak values are unchanged. Thus, we look to use the more pristine faster cadence data product where possible, however, we argue the distributions are similar, especially with regard to the data set of negative polarity, such that using the 1 min average data where necessary will not present an issue. For computational tractability, the 1 min averaged data will be used when considering the whole set of five orbits, while the higher cadence data will be used when considering one-hour or one-day intervals of data.

From the top panel, the data coverage is immediately apparent, demonstrating that there is no PSP/FIELDS data outside of the ‘encounter’ phases of encounters one, two, and three, while for orbits 3 to 5, there is continuous data coverage. At all three spacecraft, there is a striking correlation between latitude variation and dominant measured polarity: When the spacecraft are at their minimum (maximum) latitudes, they are more likely to be southwards (northwards) of the HCS and, therefore, more often measuring negative (positive) polarities for $B_{R}$. For STEREO A and Wind, whose orbits have a constant angular velocity with respect to the solar co-rotating frame, the variation from two timescales is immediately apparent: slow annual latitude variation and regular HCS crossings due to the $\sim 27$ day Carrington rotation period. For PSP, at its closest approach, it can be seen that a single dominant polarity is measured for extended periods of time while PSP is close to co-rotation with the Sun, whereas for the aphelion periods between the latter orbits, the familiar Carrington rotation pattern reappears.

A dotted line at $\pm$ 3nT AU² in the latter three panels shows the same envelope as in the top panel with the radial variation scaled out. Zooming in this way, we see the 1AU spacecraft regularly measuring values of $B_{R}R^{2}$ exceeding this value, while PSP’s measurements are generally below this value, especially during the perihelia. While it is clear $B_{R}R^{2}$ is approximately of constant magnitude and similar in value at all three spacecraft, there is significant scatter and temporal variation. The remainder of this paper is concerned with analysing these distributions of $B_{R}R^{2}$ and the full HMF vector to justify a statistical ‘background’ value and obtain a best estimate of $\Phi_{H}$. We assess the extent to which it is conserved throughout the inner heliosphere and, finally, how it compares to estimates from PFSS models.

We first consider the radial variation of $B_{R}$ directly to confirm the $1/R^{2}$ transformation is appropriate for the full data set. The measurements as a function of radius are plotted in Figure 2 on a log-log scale. Faint blue dots show the spread of the 1 min averages of $B_{R}$ from PSP’s first five orbits (the lower cadence ‘raw’ data), while orange dots show the most probable value of $B_{R}$ for each hour, demonstrating the most probable values of the raw data cluster on a 1/R² trend line. The black markers indicate the most probable values for each day and the different shaped markers represent the different orbits and are shown to include a data product which is sparse enough for the reader to see the difference between the different PSP encounters. Solid red and blue lines show $\pm 2.5nT(1AU/R)^{2}$ trend lines. While these trend lines are not fitted and serve here only as a visual aid, as we show in Section 3, $\Phi_{H}=B_{R}R^{2}=2.5\text{nT AU}^{2}$ actually turns out to be the best estimate for this data set. It is worth noting that while the 1AU data appears to be significantly more spread than the PSP data, this is mostly a distortion due to the combined effects of the higher data density resulting from the limited radial variation of these spacecraft, the symmetric log scale, and the population of vectors where the full field vector fluctuates into the radial direction, such that $B_{R}=|B|$ instantaneously; at 1AU, $|B|$ is significantly larger than the background $B_{R}$ value.

This effect aside, we observe: (1) the radial component of the field indeed varies as 1/$R^{2}$ for all radii probed by PSP to date. (2) the scatter around the trend line is proportionally smallest at the closest in heliocentric distances. It is also systematically skewed towards zero but this is not trivial to observe on the log scale. (3) Orbits four and five show a significant departure from the trend line towards weaker field strengths, particularly on the positive field branch around 0.2AU. This is likely due to PSP running very close to a very flat heliospheric current sheet during these latter orbits (Chen et al., submitted to this issue), and thus sampling more weakly magnetised streamer belt plasma. We conclude using (1), that the quantity $B_{R}R^{2}$ is a useful quantity which is, at least approximately, conserved throughout PSP’s mission so far. There is no suggestion of a break down of the Parker expectation or the Ulysses result of latitudinally isotropised radial magnetic field and suggests the dipole-dominated coronal field has fully relaxed into this isotropic state, which is well within PSP’s closest heliocentric distance of $28R_{\odot}$.

With the investigation of the quantity $B_{R}R^{2}$ remaining well-motivated, we turn our attention to measuring it directly. A key assumption made in displaying the data above, as we have done so far, is to show the most probable value of the field. That is, we take all the data from some time interval or other binning procedure, produce a distribution of that binned data, and assume the peak of that distribution (the bin with highest number of counts) represents the value we are trying to measure.

In Figure 3, we demonstrate how the field is distributed in the R-T plane at different heliocentric distances. We very coarsely bin all the data shown in Figure 1 into four categories: PSP data for R $<$ 0.3 AU (specifically, PSP data interior to the closest perihelia of Helios), PSP data for 0.3AU $<$ R $<$ 0.6AU, 0.6 $<$ R $<$ 0.95 AU, and finally, all the STEREO A and Wind data (for which R ¿ 0.95 AU). The data shown here reflect the 1 min averages due to the aforementioned computational tractability issues when working with the whole data set. The top panel for each radial bin shows the 1D distribution of $B_{R}R^{2}$, the bottom panel shows 2D distributions of $B_{T}R^{2}$,$B_{R}R^{2}$. We note that we have made a replacement here, namely, $B_{R}=|B|\cos\alpha$ and $B_{T}=|B|\sin\alpha$. Here, $|B|=\sqrt{B_{R}^{2}+B_{T}^{2}+B_{N}^{2}}$ is the field magnitude and $\alpha=\arctan 2(B_{T},B_{R})$ is the field vector angle of rotation in the R-T plane relative to the radial direction, and will be referred to as the ‘clock angle’. This transformation is discussed in Appendix A. It is done to avoid the effect of projecting fluctuations in the normal direction onto the R-T plane, which can lead to an underestimation in field magnitude. In particular, here we are assuming the normal, tangential, and compressive fluctuations are uncorrelated, and so making this correction to the R-T components does not affect the $|B|$ and $\alpha$ distributions.

We also note that in Figure 3, the two peaks values are asymmetric. This is simply a sampling effect based on the spacecraft orbit (especially the heliographic latitude). For example, PSP’s orbit is tilted with respect to the solar equatorial plane such that at closest approach, it is also approximately at its minimum (and most negative) latitude. For this reason, PSP for R ¡ 0.3 AU is primarily southwards of the HCS and, therefore, it samples a negative-polarity magnetic field. We can confirm this via an inspection of the second panel of Figure 1 in which the dips where PSP goes into negative latitudes correspond to protracted measurements of negative polarity.

Figure 3 shows that the 1D distributions of $B_{R}R^{2}$ are the projection of a 2D bi-modal distribution which evolves with radius. This 2D distribution is aligned with a mean Parker spiral direction which becomes closer to radial with decreasing heliocentric distance. The distributions comprise of a sunward (upper left quadrant) and anti-sunward (lower right quadrant) population, and exhibiting some spread in field magnitude ($|B|$) and angle ($\alpha$) about this mean state. Empirically, these spreads decrease with heliocentric distance and this trend will be quantified in Section 2.5. Further from the Sun, an overlap of the two populations sum together to distort the 1D distribution of $B_{R}R^{2}$ by producing a large population of values with $B_{R}R^{2}$ close to zero. This effect is much weaker close to the Sun, with the two peaks in the 1D distribution much sharper and better isolated.

In this work, we use distributions of raw (meaning minimally pre-averaged) $B_{R}R^{2}$ and $B_{T}R^{2}$ measurements, such as those illustrated in Figure 3, to produce estimates of the heliospheric flux. We consider three particular methods: the mean, the mode, and the Parker spiral method. The mean and mode are performed directly on the 1D $B_{R}R^{2}$ distributions. The mode simply identifies the bin with the most counts for either $B_{R}>0$ or $B_{R}<0$. For the mean, we take the bi-modal 1D distributions and bifurcate it into two truncated distributions, $B_{R}>0$ and $B_{R}<0$, and the mean of each is taken. At certain points in this work, we consider one-hour or one-day intervals of data. For these cases, the distributions are usually single-peaked (only one magnetic sector explored) and, therefore, the procedure we follow is to identify the single most probable value for that hour or day and then we use the sign of this value to decide which side of the distribution to truncate. The Parker spiral method utilises the full 2D distribution, which we introduce and describe below.

Considering the 1D distributions we see at all radii, the shape of each half of the distribution is asymmetric (or ‘skewed’), meaning that the production of a best estimate of the quantity is not obvious given the mean and mode are non-trivially related (and the dependence of their relationship on radius). The mean in particular is sensitive to the ‘overlap’ of the two sectors. In the simplest approach (e.g. Owens et al., 2008; Linker et al., 2017; Wallace et al., 2019), the modulus of the distribution is taken before computing the mean, which is similar to the truncation procedure we described above. In this method, the larger the ‘overlap’ population relative to the peak, the higher the mean. Given that this overlap population grows with heliocentric distance (see Figure 3 above and also Erdős & Balogh (2012)), this may cause the apparent heliospheric flux to grow with radius when this method is used (Smith, 2011).

In a more careful treatment of the bi-modal population, we would need to subtract the population of the opposing sector, for example, by fitting a curve to the distribution, or to make another approximation, such as bisecting the data in the 2D distribution with an approximate Parker spiral, and to assume that no data fluctuates past 90 degrees from the mean Parker spiral direction (Erdős & Balogh, 2014). An even more sophisticated approach would be to use the electron heat flux to delineate these two populations such that intervals where the electron heat flux is parallel (anti-parallel) to the magnetic field correspond to anti-sunward (sunward) sectors (see Section 3.2 of Macneil et al. (2020)). As mentioned at several points in this work, this application of the electron strahl remains an important avenue for exploring additions to the analysis we describe here. We omit the usage of electron strahl measurements due to the more intermittent temporal coverage of solar wind plasma measurements by the PSP Solar Wind Electrons Alpha and Protons (SWEAP) instrument (Kasper et al., 2016) as compared to the full FIELDS data set analysed here.

The mode is more robust in that it is less affected by the two sectors overlapping and has the advantage that it can easily be defined for both sunward and anti-sunward sectors. It is still possible the overlapping populations could change the peak value somewhat. The mode also requires a large sample size in order to be well-defined and can have a large error when computed from a distribution with a flattened peak. In all these regards, we see the mode is better defined for the distribution of $B_{R}R^{2}$ closer to the Sun, where the peaks are narrower and the impact of overlap from the positive sector is very weak.

In light of these subtleties, Erdős & Balogh (2012) proposed utilising the apparent symmetry in the 2D distributions directly by projecting the data along the nominal Parker spiral direction, computing the mean of data along this direction and projecting this into the radial direction to obtain a best estimate of $\Phi_{H}$. As shown by the shape of the 2D distributions in Fig 3., this approach is supported by PSP data and even more so close to the sun where the peak in the 2D distributions become very sharply defined (even though, for this data set, they are predominantly in the sunward sector). Moving forward, we refer to this technique as the ‘Parker spiral method’, or ‘PSM’. In this work, we use an empirically measured Parker spiral angle rather than the ideal angle generated with a solar wind velocity as done by Erdős & Balogh (2012). This choice is made because, as mentioned above, the solar wind velocity data set is a subset of the magnetic field data set from PSP that is due to differing instrument operation schedules and constraints. Thus, it is desirable to use a method which can be performed purely with magnetic field observations. Furthermore, the Parker spiral is an idealised model and, in fact, driving a model with a varying solar wind velocity leads to some unphysical inferences, such as slower wind streams overlapping with faster ones (which form stream interaction regions in more realistic simulations). Additionally, recent observations by PSP (Kasper et al., 2019) call into question the Parker spiral assumption that the solar wind velocity flow is purely radial at distances probed by PSP so far. By empirically estimating the background vector field, we avoid the need to assume an exact Parker model or to use more sophisticated models. Finally, an empirical Parker spiral angle allows us to work directly with the vector in a spherical representation when implementing the PSM, which is preferable to the Cartesian representation, as described below.

To implement the PSM, we form 1D distributions of —B— and the clock angle ($\alpha$) which underlie the 2D distributions in the R-T plane shown in Figure 3 and measure the most probable value of these quantities, $B_{0}$, $\alpha_{0}$. When the distributions are bi-modal, a most probable sunward and anti-sunward value of $\alpha_{0}$ is found. The most probable values are computed by identifying the bins of these 1D distributions of $|B|$ and $\alpha$ with the most counts. Then the estimate of $\Phi_{H}$ for the PSM is derived as $\Phi_{H}=B_{0}\cos\alpha_{0}R^{2}$ where $R$ is the spacecraft heliocentric distance relevant to that distribution. The separation of these 2D distributions into 1D magnitude and angle distributions are illustrated in Figures 12 and 13 in Appendix B.

Is is important to discuss the effect of pre-averaging on our different estimation methods. When the HMF vector field is pre-averaged in Cartesian coordinates, the whole 3D distribution is deformed by more and more aggressive time-averaging and this procedure affects all the methods, including the PSM and the mode. In particular, the effect on bi-modal distributions is to shift the peaks towards zero, causing estimates of the heliospheric flux to decrease, as observed for all methods by Owens et al. (2017). In terms of the 2D distributions depicted in Figure 3, the effect of longer pre-averaging on the raw data would be to shrink the whole annulus of data towards the centre (see Appendix C). In terms of quantifying this, Owens et al. (2017) showed one hour and one day averages produced significantly lower values of $\Phi_{H}$ (a reduction of $\sim 20\%$, although this only applies to 1AU data) and we should note that Erdős & Balogh (2012) used a six-hour average for their PSM implementation. In this work, we attempt to avoid pre-averaging issues primarily by applying our estimation methods to distributions of 1 min averaged or higher cadence data such that our vector distributions are minimally distorted. However, we also note one very important characteristic of our PSM implementation is that we parameterise our vector field in spherical coordinates at very high time resolution. The resulting parameterised distributions are evenly distributed about a central value for most radii. This means that as long as the field stays in one sector (sunward or anti-sunward) for most of each averaging window, these time series of $|B|$ and $\alpha$ should be able to be arbitrarily ”pre-averaged” without affecting the results. This is demonstrated with STEREO A data in Appendix C. As mentioned above, use of the electron strahl would be a robust way to delineate the sectors a given HMF measurement belongs to.

We reiterate here that in certain parts of this work, our flux estimation methods are applied: (1) to data over one-hour intervals; (2) data over one-day intervals; (3) data from the full mission binned by heliographic location. For case (3), there are sufficient data points using 1 min averages to establish robust histograms to study. For cases (1) and (2), we utilise the four-samples-per-cycle PSP/FIELDS data product, providing N = 16500 data points for each hour and over 400,000 for each day, ensuring well-formed histograms for these cases too. We next move on to compare and contrast our different estimation methods introduced above, starting by applying them to an idealised synthetic data set.

To investigate how the underlying 2D distribution of magnetic field data affects the statistics of the 1D distributions of $B_{R}R^{2}$ and how this affects our different ways of measuring $\Phi_{H}$ introduced in the previous section, we begin by producing a synthetic data set. Motivated by Erdős & Balogh (2012), we suggest that the 2D field be modelled in a polar representation as a set of normally distributed fluctuations in vector magnitude ($|B|$) and R-T clock angle ($\alpha$) with standard deviations $\sigma_{|B|}$ and $\sigma_{\alpha},$ respectively, which may be thought of as encoding the relative balance of compressive and rotational (or Alfvénic) fluctuations respectively. It is important to state we are making a strong assumption that fluctuations are evenly distributed about an average value which also corresponds to the mode and the ‘central value’. This assumption allows us to establish a set of quantities ($|B_{0}|$, $\alpha_{0}$, and $B_{R0}=|B_{0}|\cos\alpha_{0}$), which can be defined as a ‘ground truth’, and which the PSM measures directly by our method of constructing this synthetic data.

In Figure 4, we generate ensembles of synthetic measurements by drawing from normal distributions using the random module in the numpy library (Harris et al., 2020) in python over a range of background clock angles and fluctuations, and project them into synthetic 1D $B_{R}R^{2}$ distributions as in Figure 3. Because we generate the sunward and anti-sunward sectors separately, we can keep track of the individual distributions (red and blue curves) and see how they add up to give the overall distribution. It is immediately apparent that the overlap grows with radius simply because of the increase in the Parker spiral angle (c.f. Erdős & Balogh, 2012, 2014).

We note that for a central value of $|B_{0}|$, the dimensionless quantity $\sigma_{|B|}/|B_{0}|$ is a useful figure of merit which may be compared to $\sigma_{\alpha},$ as expressed in radians to describe the balance of fluctuations. In the top row of Figure 4, distributions dominated by rotational fluctuations at various background clock angles ($\alpha_{0}$, as indicated by the white dashed line) are shown, while in the bottom row, compressive and rotational fluctuations are balanced. By comparing Figures 3 and 4, we can see that the ‘boomerang’ shapes of fluctuations dominated by rotations are heuristically more similar to the real PSP data. This will be established quantitatively in the next section.

As expected, the Parker spiral method exactly reproduces the ‘ground truth’ as indicated by the blue and red circles aligned with the spiral direction in the 2D distributions. The solid (red or blue) vertical lines indicate the value of $\Phi_{H}$ that would be measured by the Parker spiral method, and that link the 2D and 1D distributions. Dotted and dashed vertical lines respectively respectively show the measured mode and mean based on the 1D summed distribution (black curve). The mean for both positive and negative sectors is computed by truncating the distributions at $B_{R}R^{2}=0,$ as discussed in Section 2.3. To be clear, these operations are applied directly to the distributions depicted here which represent the ‘pre-averaged’ or base data product. At all points in dealing with real data, the maximally pre-averaged data we use are the 1 min averages.

By comparing the vertical lines in Figure 4, we see that our three methods give slightly different results and the specific differences depend on both the nature of the fluctuations and the mean spiral angle. In particular, when rotational fluctuations dominate, the mean value is an underestimation, especially when the Parker spiral is close to radial (left-most column). At these spiral angles, the mode is also a slight underestimation; however, at a higher Parker spiral angle (right-most column), the mode appears to overestimate the true value. The skewed 1D distribution is being primarily shaped by the geometry of projecting angular fluctuations and so, the distribution peak does not in general correspond to the 2D background value about which the field vector is fluctuating. Thus, the assumption that the mode is ‘more representative’ than the mean is not yet justified.

When the compressive and rotational fluctuations are balanced, all three values are closer together (the mean is a better approximation), especially at a higher Parker spiral angle. In this case, for low Parker spiral angles, the mean is, in fact, higher than the mode.

Overall, in this section, we show that the nature of the 2D fluctuations of the vector and the inclination of the Parker spiral angle can affect our measured value of the heliospheric flux. Therefore, if the vector fluctuations can indeed be modelled as normally distributed about a mean magnitude and angle – especially if rotations dominate over vector magnitude fluctuations and then the Parker spiral method is clearly the most robust and physically motivated method considered here. The fluctuation balance would also determine the appropriateness of each method. Next, we seek to measure the background vector and fluctuation statistics as a function of radius from the real spacecraft data and we also check the extent to which fluctuations are evenly distributed.

In Section 2.4, we found that for idealised synthetic data with symmetrically distributed 2D vector fluctuations, the nature of these fluctuations and the large-scale inclination of the Parker spiral could affect our estimates of the heliospheric flux. In particular, how the fluctuations were partitioned into angular and compressive (—B—) fluctuations had an impact. In this section, we seek to measure the inclination of the Parker spiral and the distribution of 2D fluctuations to assess how these effects manifest in estimates of the heliospheric flux with real data.

In Figure 5, we compute the statistics of the field magnitude (scaled by $R^{2}$) and clock angle and then we plot the results as a function of radius. We do this by binning all the data by radius (large squares) as well as by binning in time at a daily and hourly cadence (faint and fainter circles). The top row shows the most probable value, the second row shows the standard deviation (for the field magnitude, scaled by the mode to make it dimensionless), and the bottom row shows the skew, defined as the (mean - mode)/standard deviation.

In the top row, we also plot the expectation of these quantities for a Parker spiral model as a function of radius. The analytical expressions (Parker, 1958) are :

Here, $v_{SW}$ is the solar wind velocity, $\Omega_{\odot}$ is the solar equatorial rotation rate (14.17 deg/day), and $+(-)$ indicates the anti-sunward (sunward) sectors. The models shown use $B_{0}=2.2nT$ and $v_{SW}=300km/s$. As can be seen, in all binning schemes, the most probable value of the data is in good agreement with the model across different radii, with larger scatter for the shorter binning intervals.

For the standard deviations (middle row, Figure 5), which describe the relative strength of fluctuations, we see both compressive and rotational fluctuations grow in amplitude with radius, but we also see that at all radii, the amplitude of rotational fluctuations is larger than compressive fluctuations. Binned by radii (such that the distribution is formed over many streams), rotational fluctuations are about three times as large as the compressive fluctuations at closest approach and slightly less than two times as large at 1AU.

The pattern in the skew is less clear. At a one-hour or one-day cadence, there is no trend: at all radii, there are as many distributions skewed positive as negative. For vector magnitude, averaged over many streams, there is a hint of negative skew at perihelion and a small positive skew at 1AU. The negative skew at low heliocentric distances arises from the near-HCS weaker field mentioned in Figure 2; many single-day intervals at these radii have very low skew. For the clock angle, the skew at perihelion is negligible, but trends positive at 1AU, which means the fluctuations are biased in favour of rotation towards the radial direction.

The relationship of the different timescales considered here and the underlying data products warrants some discussion in passing. Firstly, we note explicitly that the one-hour and one-day cadence data points are computed from the higher cadence four-samples-per-cycle data product, while the data binned by radius (white squares) are based on the 1 min averaged data. The fact that the white squares lie in the middle of the distribution of one-hour and one-day measurements fuels our confidence that the two data products are largely interchangeable. This also shows that at a given radius, non-consecutive hour or day intervals from different orbits with no temporal correlation fill out a well-defined distribution around the relevant white square. This is useful evidence that the central value is not tied to a specific time interval and it is, therefore, real spatial information about the HMF, meaning it is justified to use different timescales for estimating the heliospheric flux at various points in this paper.

To summarise, in this section, we successfully measured the 2D vector fluctuations and underlying large-scale radial trends in the HMF that, as we showed in Section 2.4, could impact our measurements of the heliospheric flux. In addition, we characterised the skew of the distributions, broadly assessing the requirement for the vector fluctuations to be evenly distributed about a central value to justify the use of the Parker spiral method. Additionally, based on this, we infer that: (1) the Parker spiral is a good model for the underlying large-scale structure; (2) at all radii from 0.1-1AU, the amplitude of rotational fluctuations are stronger than compressive ones ($\sigma_{|B|}/|B|<\sigma_{\alpha})$, but this predominance decreases with heliocentric distance; (3) the Parker spiral method for computing $\Phi_{H}$ (which relies on having an unskewed distribution in magnitude and clock angle) may be less applicable at 1AU, where the distributions averaged over streams appear to have a positively biased skew in both polar components (and the central assumption of evenly distributed fluctuations breaks down). The field magnitude over multiple orbits also showed some signature of skew close to the sun, but this can be attributed to distortion by the weaker field near the HCS during encounters four and five. This is accounted for by taking the most probable value of the magnitude rather than the mean. For reference, the radially binned distributions which give rise to the statistics shown in this section are shown in Appendix B. Based on Figure 4 and the above observation that the rotational fluctuations are dominant, we expect the mean to be a systematic underestimation of the background value of $B_{R}R^{2}$, especially close to the sun.

Having measured the large scale variation of the HMF and the characteristics of the vector fluctuations, we next seek to quantify the expected and measured variation in the statistical methods of estimating $B_{R}R^{2}$ as a function of radius. In Figure 6, we compare measurements of the mean (diamonds), modes (squares), and the Parker spiral method (circles) from purely synthetic distributions (top panel) and directly from measurements (bottom panel).

For the synthetic measurements, we take the Parker spiral model values for the background values of $|B_{0}|(R)$ and $\alpha_{0}(R)$ from Figure 5 (yellow curves, top panels) and for the fluctuation amplitudes, we take the linear fits to the standard deviations from the middle panels of Figure 5. As earlier, even though we compute the two sectors separately for the synthetic data, we use the joint distribution projected into 1D and split at $B_{R}R^{2}=0$ to compute the mean and mode. This method simulates the process of obtaining these measurements from the real data where the components of the two magnetic sectors cannot be distinguished in the 1D distribution. We do actually allow the synthetic distribution in magnitudes to have a non-zero skew further from the sun (see Appendix B since fitting a standard deviation from these skewed distributions would produce an nonphysical Gaussian distribution, where a significant population of magnitude values would go through zero. In these cases, we are modelling our ‘ground truth’ as the peak of the skewed distribution of magnitude.

As expected, for the purely synthetic data, the Parker spiral method produces the same value of magnetic flux for all radii. While this background value remains constant, the estimated mean and modes evolve with radius due to the changing impact of the fluctuations and background Parker spiral angle. For low radii, the mean produces a constant value but significantly underestimates the magnetic flux. Around 1AU, the mean systematically increases relative to the background value, and in this model becomes an overestimation. Closer than about 0.8 AU, the mode is a slightly better estimate of the background value but it fluctuates, at some points, resulting in either an overestimation – or sometimes an underestimation. At 1AU, the mode also systematically increases and overestimates the ‘true’ value.

For the real data, the PSM is performed directly with the most probable values of —B— and $\alpha$ taken from the white squares in Figure 5 to show representative values integrated across multiple orbits. These are derived with the 1 min averaged data products. We note that the yellow and magenta filled data points are from STEREO A and Wind respectively, while the remainder of the data is from PSP. For completeness, we show both the positive and negative sectors measured by PSP, although the trends discussed below are primarily observed in the negative polarity (blue data points). This is due to the limited time intervals where PSP has been above the HCS, particularly at the closest approach. Furthermore, the dominant data contribution above the HCS for radii closer than 0.3AU comes from encounters four and five, which were distorted for extended periods by PSP skirting very close to the HCS (Chen et al., submitted to this issue).

As faint data points in the background, PSM estimates of $B_{R}R^{2}$ for one day time intervals are shown from the four-samples-per-cycle data, demonstrating the variability of the estimate from one day to another and that the values of the estimate integrated over multiple orbits lie in the centre of this distribution. This distribution is more spread than the difference in the central value from one radial bin to the next.

Although the trends with real data are unavoidably noisier, when focusing on the negative sector, the mean can be clearly seen to exhibit the same behaviour as in the synthetic data, producing a systematic underestimation close to the Sun, but growing and exceeding the Parker spiral method at 1AU. In the positive sector, the same behaviour is seen, but less clearly for R ¿ 0.3AU, with a flat mean eventually trending upwards for R ¿ 0.8AU. Importantly, the trends implied by the PSP data converge nicely with the 1AU data. For the negative polarity data trends the mode is also in good agreement with the synthetic observations, providing better agreement with the Parker spiral method close to the Sun, but also systematically increasing near 1AU and at intermediate radii overestimating the PSM. The PSM, while not perfectly flat, is generally quite consistent with the radius, but not dramatically more so than the other methods. For $B_{R}R^{2}<0$ near 1AU, there is a significant drop-off in the estimate of $B_{R}$ using the Parker spiral method. This may be related to the growing skew exhibited in the bottom panel of Figure 5 and suggests the method may be less robust near 1AU. As mentioned previously, our implementation of the PSM relies on the fluctuations (particularly in the angle) being evenly distributed about a mean. The bottom panels depicting the skew in Figure 5 suggests this assumption starts to breaks down close to 1AU.

Overall, we conclude that this model of fluctuations in magnitude and angle is a good approximation for understanding the dependence on heliocentric distance of the different methods of computing $\Phi_{H}$ considered here, at least for R ¡ 0.8 AU. In particular, Figure 6 suggests that fluctuations cause the mean of PSP data close to the Sun to constitute an underestimation – and to be an overestimation beyond 0.8AU. Below 0.8AU, at least for the negative polarity measurements, the PSM and mode estimates are comparably variable with radius but show the ordering expected from our model of fluctuations – with a higher value of mode as compared to PSM and a convergence of the values close to the sun. Comparisons of the synthetic and measured 2D distributions and distributions of $|B|$ and $\alpha$ in 0.1AU bins may be found for reference in Appendix B.

As a result of the findings in this section, we conclude that it is necessary and important to utilise the 2D distribution of the magnetic flux in order to compute the background value of $B_{R}R^{2}$; without this step, the radius-dependent fluctuations can lead to systematic biases in the mean or modal statistics value. Thus, of the methods considered here, the Parker spiral method presents the most robust estimate, although its estimate appears interchangeable with the mode close to the sun and may be less reliable near 1AU. In the next section, we present our results, using the Parker spiral method (and the other statistical methods where justified) to measure $B_{R}R^{2}$ throughout the inner heliosphere and, finally, we compare it to expectations from coronal models.

In Figure 6, near 1AU, at least for the positive polarity data, all three metrics appear to climb to near 3nT AU², possibly indicating a physical enhancement in the heliospheric flux at 1AU. To probe this observation further, as well as to develop a ‘best estimate’ of the heliospheric flux closer to the sun, we sum together all the $B_{R}R^{2}$ distributions over radius and show their bulk statistics. For every hour of data, we compute the mean, mode, and Parker spiral method estimates as described in Section 2.3 and show how these estimates are distributed over the studied time interval. A one-hour timescale is chosen here to produce a sufficient number of individual ‘measurements (on the order of 15000) with each method across the whole data set expected to produce a robust histogram and define a peak value. A longer integration timescale would sharpen these distributions, but limit the statistics. This means our computations of error are likely conservative.

The results are shown in Figure 7. We show the three spacecraft separately – for PSP we include a panel showing only data from close to the sun (R ¡ 0.3 AU) – and overlay the histograms of hourly measurements with the three different methods. The information contained in these histograms is compressed into a most probable value and upper and lower values defining the full width half maximum (FWHM) for both peaks. These data are shown as text in the corresponding colour in the same plot.

We first observe that in all cases, the distributions of measurements using the mode and Parker spiral methods are quite similar to each other, whereas the distribution of mean measurements is strongly distorted: while the truncation procedure (similar to taking the modulus, see Section 2.3) successfully keeps the sunward and anti-sunward peaks distinguishable, the overall distribution is broadened and shifted towards $B_{R}R^{2}=0$, producing a systematic underestimation compared to the other methods and with a larger uncertainty as encoded by the FWHM. We see by restricting the data to PSP measurements close to the Sun (top left panel in Figure 7), the distributions for all methods are the sharpest, although the systematic underestimation between the mean and the other methods persist, as expected from our synthetic data model. The broadening of all distributions with distance from the sun shows the increasingly prominent role of fluctuations and other distorting effects with increasing heliocentric distance.

Distinguishing the mode and PSM is slightly more subtle. With regard to the near-1AU distributions, the Parker spiral method is seen to produce distributions that are more similar to the close approach distributions further from the Sun, doing a better job of removing the population of measurements close to zero and separating the two peaks, as compared to the mode. While this is most clearly seen at 1AU where there is a larger population of $B_{R}R^{2}$ measurements close to zero in the red curve (mode) compared to the black curve (PSM), integrating the PSP data over all radii, (top right panel in Figure 7), the difference is slight but still observable. Once we restrict the data below 0.3 AU (top left panel in Figure 7) both the mode and PSM distributions are close to indistinguishable. This supports our inference, presented in Section 2.6, that the mode is distorted more than the PSM by fluctuations as distance from the sun grows, although for practical purposes, for most of the PSP data, the noise in these distributions means the mode and PSM measurements derived are consistent within the error; the growth of the difference between these methods is most pronounced between the outer ranges of PSP’s orbit and 1AU.

We note for the PSP data, the negative peaks are much better resolved than the positive peak due to the same latitudinal sampling effects and near-HCS intervals which caused the trends to be more apparent with the negative polarity data in Figure 6. Integrated over all radii, the negative polarity PSM peak gives a bulk value of $\Phi_{H}=2.5^{+0.3}_{-1.0}$nT $AU^{2}$ where, as mentioned above, we have computed the error based on the FWHM of this profile. Confined to just R ¡ 0.3 AU, the measurement tightens to $\Phi_{H}=2.5^{+0.3}_{-0.6}$nT AU².

We also see the most probable PSM-derived estimate of $B_{R}R^{2}$ at 1AU is significantly larger ($\sim$2.8 nT $AU^{2}$) than the value derived above from the negative peak of the PSP data, although just within the error. The enhancement is apparent in both Wind and STEREO A data. The trend is also seen in the most probable value of modal values, but not here for the mean. For the mean, this seems to be an issue with defining the mode for a broad peak, looking at the upper FWHM value, this can be seen to clearly increase for all methods in comparing the top two panels to the bottom two. While we expected from the analysis of fluctuations that the mean and mode of $B_{R}R^{2}$ will grow with radius, this cannot explain why the value would be significantly higher than computed with the Parker spiral method. Although we expect the PSM results to be less robust in this vicinity due to apparent skew of the distributions and associated breakdown of the even distribution assumption (Section 2.5), it is not immediately apparent why this should result in a systematic increase of $B_{R}R^{2}$ as opposed to a decrease. Therefore, it appears possible that this flux enhancement could indeed be physical. We investigate a possible explanation for such an enhancement in the next section.

In the heliosphere, numerous physical processes such as waves, turbulence and stream interactions perturb the quiescent picture implied by the Parker model. One particular consequence of such effects is that local field inversions develop wherein the field lines connecting back to the corona warp into an S shape. As argued by Owens et al. (2017), such inversions entail the same field line threading a spherical surface at a fixed radius multiple times (see also Appendix 15 for an illustration) such that a given field line would contribute triple the flux at this radius as compared to its contribution leaving the corona. Therefore, they could contribute to a physical excess flux. The cited authors identify such inversions using the electron heat flux (aka strahl) which follows the HMF topology and carries energy away from that field line’s point of origin in the corona. Therefore, when the strahl is observed to be sunward, it is inferred that the electrons must be travelling along a kinked field line such that they escaped and moved outwards from the sun but have then been guided sunward by a field inversion. This approach was recently applied (Kasper et al., 2019) to probe the topology of the switchbacks observed by PSP. Owens et al. (2017) used this method with ACE data over a solar cycle to derive a correction factor to the heliospheric flux. They found this estimate to be consistent with the values obtained by pre-averaging signed $B_{R}R^{2}$ data over one day intervals before taking the modulus, and slightly weaker than the kinematic correction factor proposed by Lockwood & Owens (2009); Lockwood et al. (2009a, b). Furthermore, Macneil et al. (2020) showed with Helios data that the frequency of such local inversions increased with radius from 0.3 to 1.0 AU, suggesting that this correction factor would grow with radius.

In Figure 8, we investigate this same radial trend with PSP $B_{R}$ data from 0.13 to 0.9 AU. For each day of PSP data, we compute the fraction of the distribution which is of the opposing sign to the most probable value, making the assumption that over each 24-hour period, PSP remains on one side of the HCS and that all zero crossings in this interval are due to fluctuations. This is slightly less robust than the Owens et al. (2017) strahl method which can distinguish HCS crossings unambiguously and account for transient structures such as coronal mass ejections (CMEs) but still carries useful information about the trend. A one-day interval is chosen as the trend is readily apparent with no further processing; computing it with a shorter time interval, such as one hour, is much noisier. In Figure 8, we plot this fraction of inverted flux as a function of time and also show the times of PSP’s perihelia and its radial variation. While the metric is noisy, as might be expected given the above assumptions, there is a clear correlation with radius, in which the fraction of inverted flux grows with distance from the Sun, which is consistent with Macneil et al. (2020). The fraction has a well-defined floor which varies from $<3\%$ at perihelion to as much as 20% at 1AU. It is therefore plausible that if such fluctuations can contribute to excess flux, their impact will grow with radius in the inner heliosphere. However, conversely, we see that at perihelion the contribution is extremely small (¡ 3%). This is interesting given the prominence of the switchbacks observed by PSP. This implies that while these switchbacks are very striking given the large amplitude and sharpness of the rotations (e.g. Horbury et al., 2020), the population of measurements in which the field actually reverses is a very small fraction of the total. They are transient impacts perturbing a quiescent background state.

To numerically estimate the possible impact of the flux, we make a very simple construction, which is as follows: if there are $N$ open field lines emerging from the corona and a fraction $f$ of them are locally inverted at a given radius, then the true open flux per unit area is $d\Phi_{open}/dA=N;$ but the actual measured flux at that radius will be $d\Phi_{H}/dA_{m}=(1-f)N+3fN$. This can be rearranged to relate the true and measured flux : $\Phi_{open}=1/(1+2f)\Phi_{H}$. For $f=3\%$, this simple heuristic implies $<\%5$ open flux estimate error at perihelion and up to a $30\%$ correction at aphelion, which is similar (slightly larger) than the factor derived by Owens et al. (2017), as judged from their Figure 5, and significantly larger than the fraction implied by Figure 7 in the present work, which suggests a $\sim 20\%$ reduction in flux between PSP and 1AU. This metric is likely an overestimation since it assumes the switchback fraction measured near the ecliptic plane, where the solar wind is predominantly slow, can be applied uniformly in latitude but, as shown with Ulysses data (Erdős & Balogh, 2012), fast polar coronal hole wind for most latitudes has much lower fluctuation levels than the slow wind.

Regardless, in terms of contribution to the excess flux, near PSP’s perihelia flux inversions are largely negligible. Combined with the much reduced relative amplitude of fluctuations at closest approach, we expect that the value of heliospheric flux, $\Phi_{H}$, at least at PSP’s closest approach, should robustly correspond to the true open flux escaping the corona, $\Phi_{open}$, with $5\%$ as an upper bound to the possible deviation.

We show here how the Parker spiral method determination of heliospheric flux varies with heliographic position, while examining all three dimensions of radius as well as heliographic latitude and longitude. The results are summarised in Figure 9. In the background, we show 2D histograms of the full data set of 1 min averages (again due to tractability issues of working with the full four-samples-per-cycle data set). In the foreground, we plot the Parker spiral method results computed from each day of data and colour them by the polarity. One-day intervals are used to reduce the scatter and not overcrowd the plot as opposed to one-hour intervals, while spatial binning is avoided so that each data point may simply be assigned a 3D spatial coordinate so that it may easily be plotted against each coordinate (radius, latitude, and longitude). The top row shows results against radius, the middle against latitude, and the bottom row against longitude. Columns from left to right shows PSP, STEREO A and Wind results.

First we see that in all cases, the one-hour estimates trace the regions of highest data density (dark colours in the histograms) well, indicating the close relationship between the Parker spiral method and the mode. Second, plotted against all the coordinates, the PSP data shows the tightest confinement and least scatter due to the much reduced fluctuation levels (see Figure 5, middle row). We confirm the heliospheric flux over these two year intervals is largely independent of radius and longitude, but most clearly shown by the PSP data for $R$ ¡ 0.5AU and for negative latitudes. To emphasise the radial trend, we make linear fits to the data shown (with negative and positive data points fitted separately). These fits are shown as solid black lines in the top left panel of figure 9 and their fitted equations are printed within the same panel. The fit to positive polarity shows a substantial slope due to the positive polarity of the weakly magnetised plasma PSP sampled at it running close to the HCS at closest approach during encounters 4 and 5 (see Figure 2 and Chen et al., submitted to this issue). The negative polarity data is much less perturbed by these effects and benefits from a larger sample size (three whole encounters predominantly below the HCS). For this case, the trend is well fit by a constant value indicating constant heliospheric flux as a function of radius and limiting the variation in the flux with radius to 0.1 nT AU² per AU. In fact, this fit suggests a weak decrease in flux with radius. We note that as suggested by Figure 6, the onset of increasing flux with distance may be non-linear and may inflect close to 1AU.

Over the limited range of latitudes probed by PSP and the 1AU spacecraft, the heliospheric flux is also seen to be constant. It is no coincidence that the latitude and radius panels show very similar trends since these coordinates are strongly correlated due to the tilt of the spacecraft orbits relative to the solar equatorial plane. The high data scatter at high PSP latitudes is simply a result of PSP being at aphelion during the sampling at these high latitudes, where the fluctuation environment is much noisier.

Comparing the longitude plots, we see a similar sector structure pattern measured by all three spacecraft over the full two-year data set, indicating the dominant warps in the HCS were consistent throughout. The predominantly positive sector in the PSP panel (top right) is shifted to a slightly higher longitude than on the 1AU panels, which is consistent with the Parker spiral shift from PSP’s closer heliocentric distance out to 1AU.

The STA and Wind data show a slight asymmetry in field scatter by latitude. At highest latitude, they measure slightly higher positive flux, and at lowest latitude they measure slightly stronger negative flux. The reason is not perfectly clear. From Figure 1, it is clear that the latitude and field values are correlated. Extremes in latitude means the spacecraft is spending less time close to the HCS in general. This might mean faster solar wind with lower fluctuation levels and therefore a measured value which is closer to the ’true’ open flux and would suggest the excess flux is physical.

Lastly, we compare the Parker spiral method measurement results to estimates of the open flux from Potential Field Source Surface models (PFSS Altschuler & Newkirk, 1969; Schatten et al., 1969; Wang & Sheeley, 1992). The specific model implementation is the same as described in Badman et al. (2020) and makes use of pfsspy (Yeates, 2018–; Stansby, 2019–; Stansby et al., 2020), a python implementation of the PFSS model. The inputs of the model are 2D magnetogram maps of $B_{R}$ at the solar photosphere, and some choice of source surface height ($R_{SS}$). The output of the model is the vector magnetic field, $\boldsymbol{B(\boldsymbol{r})}$, for the annular region between the photosphere and the source surface ($1R_{\odot}<R<R_{SS}$). The open flux is estimated from this model by integrating the field vector over the source surface, where all field lines are radial by construction and interpreted to be open to the solar wind. This integral is given by Equation 1, which is repeated here for convenience: :

As a reminder, to make this quantity comparable to the quantity $B_{R}R^{2}$ measured in situ, we normalise it by $4\pi$ (1AU)², such that is is expressed as a field strength in nT at 1AU.

We used the Air Force Data Assimilative Photospheric Flux Transport (ADAPT, Arge et al., 2010) magnetograms ²²2ADAPT maps are accessed from https://gong.nso.edu/adapt/maps/., which assimilate the most recent available photospheric data (from the visible part of the Sun) into a surface flux transport model which forward-models the magnetogram from the previous time-step, taking into account differential rotation and meridional flux transport. In this way, it seeks to model the magnetic field of the entire solar photosphere ‘synchronically’ (meaning all longitudes are modelled at the same time), as opposed to synoptic magnetograms, which simply merge together old and new data. These maps are produced multiple times per day (much more frequently than a Carrington rotation). In this work, we access one magnetogram per day. We utilise ADAPT magnetograms, which assimilate photospheric data from the ground-based Global Oscillation Network Group (GONG; Harvey et al., 1996) and space-based Heliospheric Magnetic Imager (HMI; Scherrer et al., 2012) which are both available at the above URL for the PSP mission duration. The HMI data product is only available for about $\pm 1$ month about each perihelion, while the GONG data product is available for the full duration studied.

In Figure 10, we compare the open flux estimated from the PFSS models to that measured in situ by PSP as a function of time. Results for two different values of source surface heights ($2.0R_{\odot},2.5R_{\odot}$) are shown for ADAPT-GONG and ADAPT-HMI maps. The lower value is chosen based on better agreement with the current sheet crossings measuring during PSP’s first perihelion (Badman et al., 2020; Réville et al., 2020; Szabo et al., 2020), while the higher value is the widely accepted canonical value of $R_{SS}$ (Hoeksema, 1984), and which appears a better fit for PSP’s second orbit (Panasenco et al., 2020).

To compare to the in situ measurements, we here present seven day averages of the hourly Parker spiral method results from Figure 9. This is chosen to avoid overcrowding the plot and emphasise where the data is most self-consistent and why the most probable trend line is located where it is. Error bars show the standard deviation of the hourly measurements within these seven days. In addition, the most probable value and error generated from the bulk distribution of the PSM results (histogram peak and upper and lower FWHM, top right panel in Figure 7) is shown as a solid black horizontal line and light grey region. In addition, the darker grey region shows the smaller error region of the distribution for R ¡ 0.3 AU, which has the same central value (top-left panel, Figure 7). This error region captures all the data points for the first three perihelia well. As noted in Section 3.1, these errors were generated using the distribution of one hour PSM measurements. This short time interval was necessary to produce enough statistics to produce a robust distribution, but as we see in the discussion in Section 2.5, longer time intervals reduce the scatter further and, thus, the error bounds here are likely to be quite conservative. We see that the most reliable measurement of the open flux occurs southwards of the HCS, around perihelion of the first three orbits (November 2018, April 2019 and September 2019). This ‘best estimate’ is shown only for negative polarity since PSP has been predominantly south of the HCS during each perihelion pass when it measures these low error values and may be summarised as $\Phi_{H}=\Phi_{open}=2.5^{+0.3}_{-0.6}$ nT AU². We note the $5\%$ error implied by the excess flux discussed in Section 3.2 is $0.125$ nT AU² and is sub-dominant to the statistical uncertainty.

Comparing the PFSS and in situ results, we see regardless of time evolution, distribution time interval, source surface height, or choice of magnetogram, the best estimate of open flux measured by PSP exceeds that expected by PFSS estimates. As is to be expected, lowering the source surface height increases the estimated flux. Continuing below 2.0 $R_{\odot}$ would eventually enhance the flux sufficiently to match the in situ data, but this is at the expense of realistic coronal hole structure (Linker et al., 2017). The PFSS estimates shown here do exhibit some limited time evolution over these two years of data but the fluctuations are in general much less than the difference (1) between different PFSS models (source surface height and magnetogram source) and (2) the difference between the model values and the in situ values. Overall, we observe that the open flux measured by PSP for the time interval from 2018-2020 (deep solar minimum) remains too high to be explained by standard PFSS estimates.