Global solar wind variations over the last four centuries

The most recent “grand minimum” of solar activity, the Maunder minimum (MM, 1650–1710), is of great interest both for understanding the solar dynamo and providing insight into possible future heliospheric conditions. Here, we use nearly 30 years of output from a data-constrained magnetohydrodynamic model of the solar corona to calibrate heliospheric reconstructions based solely on sunspot observations. Using these empirical relations, we produce the first quantitative estimate of global solar wind variations over the last 400 years. Relative to the modern era, the MM shows a factor 2 reduction in near-Earth heliospheric magnetic field strength and solar wind speed, and up to a factor 4 increase in solar wind Mach number. Thus solar wind energy input into the Earth’s magnetosphere was reduced, resulting in a more Jupiter-like system, in agreement with the dearth of auroral reports from the time. The global heliosphere was both smaller and more symmetric under MM conditions, which has implications for the interpretation of cosmogenic radionuclide data and resulting total solar irradiance estimates during grand minima.

Scientific RepoRts | 7:41548 | DOI: 10.1038/srep41548 found to be a function of only solar cycle phase. Thus OSF can be reconstructed by sunspot number alone. Such estimates of OSF have been well validated against geomagnetic reconstructions back to 1845 12 and against 10 Be, 14 C and 44 Ti in meteoritic material back to 1610 15,26 .
New OSF, in the form of closed loops, will enter the streamer belt. As these loops propagate further into the heliosphere, they will ultimately add to the open corona-hole flux. By optimising the time constant for conversion between streamer belt and coronal hole flux, Lockwood and Owens 27 reconstructed the fraction of OSF contained within the streamer belt. From this, they estimated the latitudinal streamer belt width (SBW) as SBW = sin −1 (F SB / OSF), where F SB is the magnetic flux contained in the streamer belt, back to 1617. The observed correspondence between slow wind and the streamer belt means SBW is expected to be closely related to the latitudinal extent of slow solar wind. Figure 1 summarises the pertinent results of Lockwood and Owens: panel a shows the two sunspot data composites considered; R G , the group sunspot number record 19 , and R C , a composite of group and international records with various calibrations/corrections applied 28 . Panel b shows the estimate of SBW from the sunspot records and their agreement with SBW estimates manually scaled from available eclipse images (detailed in the Figure caption and Appendix A). It can be seen while the two sunspot records differ considerably, the resulting SBW is not sensitive to the choice of sunspot record. In particular, the change in SBW between the modern era (e.g., 1960-present) and the MM is very similar for the two sunspot number records. Thus while there is currently much on-going work to recalibrate and correct the sunspot records [29][30][31] , the outcomes are not expected to greatly influence the results presented here.

