Reconstruction of Helio-latitudinal Structure of the Solar Wind Proton Speed and Density

The modeling of the heliosphere requires continuous three-dimensional solar wind data. The in-situ out-of-ecliptic measurements are very rare, so that other methods of solar wind detection are needed. We use the remote-sensing data of the solar wind speed from observations of interplanetary scintillation (IPS) to reconstruct spatial and temporal structures of the solar wind proton speed from 1985 to 2013. We developed a method of filling the data gaps in the IPS observations to obtain continuous and homogeneous solar wind speed records. We also present a method to retrieve the solar wind density from the solar wind speed, utilizing the invariance of the solar wind dynamic pressure and energy flux with latitude. To construct the synoptic maps of the solar wind speed we use the decomposition into spherical harmonics of each of the Carrington rotation map. To fill the gaps in time we apply the singular spectrum analysis to the time series of the coefficients of spherical harmonics. We obtained helio-latitudinal profiles of the solar wind proton speed and density over almost three recent solar cycles. The accuracy in the reconstruction is, due to computational limitations, about 20%. The proposed methods allow us to improve the spatial and temporal resolution of the model of the solar wind parameters presented in our previous paper (Sok\'o{\l} et al. 2013) and give a better insight into the time variations of the solar wind structure. Additionally, the solar wind density is reconstructed more accurately and it fits better to the in-situ measurements from Ulysses.


Introduction
The heliosphere is a cavity in the local interstellar medium created by the solar wind (SW thereafter), the plasma that originates from the Sun. For modeling the structure of the heliosphere and ionization conditions of the interstellar neutral gas inside the heliosphere, a three-dimensional (3D) structure of the SW parameters has to be known. The required dimensions are: time, latitude, and distance to the Sun. In-situ measurements of the SW are available from the mid-1960s, but mostly in the ecliptic plane and at 1 AU (for recent review, see Bzowski et al., 2013). In-situ measurements of the SW out of the ecliptic plane were only carried out by Ulysses from 1990 to 2009 (e.g. Bame et al., 1992;McComas, Gosling, and Skoug, 2000;. However, they were point measurements, and a time series of these observations is in reality a convolution of variabilities in time, distance, and helio-latitude. The information on continuous variations in the global (particularly spatial) structure of the SW is still missing in the in-situ measurements. On the other hand, the SW is also observed by remote-sensing methods, by interplanetary scintillation (IPS thereafter; e.g. Hewish, Scott, and Wills, 1964;Houminer, 1971;Coles et al., 1980;Kojima et al., 2004) and by Lyman-α helioglow (e.g. Lallement, Bertaux, and Kurt, 1985;Bzowski et al., 2003) observed by the Solar Wind ANisotropy (SWAN) instrument onboard the SOlar and Heliospheric Observatory (SOHO). Sokó l et al. (2013) used Ulysses measurements, the OMNI in-ecliptic SW database (King and Papitashvili, 2005), and the SW speed derived from IPS observations conducted at the Solar-Terrestrial Environment Laboratory (STEL), Nagoya University  to reconstruct the SW speed and density evolution from 1990 to 2011. They constructed SW structures in a grid of 1 year spacing in time and 10 • in helio-latitude. The choice of the grid steps was based on the available time coverage of the SW speed obtained from the IPS data, which have systematic yearly breaks, typically from December to March each year . Moreover, all heliolatitudes are not continuously observed; see the example SW speed map from IPS data in Figure 1. A more detailed description of the SW speed derived from IPS observations is given in Section 2.
In this paper we eliminate the limits in the analysis of Sokó l et al. (2013) and fill the gaps and breaks in the SW speed data from IPS observations to obtain maps of the SW speed that are continuous in time and helio-latitude. We propose a method of filling the temporal and spatial gaps by adopting the techniques used with success in the analysis of geophysical, solar, and in-ecliptic SW data (e.g. Kondrashov and Ghil, 2006;Kondrashov, Shprits, and Ghil, 2010;Dudok de Wit, 2011;Kondrashov et al., 2014). We decompose the SW speed maps from IPS observations into a set of spherical harmonics and perform a singular spectrum analysis (SSA; based on Ghil et al., 2002) of the time series of the coefficients of spherical harmonics to fill the data missing in time. Sokó l et al. (2013) proposed a method to calculate the SW proton density based on the SW speed using an empirical relation between the proton speed and density retrieved from Ulysses measurements. We revise this method by comparing its results with the SW proton density calculated using the SW heliolatitudinal invariants, which are the SW dynamic pressure (McComas et al., 2008) and the SW energy flux (Le Chat, Issautier, and Meyer-Vernet, 2012).

Solar Wind Speed from IPS
IPS is the phenomenon due to diffraction patterns on an observer's plane produced by interference of radio waves coming from a remote compact radio source (e.g. quasar), scattered on electron density irregularities (fluctuations) in the SW (e.g. Hewish, Scott, and Wills, 1964;Houminer, 1971;Coles and Maagoe, 1972;Kakinuma, 1977;Coles and Kaufman, 1978;Coles et al., 1980;Kojima and Kakinuma, 1990;Jackson et al., 1997Jackson et al., , 1998Jackson et al., , 2003. IPS observations are ground-based, line-of-sight (LOS) integrated measurements, which limit the availability of the data in the case of adverse weather conditions or insufficient elevation of the Sun above the horizon.
In our analysis we use the SW speed data derived from the IPS observations of STEL . The IPS data were deconvolved using the computer-assisted tomography method (CAT; Asai et al., 1998;Jackson et al., 1998;Kojima et al., 1998) to provide the reconstructed SW speed. This analysis provides Carrington rotation 1 (CR) maps of the SW speed, which are available since the 1980s. In this method, spatial and temporal gaps are inevitable in the derived SW maps. Periodic gaps appear as a break in the SW speed because all or a portion of the STEL IPS system is closed during the winter months. Additionally, spatial gaps are present due to the reduction in valid observations at higher helio-latitudes, particularly in southern polar regions in winter. Figure 2 shows the time breaks from 1985 to 2013, and the left-hand panels in Figure 3 show spatial gaps in the CR maps. The SW speed retrieved from IPS observations does not exceed 800 km s −1 , because the SW model used in the CAT analysis has an upper bound at 800 km s −1 established based on the Ulysses measurements of the fast SW.
The retrieval of the SW speed from the IPS observations is possible owing to a relation between the electron density fluctuation level and the SW speed as a function of heliocentric distance  as where r is the heliocentric distance, N e is the electron density, V is the SW speed, and index γ (r) varies with solar distance and must be determined experimentally (e.g. Coles et al., 1995;Manoharan, 1993;Tokumaru, Kojima, and Fujiki, 2012). The relation used to connect the fluctuation level with the SW speed had been established before the secular changes in the SW data were observed both in the ecliptic plane and out of the ecliptic by Ulysses (see the discussion in Sokó l et al., 2013). This relation may be a possible source of systematic bias in the reconstruction of the SW speed. The relation between ∆N e and V was used to analyze the IPS data before 1997 using the CAT method. After 1997 the quantity called the g-value (Gapper et al., 1982) was introduced to the IPS analysis performed by STEL (Tokumaru et al., 2000). It represents the relative variation of the scintillation strength with respect to the mean level: g = ∆I/∆I (R), where ∆I and ∆I (R) are the observed instantaneous and average scintillation levels, respectively. ∆I is computed from the observed power spectrum of the fluctuations P (f ) as ∆I (R) = aR −1.5 is a function of the solar offset distance R, by assuming a spherically symmetric distribution of the SW density fluctuations. The value of a is derived by fitting the observations of ∆I at low latitudes (see more details in Tokumaru, Kojima, and Fujiki, 2012). Since 1997, therefore, two alternative reconstructions of the SW speed have been available from STEL: one obtained using the g-value and ∆N e , and the other one using only the relation of Equation (1). Some revision of the ∆N e ∼ V relation for the inner heliosphere was published e.g. by Hick and Jackson (2004) and Jackson et al. (2010); they also found little variation in the level of g-values with solar distance beyond that imposed by an R −2.0 radial ∆N e fall-off.
As shown in Asai et al. (1998) and Tokumaru, Kojima, and Fujiki (2012), the SW speed derived from the relation ∆N e ∼ V −0.5 fits well the results of analysis with additional g-value data, except for the years of solar maximum. Therefore, the results from the CAT analysis using only the ∆N e data can be used to discuss long-term changes in the SW over solar cycles.
In this paper we focus only on the SW speed data retrieved from the IPS observations of STEL. In Section 5 we present a simple method of calculation of the SW density based on the SW speed and helio-latitudinal invariants.

Methods
The IPS data are extensively used (e.g. recent works by Fujiki et al., 2014;Kim et al., 2014;Jackson et al., 2015) and widely incorporated into the MHD codes for the modeling of the near-Sun environment as well as the far regions of the heliosphere. In our analysis we want to improve the simple model proposed by Sokó l et al. (2013), which reconstructed the structure of the SW proton speed and density at 1 AU averaged on a time scale of ≈ 1 year, by increasing the temporal and spatial resolution.
As the starting point, we use the SW speed from the CAT analysis of STEL IPS observations, which are yearly grouped CR maps reduced to the distance of 2.5 solar radii, with a 1 • resolution in longitude and latitude, from 1985 to 2013. There was a one-year break in 2010, when the number of observations was insufficient for a complete analysis. For the period of 1985 to 2007 we adopt the maps prepared using the CAT analysis based on ∆N e ∼ V only, and for 2008-2013 we use the data prepared with the additional information of the g-value.
There were 389 CRs from the beginning of 1985 to the end of 2013. During this period we have 298 CR maps of the SW speed from IPS observations, among which 88 have more than 50% of surface empty. This means that we need to fill the spatial gaps in 210 maps and must fully reconstruct 179 CR maps. All calculations were done using Mathematica 2 (Wolfram Research). The details of the algorithms of the decomposition into spherical harmonics and singular spectrum analysis are described in Appendices A and B, respectively.

Reconstruction of the Helio-latitudinal Structure by Decomposition into Spherical Harmonics
The IPS observations give information about the SW speed which originates from the surface of the Sun. Therefore the spherical frame is the most obvious coordinate system to use in our analysis. In the first step of the process of filling the SW speed structure we decompose the original maps into spherical harmonics which, due to their characteristics of orthogonality on a sphere, are a common tool used in similar studies. Because of the computational limitation (the analysis was performed on a personal computer) we reduce the spatial resolution to 3 • × 3 • in longitude and latitude. Next we calculate the fraction of empty cells (3 • × 3 • in size) per map, and pick up the CR maps with less than 50% data missing. The SW speed maps obtained from the IPS observations were reconstructed by the decomposition into spherical harmonics, where the spherical harmonics Y ℓm (θ, φ) are defined in Equations (6) and (7) as real functions in Appendix A, W ℓm are the expansion coefficients, θ ∈ (0, π) is the colatitude, and φ ∈ (0, 2π) is the longitude. We set the limit ℓ ≤ 12, which allows us to reconstruct only structures larger than approximately 15 • . Each CR map of the SW speed was decomposed separately. The reconstruction of all maps for a given year was done with the information from the adjacent CR maps, but without information from the adjacent years. This approach is justified by the distribution of a few month breaks in the IPS data at the end and beginning of subsequent calendar years. Also, it is less probable that the SW source structures in the solar corona last longer than a few full rotations.
For each year we select the CR map with the smallest fraction of gaps (CR * 0 ). The best possible case is the period without gaps in longitude and latitude, e.g., CR 1938 in 1998, shown in the top-left panel of Figure 3. If there is a single map with the minimal fraction of gaps for the given year, then we use it as the start map for the analysis of this year. But if there are multiple maps with high coverage and they are not the consecutive CRs, and additionally each map has at least one non-empty cell in longitude for a given latitudinal 3 • -bin, then we set up the initial conditions for the neighboring maps using the W ℓm coefficients retrieved from the nearest start map. For most of the analyzed years, there are two distinct CR maps that can be used to set the initial conditions for the decomposition (e.g. CR 1 0 and CR 2 0 ). When we start the decomposition of the CR maps for a given year, we fill the empty 3 • × 3 • cells in the first step. For the start maps with the smallest fraction of gaps, we fill the gaps by an average value for each 3 • latitudinal bin using the information from all longitudes. The gaps in the next neighboring map are filled at the first step of the iterations by using the spherical harmonic decomposition of either the start map for a given year (CR * 0 ), if it is the nearest, or by the nearest completed map. The same process is applied to all CR-maps for the given year. This method makes effective use of long-lived, large and unambiguous structures at the source surface of the SW.
We calculate the spherical harmonic coefficients W ℓm as the scalar product of the filled map with spherical harmonics with appropriate weight function on the sphere given by Equation (6) in Appendix A. For each map the procedure is repeated until the mean relative difference of speed values for the initially empty cells (which were the gaps) stop to change by more than 0.01 or start to increase (see Appendix A for more details).
The decomposition into spherical harmonics was carried out to fill the spatial gaps for all the CR maps analyzed in this study. As shown in Figure 3, this method works very well compared to the original data for all kinds of spatial distribution of gaps both during the solar minimum and solar maximum. The left-hand column in Figure 3 presents various completeness levels of the original IPS SW speed CR maps (i.e., complete, with gaps only on one hemisphere, with gaps distributed only on one half of the map in longitude, and with a random gap distribution), and the right-hand column illustrates the final reconstruction by decomposition into spherical harmonics.

Filling in the Temporal Breaks by Singular Spectrum Analysis
Having the W ℓm coefficients for each CR map, we have filled the empty cells in maps, but we still have the breaks in time. As shown in Figure 2, the time breaks occur mostly at the beginning and end of the year. We have also one longer gap in 2010. We treat the coefficients W ℓm of spherical harmonic decomposition as a time series and we use the singular spectrum analysis (SSA) to fill the time breaks. We only focus on the monthly average 3 of the helio-latitudinal profiles of the SW. Therefore, we reproduce the maps that are averaged over longitude, namely we only retain coefficients W ℓm with m = 0 4 .
For the SSA to reconstruct the solar cycle (SC) variations of the SW, 2.5 SCs may not be sufficient to provide a valid reconstruction of this periodicity. Thus, we decided to use the solar F10.7 radio flux 5 (Tapping, 2013; see Figure 4) as an indicator of the SC variations because this record is being measured longer than three SCs. We correlate and normalize the W ℓm coefficients with the F10.7 flux, smoothed by a 13-CR running average as illustrated in Figure 4 (see Appendix B for more details).
We use the SSA algorithm proposed by Ghil et al. (2002). The window width M , which determines the longest periodicity captured by SSA, is in our analysis adopted as constant (more details in Appendix B). We search for the best representation iteratively up to the moment when the root-mean-square difference between the computed coefficients and the known ones is smaller than 0.01. The procedure is repeated for each ℓ of W ℓm coefficients separately. The results are shown in Figure 5, which presents the time series for some example numbers ℓ. The values found for the missing elements in the time series of coefficients fit very well the general shapes of the time series and trace their variations in time.
3 By monthly we actually mean the average over one CR period.
4 Note that for m = 0, the longitudinal average of the spherical harmonics vanishes.

Results of Reconstruction of Gaps in the Solar Wind IPS Data
The procedures described in the previous sections allow us to fill the spatial and temporal gaps in the SW speed maps from 1985 to 2013. Figure 6 presents the helio-latitudinal profiles of the SW speed for all the CR periods under study. Shown in gray are the profiles by filling the spatial gaps, and in red are shown the profiles after filling the time gaps. The reconstructed profiles trace the general shape of the latitudinal dependence of the SW speed very well for all phases of solar activity, from almost flat profiles during the solar maximum (e.g., 1990 and  2002) to a bi-modal structure close to the solar minimum (e.g., 1998, 2006). The complete map of the 3D structure of the SW proton speed at 1 AU is presented in Figure 7.

Solar Wind Density
The heliosphere is shaped by the SW ram pressure, which is proportional to the product of the SW density and squared speed. With the time-and latitudinal variations of the structure of the SW speed available, one needs information on these variations in the SW density for complete modeling. The SW density in Figure 8. Relation between daily SW proton density at 1 AU and speed in the ecliptic plane (source: OMNI dataset) shown as a 2D histogram (100×100 bins) with the count scale illustrated by the side bar.
the ecliptic plane has been measured in-situ by a fleet of spacecraft since mid-1960s, and the data are now available in a compiled and meticulously curated form; the OMNI database 6 ( King and Papitashvili, 2005). The out-of-ecliptic insitu observations are solely available from Ulysses, but only for years 1990-2009. The out-of-ecliptic remote-sensing observations of the heliospheric Lyman-α glow from SOHO/SWAN can provide the required information about the structure of the SW density after complex modeling (see a sketch of the algorithm of reconstruction of the SW density from the SWAN and SW speed from IPS data in Bzowski et al. (2013)), but this procedure is not free from yet unanswered questions (Katushkina et al., 2013). Also the g-level values from IPS observations can serve as a proxy to the out-of-ecliptic SW density (e.g. Houminer andHewish, 1972, 1974;Tappin, 1986;Hick and Jackson, 2004;Jackson and Hick, 2004). There exists a nonlinear relationship between the bulk density and the IPS scintillation level that has been widely and successfully used to study the corotating structures and coronal mass ejections (e.g. Jackson et al., 1998;Tokumaru et al., 2007;Bisi et al., 2009Bisi et al., , 2010Fujiki et al., 2014; a comprehensive review by Jackson et al., 2011), but in this study we focus only on the use of the SW speed determined from IPS observations. As is already known, there is no clear correlation between the SW proton speed and density in the ecliptic plane at 1 AU based on the OMNI dataset ( Figure 8). It seems that the crescent-like relation may be expected, but the shape of the crescent is too wide to provide a clearly defined correlation function.
We propose a simple approach with the use of the SW invariants to determine the SW density from the SW speed. Based on the Ulysses measurements Figure 9. The CR-averaged SW dynamic pressure (blue; Equation (3)) and energy flux (green; Equation (4)) in the in-ecliptic plane at 1 AU based on OMNI. The latter was multiplied by 1000 to roughly match the two scales.
two alternative SW quasi-invariants were inferred: the SW dynamic pressure (McComas et al., 2008) and the SW energy flux (Le Chat, Issautier, and Meyer-Vernet, 2012) where n p is the SW proton density, V is the SW proton speed, m p is the proton mass, ξ α is the abundance of α-particles (≈ 4% after Kasper et al. (2012)), m α is the mass of the α-particles, M ⊙ is the mass of the Sun, R ⊙ is the solar radius, and G is the gravitational constant. As shown in Figure 1 of Le Chat, Issautier, and Meyer-Vernet (2012) and in Figure 3 of McComas et al. (2008) for the Ulysses measurements, both the SW energy flux and the dynamic pressure are almost constant, independent of latitude. Figure 9 presents two quantities as a function of time, calculated for the in-ecliptic SW data at 1 AU from the OMNI database. Both quantities show very similar variations as a function of time and they can serve as the SW helio-latitudinal invariants to calculate the out-of-ecliptic structure of the SW density from the latitudinal variations of the SW speed as is discussed in Appendix B of McComas et al. (2014).
In the analysis below we calculate the SW proton density as a function of helio-latitude based on the SW proton speed and the SW invariants in latitude. We calculate the SW dynamic pressure from Equation (3) and the SW energy flux from Equation (4) for the in-ecliptic SW at 1 AU from the OMNI database as averages of the daily data over each CR period for the time ranges in question. Next, we calculate the density profiles (n p ) as a function of helio-latitude based on the SW speed profiles presented in Figure 6 and the invariant calculated from the in-ecliptic measurements (Figure 9). Results for the corresponding years are presented in Figure 10. Because the SW energy flux and dynamic pressure are invariant in latitude, the difference between the ecliptic plane and the solar equatorial plane does not matter. The invariance of the SW dynamic pressure and energy flux derived from Ulysses is also seen in the in-ecliptic values measured by spacecraft collected by OMNI, as both datasets agree quite well for the ecliptic values, as shown by Sokó l et al. (2013) in their Figures 11 and 12, over the yearly and hourly time resolution.
The complete map of the 3D structure of the SW proton density at 1 AU is presented in Figure 11. The approach we adopted has led to the confirmation of a general dichotomy of the SW, i.e., dense and slow wind vs. low density and fast-wind, as is seen in Figures 7 and 11.

Comparison with OMNI and Ulysses
SW measurements made in-situ give the reference for all studies of the SW structure. The in-ecliptic measurements at 1 AU have been gathered by various spacecraft since 1960s and are compiled in the OMNI database. The only available in-situ measurements out of the ecliptic plane come from the unique Ulysses mission, but as it samples a single point at a time in a wide range of radial distance (from ≈1.2 to ≈5 AU), latitude (≈ ±80 • ), and time (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), they cannot answer all questions about the latitudinal structure of the Figure 11. Helio-latitude vs. time maps of the SW density at 1 AU calculated from the SW speed and SW energy flux (Equation (4)).
SW and its evolution over the SC. Nevertheless, to be credible, all attempts of reconstruction of the SW must give results in agreement with these two datasets.
In the first step, we compared the results of our analysis with the in-ecliptic SW parameters at 1 AU from the OMNI database. We compared the CRaveraged values of the SW proton speed and density at 1 AU with those obtained from our model. The comparison for the speed is shown in the top-left panel of Figure 12; the blue line is for the time series from OMNI and the orange line is for the time series from the SW reconstructed from the IPS data. The agreement between the two is very good until the end of 2008; afterwards a systematic deviation appears, with the lower values from OMNI. Before this time, which coincides with the deepest solar minimum since the beginning of the space-age, the SW reconstruction traces the variations in the in-ecliptic SW very well. The middle left panel of Figure 12 presents the ratio between the SW speed from OMNI to the SW speed reconstructed in this analysis. For almost all years, the ratio is constant with a mild systematic upward trend. The source of this trend may be in the assumed relation between the density fluctuations and the SW speed from Equation (1) and may not be related to the adopted version of analysis of the IPS data (either ∆N e or ∆N e and g-value relation). This difference after 2008 suggests that either the ∆N e ∼ V relation needs verification considering the secular changes in the SW or, although highly unlikely, the OMNI data are in error.
The right-hand column of Figure 12 presents the results for the SW proton density adjusted to 1 AU from the Sun. The OMNI values (blue line) are compared with the density calculated using the SW latitudinal invariants; the SW dynamic pressure (green line, Equation (3)) and the SW energy flux (gray line, Equation (4)). Similar to the case of the SW speed, the agreement between the two is very good up to 2009. Then the discrepancy increases, with the density from our analysis underestimated. This difference is understandable, because in our method of calculation of the SW density, any inaccuracy in the SW speed affects the derived density. As shown by the histograms of differences of Figure 12. CR-averaged SW proton speed (left-hand column) and density (right-hand column) derived in our model are compared with the values from the OMNI database (blue lines in the top-row panels) in the ecliptic plane at 1 AU. In the top-left panel the orange line shows the SW speed reconstructed in this study; in the top-right panel the gray and green lines show the SW density calculated from the SW energy flux or dynamic pressure, respectively. The middle-row panels show the ratios of the SW quantities obtained in our analysis with respect to the OMNI data. The bottom panels show the histograms of the differences between our model and the OMNI values; the color coding is the same as in the middle-row panels.
the absolute values in the bottom panels in Figure 12, the mean accuracy of reconstruction is about 50 km s −1 for the SW speed and 1.5 cm −3 for the SW density.
We also compared the results of the present method with the out-of-ecliptic measurements from Ulysses 7 and our previous model   (Figure 13). The Ulysses data are scaled to the distance of 1 AU. The present model traces the short time-scale variations of the SW speed much better than the previous model by Sokó l et al. (2013). However, for some years it does not precisely reproduce the SW speeds measured by Ulysses, where the speeds faster by more than 50 km s −1 (1993-1994) or slower (1999-2001) are measured. The comparison of SW density looks much better; the present model reproduces 7 Ulysses Solar Wind Observations Over the Poles of the Sun (SWOOPS) experiment, data source: ftp://nssdcftp.gsfc.nasa.gov/spacecraft data/ulysses/plasma/swoops/ion/. the SW density measured by Ulysses with higher precision except for the solar maximum of SC 23, when Ulysses sampled higher values.
The results of the reconstruction of the SW structure proposed by Sokó l et al. (2013) and in this paper are, for both speed and density, similar in general, but the present one is much better in detail. The blue line in the top-left panel of Figure 13 oscillates around the black line, which means that the present model gives more precise information on time variations of the SW speed. Furthermore, the new method based on the SW helio-latitudinal invariants reproduces the SW density better than the model proposed by Sokó l et al. (2013) which was based on the relation between the SW speed and density measured by Ulysses during its fast latitudinal scans (i.e., part of the orbit closest to the Sun). The higher resolution in latitude and in time in the present model has provided us with a better reconstruction of the monthly variations of the SW data, compared to the smooth time series in the yearly-averaged series from Sokó l et al.  Figure 13 where the present model fits much better the in-situ measurements from Ulysses.
The previous  and the present models of reconstruction of the SW structure, based on the IPS observations, indicate deviations from the measured Ulysses values at the same places (in time and latitude). Since both models are based on the same SW speed data from IPS observations, these differences may be real; the Ulysses data are only point measurements while the SW IPS data analysis gives information from a broader region of the sky. On the other hand, these differences may also point out that the IPS observations underestimate the SW speed at higher latitudes during the years with high solar activity.
Our final results for the SW speed and density, when compared to the insitu values from in-ecliptic and out-of-ecliptic spacecraft, show that typically the reconstructed values differ by about 20% from the measured values, except for solar maximum years of SC 23, when the discrepancies in the SW density increase to 40%.

Discussion
The time vs. latitude maps of the SW proton speed and density at 1 AU (Figures 7 and 11) illustrate the variation of the SW structure over almost three SCs. Our analysis begins in 1985, during the low activity level of SC 22 and ends in 2013 in the middle of the maximum of SC 24. During the activity maximum the SW is almost homogeneous in helio-latitude, with slow and dense streams distributed from the equator to both poles. Inspection of the density map ( Figure 11) shows a clear difference in the SW density from the maximum of SC 23. The SW is less dense during the maximum close to the solar equator than it was at the end of 1980s; the density is reduced from about 10 cm −3 around 1991 to ≈5 cm −3 in 2013. This confirms the secular changes in the SW observed in the in-ecliptic measurements (see Figures 1-3 in Sokó l et al. (2013) and also McComas et al. (2013)). The maps show also very distinct solar minima (≈1985, ≈1996, and ≈2009) with a slow (about 500 km s −1 ) and dense (about 8 cm −3 ) SW at low latitudes around the solar equator, and a fast (about 800 km s −1 ) and dilute (about 3 cm −3 ) flow at higher latitudes (i.e., outside the ±20 • − 30 • equatorial band). But the structure of the SW during these two solar minima is quite different; the slow SW extended much higher in helio-latitude (to about ±30 • ) during the last SC than during the previous ones, when it reached to about ±20 • . These differences were observed also by Ulysses (see Figure 6 in Sokó l et al. (2013) and also McComas et al. (2006McComas et al. ( , 2008) and Figure 13 in this article. The slow SW prevails in mid-latitudes much longer after the maximum of SC 23 than after the maximum of SC 22, as can be seen from the latitudinal structure of the SW speed and density for the years 2002-2005 in comparison with those for the years 1991-1993 (Figures 7, 11, and 13).
Our method also correctly reproduced the asymmetry between the north and south hemispheres, manifested mainly during the maximum of SC in the form of phase shift between the hemispheres. This aspect cannot be retrieved from the point measurements from in-situ measurements. The phase shift, noticed already by Sokó l et al. (2013), has been observed, e.g., in the polar coronal holes by Hess Webber et al. (2014) and in the analysis by Tokumaru, Fujiki, and Iju (2015).

Summary and Conclusions
The SW is one of the main factors that create the heliospheric environment from distances within a few solar radii of the Sun to beyond the edge of the solar system, where it interacts with the surrounding interstellar medium. To correctly account for its role in the heliospheric models, one needs to include time variations of its latitudinal structure. The distribution of the SW speed has a significant imprint on the production of energetic neutral atoms (ENA) of various energies (McComas et al., 2012(McComas et al., , 2014. The Interstellar Boundary Explorer (IBEX; McComas et al., 2009b,a) probed the structure of the boundary layers of the heliosphere and its interaction with the surrounding interstellar medium through ENA observations (e.g. Zirnstein et al., 2013;Heerikhuisen et al., 2014;Czechowski, Grygorczuk, and McComas, 2015). Additionally, the 3D structures of the SW speed and density are also needed for the modeling of He ENA (Grzedzielski, Swaczyna, and Bzowski, 2013;Grzedzielski et al., 2014;Swaczyna, Grzedzielski, and Bzowski, 2014), to understand the measurements from the Voyager spacecraft carried out in the inner heliosheath (e.g. Stone et al., 2005Stone et al., , 2008Richardson et al., 2008;and recent works by Provornikova et al., 2014;Grygorczuk, Czechowski, and Grzedzielski, 2015), and the origin and structure of the IBEX ribbon (e.g. Schwadron and McComas, 2013;Funsten et al., 2015). The continuous information about the latitudinal structure of the SW parameters is also essential for modeling of interstellar neutral (ISN) hydrogen (e.g. Katushkina et al., 2013Katushkina et al., , 2014Katushkina, Izmodenov, and Alexashov, 2015) and explanation of the difference between the inflow direction derived from the observations of ISN He and H and its imprint on the asymmetry of the front structure of the heliosphere. Furthermore, knowledge of the variation of the global structure of the SW provides us with information about the long-term processes on the Sun, evolution of its cycle of activity, and differences between the north and south hemispheres.
In this paper we focused on variations of the helio-latitudinal profile of the SW proton speed and density from 1985 to 2013, i.e., during almost three SCs. We are interested in the long-term and large-scale variations of the SW for the purpose of modeling of the shape of the boundaries of the heliosphere and the ionization conditions in the heliosphere from the SW termination shock to 1 AU. For these studies we do not concentrate on the short time and small spatial scale details of SW variations that are necessary in the modeling of space climate or the MHD modeling of boundaries of the heliosphere (e.g. van der Holst et al., 2010;Washimi et al., 2011;Fujiki et al., 2014;Kim et al., 2014).
Our aim was to increase the temporal and spatial resolution of the model developed by Sokó l et al. (2013). We used the SW speed data retrieved from IPS observations from STEL (Tokumaru, Kojima, and Fujiki, 2012) and we focused on reconstruction of the data missing in time and helio-latitude in the CR maps. We adopted the methods commonly used in studies of geophysical data by decomposing them into spherical harmonics to reconstruct the missing regions in longitude and latitude, and the singular spectrum analysis to fill in the breaks in time. We left out the longitudinal variations of the SW speed and restricted the results only to the helio-latitudinal profiles. With this twostep reconstruction algorithm we were able to obtain a continuous time series of the SW speed from 1985 to 2013 with more detailed information about the variations at a given latitude over the whole CR. The details of the reconstruction are described in Appendices A and B. The final time versus helio-latitude map of the SW speed is presented in Figure 7.
After the SW speed was reconstructed, we calculated the helio-latitudinal variations in the SW density using the method that is based on latitudinal invariants of the SW: the SW dynamic pressure (Equation (3)) and the SW energy flux (Equation (4)). The invariants were calculated from the in-ecliptic data taken from the OMNI database, and densities at various helio-latitudes were computed using the SW speeds that we obtained in the present analysis (see Figures 6 and 10). The resulting time versus helio-latitude map of the SW density at 1 AU, complementary to the analogous map of the SW speed, is presented in Figure 11.
The SW speed and density profiles presented in Figures 7 and 11 are available via the online supplementary material.
We cross-validated our resulting SW speed and density data with other available observations, i.e., with the OMNI database for the ecliptic plane ( Figure 12) and with Ulysses measurements out of the ecliptic plane ( Figure 13). The agreement in both cases is very good, except a small deviation in the ecliptic after 2008, which is probably related to the assumed ∆N e ∝ V γ relation used in the CAT analysis of IPS data (see Equation (1)). The agreement with Ulysses measurements is also very good. Additionally, the retrieval of the SW density based on the new method is better than in the model presented by Sokó l et al. (2013). From the differences between our results and the measured in-situ values, we conclude that the average goodness of the reconstruction is ≈ 20%.
set A includes all points on the sphere belonging to one CR map for which we have available SW speed values retrieved by the CAT from the STEL IPS observations (with the spatial resolution reduced to 3 • × 3 • in θ and φ), set B is complementary to set A and contains the coordinates of data gaps: where N A and N B are the numbers of elements in sets A and B, respectively. By f (θ i , φ i ) we denote the SW speed values for the elements of set A, and bŷ f (θ j , φ j ) the estimates of the SW speed for the elements of set B.
The coefficients W ℓm of the decomposition are given by Here Y ℓm (θ, φ) is a real representation of the spherical harmonics. We have The spherical harmonics Y (ℓ, m, θ, φ) were calculated using the built-in function SphericalHarmonicY in Mathematica 8 using the default precision and options.
The assumed maximal order of ℓ gives the maximum spatial size s of the structures that we are able to reconstruct by decomposition into spherical harmonics, according to the relation This means that with a limiting value ℓ max = 12 we are able to reconstruct the spatial structures larger than 15 • . Additionally, the condition ℓ 2 max < N A needs to be fulfilled.
In the first step of analysis we select the maps with the minimal fraction of gaps in the CR maps for the given year, and we replace the elements of set B with the average values over longitude for the given 3 • -latitude according to the description in Section 3.1. The initial elements of set B are replaced with values calculated using the coefficients of the decomposition into spherical harmonics, In the first filling of elements in set B, the gaps in a given map for a selected year are approximately filled. In the next iteration step we search for the W ℓm coefficients with the elements of set B filled with the values calculated from Equation (8). We iterate this process, calculating new sets of coefficients of the decomposition in spherical harmonics, until the average of for all elements in set B is smaller than 0.01, or if the (δ k+1 − δ k ) start to be greater than zero.

B. Algorithm of Singular Spectrum Analysis for Filling the Gaps in the Time Series of Coefficients of Spherical Harmonics
In this part of the analysis we introduce set A of the elements for which we have the coefficients of spherical harmonics, and set B of the elements for which we want to find the values (i.e., the gaps); N is the total number of elements in both sets A and B. The window width M of the maximal periodicity captured must fulfill the requirement of N /M ≈ a few (for discussion see Ghil et al., 2002). We take M = 60, which means that the longest periodicity is determined at most by 60 consecutive points in the time series. The maximal order of iterations J is set to M/2, and the condition for the convergence of the iteration adopted is ǫ = 0.01.
In the SSA analysis we use only the coefficients of spherical harmonics for those CR maps of the SW speed for which the fraction of gaps is smaller than 50%. Additionally, we focus only on the helio-latitudinal profiles, so we limit the A set to the coefficients with m = 0.
Before we start the SSA we correlate the time series W ℓm (t) with a proxy η(t) of the solar activity cycle for elements from set A (following the rationale given in Section 3.2), with t = 1, . . . , N A , where N A is the length of set A. We use the CR averages of solar-radio 10.7 cm flux (F10.7) smoothed by a 13 CR running average (Figure 4). We look for α ℓm and β ℓm from the minimization of the expression t∈A |W ℓm (t) − α ℓm η (t) − β ℓm | 2 . (9) Next we normalize κ ℓm (t) dividing it by the standard deviation (stdDev), X ℓm (t) = κ ℓm (t) stdDev (κ ℓm ) .
Then we apply SSA to the time series X ℓm (t). The SSA analysis of the time series X(t) of the normalized coefficients W ℓm of the elements t ∈ A starts with the initial conditions of the iterators p = 1, r = 1, a normalization parameter S (p, r) = 1, and ∀ t∈A X (p=1,r=1) (t) = X(t), ∀ t∈B X (p=1,r=1) (t) = 0. In the first step, we calculate the matrix C (p,r) with i = 1, . . . , M , j = 1, . . . , M , C (p,r) ij = 1 N − |i − j| N −|i−j| t=1 X (p,r) (t)X (p,r) (t + |i − j|), where X (p,r) is a set with a sum of elements from sets A and B. Next, the eigenvectors and eigenvalues of the matrix C (p,r) ij are found E (p,r) k , λ (p,r) k , respectively, decreasingly ranked with respect to the eigenvalues. In the following step, we calculate and (M t , L t , U t ) are defined after Equation (12) of Ghil et al. (2002) as follows 8 : with N ′ = N − M + 1.
Having computed R (p,r) (t), we calculate the mean square difference Q (p,r) , and assign the new values to the elements of set B using ∀ t∈B X (p) (t) = R (p,r) (t). We normalize all elements from sets A and B in the following way: X (p,r+1) (t) = X (p,r) (t) stdDev X (p,r) (t) .
With the stdDev known for each set, we assign a new value to the normalization parameter S (p, r + 1) = S (p, r) · stdDev X (p,r) (t) . All steps described above are repeated iteratively in two loops. The internal one goes over iterator r with iterator p kept constant, and the external one over iterator p. The internal loop over r checks the value of r, and if r = 1, then increases r by one to r = 2 and starts the calculations from Equations (12) to (16). If r > 1 then the condition Q (p,r) /Q (p,r−1) < 1 − ǫ is checked, and if it is true, then r is increased by one (r → r + 1), but if not, then p is increased by one (p → p + 1), r is again set to be r = 1, and ∀ t X (p+1,r=1) (t) = X (p,r) (t). The external loop over p continues as long as p ≤ J.
The reconstructed (filled) time series of the coefficients of the spherical harmonics are built withX (t) = X (t) t ∈ A, S (p,r) · X (p,r) (t) t ∈ B.
At the end we have to reverse the normalization and go back to the real magnitudes of W ℓm . To do this we use Equations (10) and (11): for all elements fromX (t).