Methods
SBW is now calibrated in terms of solar wind speed, V SW , using a combination of in situ spacecraft observations and the "Magnetohydrodynamics Around a Sphere" (MAS 32 ) global coronal model constrained by photospheric magnetic field observations. For a given Carrington rotation (CR), the MAS model extrapolates the photospheric field distribution outward to 30 solar radii (R S ), while self-consistently solving the plasma parameters on a non-uniform grid in polar coordinates, using the MHD equations and the vector potential A (where the magnetic field, B, is given by ∇ × A), such that ∇ .∇ × A = 0 (which ensures current continuity, ∇ . J = 0, is . Circles in (b) are scaled from the eclipse images listed, with sources, in Appendix A. Coloured points and lines relate to the selected example events for sunspot minimum/maximum shown along the top/bottom of the plots. Light blue dots are SBW from photographs and detailed paintings or lithographs made from photographic images; open circles are from descriptive texts and sketches, the black dot is from a Skylab coronagraph image and the grey points are from the Loucif and Koutchmy 52 catalogue.
conserved to within the model's numerical accuracy). All model data can be downloaded from http://www.predsci.com/mhdweb/home.php. The calibration period covers 1975-2013, the timespan for which both photospheric magnetograms (and hence MAS estimates) and SBW estimates are available. A combination of Kitt Peak, Wilcox Solar Observatory, Mount Wilson Solar Observatory, SOLIS and GONG data are used to minimse data gaps and provide the longest possible time sequence. V SW is extracted at 30 R S , the interface of the coronal and solar wind models.
Comparison of MAS with spacecraft data reveals MAS reproduces the observed solar wind structure, but generally underestimates the fast wind speed. A linear scaling of 1.04 and an offset of 78.5 km/s is determined to provide the best fit of MAS 27-day means (discussed below) to spacecraft V SW observations. This corrected MAS solar wind speed is referred to as V MAS and shown in Fig. 2. To later enable comparison with the annual sunspot-based SBW estimates, which provide information about the latitudinal solar wind structure, but no longitudinal information, "zonal means" (i.e., longitudinal or 27-day averages) are used. Data are further smoothed using a 1-year rolling (boxcar) means, to mimic the annual SBW estimates while retaining rapid latitudinal variations of the observing spacecraft. Figure 2a shows sunspot number and the latitude of the observing spacecraft, with the near-Earth OMNI data 11 in white and Ulysses 33 in red. Figure 2b shows in-ecliptic V SW observations from OMNI (black) and V MAS (blue) at Earth's heliographic latitude. There is only a very weak solar-cycle variation annual V SW . Indeed, to first order, the OMNI V SW can be approximated as 430 km/s with ~50 km/s variability. In that respect, V MAS is in good agreement. Observations made by the Ulysses spacecraft, shown in Fig. 2c, provide a better picture of global solar wind structure. Ulysses performed 3 "fast-latitude scans" of the Sun, sampling all solar latitude within approximately 1 year. For the two solar minimum fast-latitude scans, in approximately 1994-1996 and 2006-2008, the global solar wind structure is fast wind from the poles and a narrow slow wind band at the equator. This is well reproduced by MAS. The broader slow wind band during the second solar minimum pass (2006)(2007)(2008), partly the result of increased pseudostreamer occurrence 2 , is also captured by MAS. During the solar maximum fast-latitude scan, 2000-2002, the agreement is good in the northern hemisphere, but V MAS is higher than observed in the southern hemisphere 34 . In general, however, there is sufficient agreement between V MAS and observed V SW that we proceed in using V MAS to calibrate the sunspot-based reconstructions.
An example of V MAS for a single Carrington rotation (CR1756, which approximately spans December 1984, the late declining phase of solar cycle 21), is shown in Fig. 3. Panel a shows V MAS at 30 R S in heliographic coordinates. The data have been reduced in resolution to 90 and 45 equally-spaced grid points in longitude and sine latitude, respectively, to enable efficient manipulation of the data. The heliospheric current sheet (HCS) is shown as the thin white line. The dashed black line in panel b shows the zonal mean V MAS for CR1756. As is typical of non-solar maximum periods, the solar wind structure can be approximated as fast wind from the poles, with a roughly equatorial "belt" of slow wind, primarily centred on the HCS. Thus the zonal mean V SW in heliographic coordinates is a product of a number of factors: 1. The inclination of the slow wind belt to the solar rotation direction. Given the association of slow wind with the HCS, this is to first order equal to the magnetic dipole tilt. 2. The angular width of the slow solar wind belt about the HCS location. 3. The higher order "waviness" of the belt. The waviness of the HCS results from quadrapolar (and higher) order components of the magnetic field. 4. Sources of slow wind not directly associated with the HCS (e.g. pseudostreamers or "S-web" 35 ).
As is shown below, factor 1, the large-scale inclination of the solar wind speed structure to the solar equator, is a strong function of solar cycle phase and thus does not need to be calibrated in terms of SBW. Removal of inclination from V MAS estimates allows the true slow wind band width to be more readily estimated from zonal means. To this end, the V MAS solution at 30 R S is put through a coordinate transform to remove the large-scale inclination of the slow wind band. As will be shown, this is approximately equivalent to transforming to heliomagnetic coordinates. Specifically, we determine the coordinate system which maximises the difference in zonal-mean V MAS between the equator and poles. For CR1756, this requires inclining the North rotation pole up by 32 degrees (through a longitude of − 74 degrees). V MAS in this inclined coordinate system is shown in panel c. While the coordinate transform is determined purely from V MAS data, it can clearly be seen that it also reduces the overall magnetic dipole tilt, producing a HCS much closer to the equator (though shorter-scale corrugation is still present). Zonal mean V MAS in the inclined coordinate system (panel d) exhibits a narrower and deeper slow wind belt than in heliographic coordinates. Note that angular width of the zonal mean slow wind band in inclined coordinates still includes the waviness of the slow wind belt (factor 3) and the non-HCS sources (factor 4), which we assume can be calibrated in terms of SBW. Figure 4 shows zonal mean V MAS for all Carrington rotations in the period 1975-2013. Panel a shows data in heliographic coordinates, while panel c shows the inclined coordinate system, as described above. The angle  36 ) with small values at solar minimum, when the rotation and magnetic axes are aligned, and a saw-tooth increase peaking at solar maximum. As can be seen from panel c, inclined coordinates do generally result in a narrower slow wind band than heliographic coordinates. In particular, slow solar wind at mid-latitudes in heliographic coordinates during the declining phase of the solar cycle (e.g., 2002-2004) is largely absent in inclined coordinates, indicating that the broadened zonal-mean slow wind band is the result of increased tilt of the slow wind band to the solar rotation direction, not a broadening of the slow wind band itself. On the other hand, the increased latitudinal extent of the slow wind during the 2007-2009 solar minimum over with the previous 2 minima appears to be the result of both increased tilt and broadening of the band itself.
The next step is to characterise zonal mean V MAS structure for each CR by a reduced number of parameters. We use a simple functional form for the zonal mean V MAS , describing it by a maximum solar wind speed (V 0 ) with a sinusoidal dip, centred on the equator, of depth dV and angular width θ V 37 . As SBW contains no hemispheric information, hemispheric averages of zonal V MAS are used, shown as solid black lines in Fig. 3b and d. V 0 is taken to be the maximum value of the hemispheric-averaged zonal mean V MAS . dV is the difference between V 0 and the minimum value of the hemispheric-averaged zonal mean V MAS . θ V is fit as a free parameter. Examples of best fits are shown as the red lines in Fig. 3b and d. In heliographic coordinates, CR1756 yields V 0 = 757 km/s, θ V = 61.1 degrees and dV = 318 km/s. In inclined coordinates it yields V 0 = 757 km/s, θ V = 34.0 degrees and dV = 373 km/s. The black lines in panels a to c of Fig. 5 show time series of annual means of V 0 , θ V and dV (from data in inclined coordinates) over the period 1975-2013. In panel a, θ V exhibits a strong, saw-tooth-like, solar cycle variation, with an increase in θ V over the last three minima. The 1980 and 2000 maxima show a slower decline in θ V than the 1990 maximum. V 0 is constant at 757 km/s for the bulk of the interval, with short-duration (1-2 years) drops at solar maximum (panel b). This is a result of the prevalence of slow solar wind at all latitudes at solar maximum, both in heliographic and inclined coordinates. dV is only weakly ordered by the solar cycle (panel c).
Panels d to f of Fig. 5 show scatter plots of V 0 , θ V and dV, respectively, with the SBW estimates from the sunspot-based OSF model. Panel d shows that SBW is strongly ordered with θ V , as expected. The red line shows the 2 nd order polynomial fit, which is used to produce the θ V estimate from SBW (shown as the red line in panel a).

Results
The three relations shown in Fig. 5 are used to reconstruct annual zonal mean V SW in inclined coordinates on the basis of sunspot-based SBW estimates. The inclination angle is assumed to be a function only of solar-cycle phase, as shown by the red line in Fig. 4b. The resulting V SW reconstruction in heliographic coordinates is referred to as V RECON . Figure 2 shows a comparison of V RECON at Earth's heliographic latitude with both in situ V SW observations and V MAS . As expected, the exact details of in-ecliptic solar wind speed (panel b) are not captured, but V RECON produces the approximately the same mean and range of variability as observed. The global solar wind structure revealed by Ulysses (panel c) is well matched, with a notable exception: As with V MAS , faster wind is present over the south pole around 2000 than was observed.
In order to investigate the global structure of V RECON further, it is compared directly with V MAS in Fig. 6. Panel a shows annual means of V RECON (top) and V MAS (bottom) as a function of heliographic latitude and time. In general, the agreement is very good, both in the solar cycle variation and the cycle-to-cycle variations, particularly the broader slow solar wind band in the 2007-2010 minimum compared with the previous minima. There are, however, a number of differences. Firstly, the slow wind band broadening at solar maximum (approximately 1980, 1990 and 2000) occurs around 0.5 to 1 year earlier in V RECON than V MAS . Secondly, during the declining phase of the solar cycle, particularly 1992-1994, the slow wind band narrows more rapidly in V RECON than in V MAS . Figure 6b shows the global mean solar wind speed from a latitudinal integration of the data presented in Fig. 6a (allowing for the reduced element area towards the poles). Again, the general agreement is very good (r = 0.87, N = 39). While the approximately 0.5-to 1-year time lag is present for much of 1975-1995, it is not obvious for 1995-2013, so no correction is made in the remainder of this study.
Having demonstrated that sunspot-based SBW estimates can be used to make a first-order reconstruction of the global solar wind structure of the period 1975-2013, we now extend the reconstruction back through the 17 th century. Because of "spin up" of the OSF reconstruction and associated SBW estimate, the reconstructions presented here starts in 1617 (rather than 1611), as results from this date onwards do not depend on chosen initial conditions. Figure 7a shows the (unsigned) open solar flux, OSF, estimated from sunspot observations 36 , from which the SBW (Fig. 1b) is subsequently computed 27 . The Maunder minimum (MM, approximately 1650 to 1710) and, to   also different during the MM, with nominal values of around 350 km/s and short-duration upward spikes to approximately 600 km/s (i.e., there is a change from a solar minimum to solar maximum dominated heliosphere). The amplitude of the solar cycle variation of in-ecliptic V RECON is also enhanced at this time, with a typical variation from 270 to 430 km/s. It's likely the change in behaviour during the MM is a result of the switch from source-driven (i.e., sunspot) variations in OSF, to loss-driven variations, which has been shown to alter phase of the solar cycle variation in OSF during the MM 39 . Perhaps most interesting is the period 1617-1650, which contains small amplitude sunspot cycles just prior to the MM (see also Fig. 1a). Here, the global V RECON drops to its lowest values, despite the OSF being enhanced above MM values. This is discussed further in the next section.
These new estimates of the global V SW enable a number of further heliospheric parameters to be calculated on both a theoretical basis and by assuming empirical relations between parameters established from spacecraft data are valid over the prior four centuries. The black line in Fig. 8a shows the in-ecliptic magnetic field intensity, |B|, at 1 AU using the radial magnetic field magnitude determined from OSF (Fig. 7a) and assuming a Parker spiral angle consistent with the in-ecliptic solar wind speed. As expected, the variation closely follows the OSF, though the fractional variation between the MM and modern era is lower in |B| than in OSF, as the reduced V SW during the MM increases the winding of the Parker spiral, adding to |B|. The pink-shaded area in Fig. 8a shows the ratio of ecliptic |B| to polar |B| at 1 AU, again assuming an ideal Parker spiral. For 1730-2013, this ratio is fairly constant at around 1.4:1, but during the MM the slower solar wind speeds cause an over-winding of the ecliptic HMF and hence an increase in the ecliptic |B|, resulting in a higher ratio of nearly 2:1.
Proton temperature, T P , can be computed from V SW using empirical relations 40 , assuming they are unchanged during grand minima. Similarly, the assumption of constant mass flux in the solar wind 41 , can be used to estimate proton density, n P . Taking 27-day means of the entire OMNI, HELIOS and Ulysses datasets, we find a good linear anti-correlation (r = − 0.69 for N = 567) between observed n p (scaled to a radial distance of R 1 = 1 AU. i.e., n P /R 1 2 ) and V SW . As the data are normally distributed, we use an ordinary least squares fit to to determine n P (R 1 ) = − 0.0146 V SW [km/s] + 13.7 cm −3 . Using this relation, Fig. 8b shows the reconstructed annual n P at 1 AU, in the same format as panel a. As expected from the linearity of the n P -V SW relation used, reconstructed n P shows the same qualitative features as V SW : relatively little in-ecliptic variation for 1730-2013, but a marked change during the MM, in this case an increase of ~40%. The n P reconstruction allow us to quantitatively investigate dynamic pressure (P D , which varies as n P V SW 2 ), Alfven speed (V A ) and Alfven Mach number (M A ) variations over 1617-2013. These are shown by panels c, d and e, respectively. During the MM, the in-ecliptic P D is reduced by more than a factor 2 over the modern era. The in-ecliptic V A variation closely follows the B variation and hence is suppressed during the MM. Finally, the in-ecliptic M A increases by up to a factor 4 during the MM compared to the modern era.

Discussion and Conclusions
In this study, we have used sunspot-based reconstructions of open solar flux (OSF) and streamer belt width (SBW) to provide the first estimate of global solar wind conditions over the last 400 years. Calibration was performed using 35 years of output from a photospheric magnetic field-constrained global magnetohydrodynamic (MHD) model of the corona and solar wind. The reconstructed annual solar wind speed (V SW ) agrees well with both in situ spacecraft observations and the global solar wind speed inferred from the MHD model. Assuming these relations hold prior to the space age, annual solar wind speed reconstructions are extended back to 1617. A dipolar solar wind structure (fast/slow wind from the poles/equator) has dominated for much of the last 300 years, with short intervals of slow wind extending to all latitudes around solar maximum. Paradoxically, during the Maunder minimum (MM, 1650-1710), solar maximum-like conditions dominated, with fast wind over the poles for shorter durations and over a reduced latitudinal range than for the last 300 years. This, however, agrees well with the "most probable" MM coronal magnetic field configuration considered by Riley, et al. 16 . In-ecliptic, the annual zonal-mean solar wind speed varies very little (~70 km/s variations about a mean value of 420 km/s) back to circa 1715, but drops to a minimum of approximately 275 km/s during and prior to the MM, in approximate agreement with estimates by Cliver and von Steiger 42 . At the times of extremely low in-ecliptic V SW , polar V SW is similarly reduced, suggesting any longitudinal structure in in-ecliptic V SW is likely to be very weak. The longest periods of sustained slow wind actually occur prior to the Maunder minimum, when there were weak sunspot cycles. This appears to be the result of critical balance between OSF source and loss terms 39 which results in maximum loss of coronal-hole flux compared to production of streamer-belt flux, maximising SBW.
From the OSF and V SW reconstructions, it is possible to further estimate: magnetic field intensity (B) by assuming a Parker spiral heliospheric magnetic field; solar wind proton density (n P ) by assuming constant mass flux; and solar wind proton temperature (T P ) using empirical relations. These data are made available as Supplementary Material. From these parameters, we further derive values for the solar wind dynamic pressure (P DYN ), Alfven speed (V A ) and Alfven Mach number (M A ). The full implications of these results will be investigated in detail in a number of follow-on studies. Here, we discuss the first-order implications. While the MM specifically is discussed below, it is expected that these findings could be broadly applied to any past or future 18 grand minima of solar activity.
Firstly, we consider the terrestrial implications. Space weather is primarily the result of rapid changes in the space environment, rather than annual variations reconstructed in this study. Nevertheless, the equilibrium state of the terrestrial magnetospheric system is expected to be very different under MM than modern conditions. This, in turn, will mean a different response to a space weather driver, such as a fast coronal mass ejection. Future work will use a global MHD model of the coupled magnetosphere-ionosphere-thermosphere system to quantitatively investigate this. But even without a numerical model it is possible to draw some qualitative conclusions. The lower P DYN during the MM would increase the average stand-off distance of the dayside magnetopause 43 . The width of the far magnetospheric tail, however, is controlled by the solar wind static pressure, P STA = n p kT SW + B 2 /(2μ o ). As the higher n p and T P have a larger effect on P STA than the reduction in B, the tail would, on average, be somewhat thinner during the MM than in modern times. Thus the magnetosphere would have presented a smaller cross-sectional area to the solar wind, reducing the electric field placed across it by the solar wind and the total solar wind energy that it intersects. A reduction in V SW and B would mean a reduction in the solar wind electric field, which in turn would combine with the smaller diameter of the magnetosphere to reduce the trans-polar cap potential and polar cap area 44 . Thus the Earth's magnetosphere would have been somewhat more Jupiter-like, with the part driven by solar wind-driven convection smaller in extent, and the part driven by internal dynamics and co-rotation larger in volume. In addition to an expected reduction in both recurrent and non-recurrent geomagnetic storms during the MM, the expected poleward motion of the nominal auroral oval position may further help explain the dearth of auroral reports from that period for all but the most northerly locations 15 . Beyond the magnetopause, the enhanced M A suggests that the bow shock strength would be enhanced, resulting in more efficient energetic particle acceleration, while the bow shock stand-off distance would be increased on average, resulting in a thicker magnetosheath 45 .
Secondly, we consider the implications for the global heliosphere. Again, a future study will use the reconstructed solar wind parameters with a MHD model of the global heliosphere, but here we consider the first-order implications. Most obviously, a drop in P DYN will result in an overall smaller heliosphere, though the contribution of pick-up ions to the total solar wind momentum budget 46 means the P DYN decrease at large heliocentric distances will be lower than the factor 2 between modern and MM 1-AU values. Any calculation of the heliopause distance under grand solar minima conditions will also need to account for the change in pick-up ion acceleration under the MM reduction in B, particularly out of the ecliptic plane. The shape of the heliosphere is also likely change under MM conditions. For the modern era, P DYN has been ~2-3 higher at the poles than the solar equator 47 , which results in latitudinal asymmetry in the heliopause stand-off distance and termination shock location 48 . During much of the MM, however, P DYN becomes almost uniform with latitude for a greater period of time, suggesting a more spherical heliosphere and termination shock.
In turn, there will also be a number of implications for cosmic ray intensity in near-Earth space, with potential knock-on effects for long-term heliospheric reconstructions on the basis of cosmogenic radionuclide records in ice cores and tree trunks 23,49,50 . The relative abundance of radioisotopes such as 10 Be and 14 C can be used to determine the effective shielding of heliosphere from the interstellar cosmic ray spectrum, referred to as the heliospheric modulation potential. Interpreting the modulation potential in terms of heliospheric parameters, such as OSF, necessitates a number of assumptions about the size of the heliosphere, the solar wind speed and the scaling of cosmic ray scattering centers with the HMF intensity 20 . During grand minima, all of these properties will change, to some degree. As already discussed, we expect a smaller heliosphere, with lower and more symmetric solar wind speeds. The lack of latitudinal solar wind speed structure suggests reduced corotating interaction region formation and hence reduced cosmic ray scattering (even for the same OSF). Furthermore, we note that enhanced V A during the MM would increase the termination shock strength and may affect the efficiency of anomalous cosmic ray acceleration 46 . While the effect of changing size/shape of the heliosphere is expected to be small on GeV (and greater) energy particles which are largely responsible for cosmogenic isotope production, and hence radionuclide reconstructions of the heliosphere and total solar irradiance 51 , it needs to be fully quantified via a galactic cosmic ray transport model and a cosmogenic isotope production model.