A new 3D solar wind speed and density model based on interplanetary scintillation

The solar wind (SW) is an outflow of the solar coronal plasma, expanding supersonically throughout the heliosphere. SW particles interact by charge exchange with interstellar neutral atoms and on one hand, they modify the distribution of this gas in interplanetary space, and on the other hand they are seed population for heliospheric pickup ions and energetic neutral atoms (ENAs). The heliolatitudinal profiles of the SW speed and density evolve during the cycle of solar activity. A model of evolution of the SW speed and density is needed to interpret observations of ENAs, pickup ions, the heliospheric backscatter glow, etc. We derive the Warsaw Heliospheric Ionization Model 3DSW (WawHelIon 3DSW) based on interplanetary scintillation (IPS) tomography maps of the SW speed. We take the IPS tomography data from 1985 until 2020, compiled by \citet{tokumaru_etal:21a}. We derive a novel statistical method of filtering these data against outliers, we present a flexible analytic formula for the latitudinal profiles of the SW speed based on Legendre polynomials of varying order with additional restraining conditions at the poles, fit this formula to the yearly filtered data, and calculate the yearly SW density profiles using the latitudinally invariant SW energy flux, observed in the ecliptic plane. Despite application of refined IPS data set, a more sophisticated data filtering method, and a more flexible analytic model, the present results mostly agree with those obtained previously, which demonstrates the robustness of IPS studies of the SW structure.


INTRODUCTION
The solar wind (SW), a supersonic stream of plasma expanding continuously from the Sun, creates a bubble-like region in interstellar matter filled with solar wind plasma, called the heliosphere. Solar wind ions and the neutral component of interstellar matter that penetrates inside the heliopause interact by charge exchange. This process operates both inside the solar wind termination shock and in the inner heliosheath, i.e., between the termination shock and the heliopause. The charge exchange reaction products are energetic neutral atoms, which are former SW ions that were neutralized due to transfer of an electron from a hydrogen atom and maintain their velocity from the moment of charge exchange. Some of these former solar wind ions run towards the Sun, where they are detected by neutral atom detectors installed onboard Interstellar Boundary Explorer (IBEX; McComas et al. 2009), thus bringing information about the physical state of the parent plasma in remote regions. However, to correctly interpret this information, one needs to know the input solar wind structure and its variation with time. Similarly, the knowledge of the solar wind structure is needed to account for ionization losses in interstellar neutral gas and energetic neutral atoms inside the heliosphere (e.g., Bzowski 2008. Knowledge of the solar wind structure is also crucial for many other areas of the solar system and space weather research.
SW properties can be determined from measurements in situ. Direct measurements of SW parameters have been performed in the ecliptic plane since the beginning of the space age and nowadays they are cross-calibrated and routinely published in the OMNI2 collection (King & Papitashvili 2005). However, only the Ulysses mission (Wenzel et al. 1989, Bame et al. 1992, Gloeckler et al. 1992) measured in situ the SW parameters in polar regions of the Sun. Ulysses operated in the years 1990-2009, so currently direct and continuous SW speed measurements at higher heliolatitudes are not available. While previously, the out-of-ecliptic structure of SW could be obtained in a few locations from images of cometary ion tails (Brandt et al. 1975), currently the SW speed measurements at high heliolatitudes are performed using interplanetary scintillations (IPS) of radio waves observed on the Earth (e.g., Jackson et al. 1997, Jackson et al. 1998, Jackson et al. 2011, Kojima 1979, Kojima et al. 2004, 2007, Manoharan 1993, Mejia-Ambriz et al. 2015, Tokumaru et al. 2000, Tokumaru 2013, Tokumaru et al. 2015, Tokumaru et al. 2021) and from spaceborne observations of the Lyman-α helioglow (Lallement et al. 1985, Summanen et al. 1993, Bzowski et al. 2003, Lallement et al. 2010, Katushkina et al. 2013, Katushkina et al. 2019. IPS is a random, sub-second fluctuation in the wave intensity of compact radio sources observed by radiotelescopes. They are caused by radio wave diffraction from small-scale (∼ 100 km) electron-density irregularities in the SW. As the irregularities drift across the line of sight, rapidly changing diffraction patterns of radio waves are observed on the Earth as scintillations (Hewish et al. 1964). IPS observations have been performed since the 1960's (Tokumaru 2013). Assuming a frozen-in hypothesis, i.e., that the density irregularities in SW are drifting without change of their properties, and that the SW flows outward radially from the Sun, a cross-correlation analysis of multi-station IPS data allows for the reconstruction of the SW speed along the line of sight. Such speeds are affected by the so called line-of-sight integration effect, which blurs the fine structures. Using tomographic analysis , Jackson et al. 1997, Jackson et al. 1998, the IPS speeds along the line of sight can be deconvolved to a spatial speed distribution, with the line-of-sight integration effect removed, providing 3D Carrington synoptic speed maps of the SW. An example Carrington speed map is shown in Figure 1.
These maps have inevitable gaps (for a broader discussion of the gaps and methods of filling them, see, e.g., Sokół et al. 2015b). However, for accurate studies of global heliospheric structure that can be extended to the outer heliospheric boundary, coverage that is homogeneous in time and heliolatitude is needed. To solve this problem in earlier analyses, the 3D Carrington speed maps are used to obtain average yearly latitudinal speed profiles on a fixed regular grid. This line of research is presented in Sokół et al. (2013), Sokół et al. (2019aSokół et al. ( , 2020. In these papers, the Carrington maps are longitudinally concatenated and averaged, and the resulting mean yearly speed profiles are approximated by an analytic function, which is subsequently mapped on a fixed heliolatitudinal grid for each year. The approximation function introduced by Sokół et al. (2013) is a piecewise first-and second-order polynomial, where approximations of each yearly profile of the SW speed were adopted as piecewise functions composed of straight-lines in polar regions, flatly con- nected to smoothly-connected parabolas between the polar regions. The heliolatitude boundaries between the pieces (the straight lines and the parabolas) are adjusted to maximize the agreement with the data. To calculate the SW speed at a given heliolatitude for time moments between halves of individual years, a linear interpolation between the values obtained for this heliolatitude from the two adjacent yearly profiles is used.
Recently, the entire set of IPS observations was reprocessed using a revised version of tomographic analysis (Tokumaru et al. 2021). The reason for this revision was that discrepancy between IPS and in situ observations became prominent in Cycle 24: the SW speeds derived from the tomographic analysis of IPS observations were systematically higher (by about 71 ± 26 km/s) than those from the OMNI measurements for the period after 2009. This discrepancy was found to be improved by changing the value of the power law index α in the empirical relation ∆n e ∼ V α that is assumed in the tomographic analysis. Here, ∆n e and V α are the SW density fluctuations and speed, respectively. In the revised version of tomographic analysis, the index α was assumed to be variable, depending on the year. The variation in α is considered to be closely linked with the long-term decline in the solar activity. By allowing this change in the IPS analysis, the discrepancy between the IPS and in situ observations was reduced to 35 ± 23 km s −1 for the period after 2009, even though the tendency to overestimate the SW speed in the reprocessed data set still remains.
Here, we overhaul the approach originally defined by Sokół et al. (2013). We take the most recent version of the IPS SW speed maps (Tokumaru et al. 2021), we introduce a novel method of filtering these maps against outliers, and we use a different approximating function, based on Legendre polynomials. The model of yearly-averaged SW speed and density is extended until 2021.
The methodology of our SW speed modeling consists of two main parts: initial processing of the IPSderived Carrington speed maps and searching for the best yearly fits to longitudinally averaged yearly speed maps. The main goal of the initial processing is to provide filtered, averaged and binned yearly latitudinal speed profiles of SW without bias. Bias removing is done by filtering out the map background before averaging the speed maps. The filtering method is based on statistical analysis to avoid any arbitrary cuts.
After filtering, the Carrington speed maps are averaged in the selected bins over CRs to provide the Carrington-mean latitudinal SW speed profiles. The binning is done for each of the CRs with IPS data. The binned data are then used for fitting an approximating function to the yearly SW speed data. The approximation model we use is more flexible than that used by Sokół et al. (2020). The speed profiles are defined as sums of Legendre polynomials with appropriate coefficients, with an additional requirement that the derivative of the profile over the pole is 0. The range of the model used varies between the years, depending on the data. The fitting is done simultaneously to all Carrington-averaged profiles from a given year. Each year is fitted separately. Here again, to eliminate arbitrary decisions, during the fitting procedure we apply a rigorous statistical analysis to identify the optimum range of the model for the given yearly SW speed profile model.
In the final step, an analysis of differences between the model and OMNI in situ measurements is performed and the yearly model speed profiles adjusted, similarly as it was done by Sokół et al. (2020). The SW density for a given year is calculated based on the adjusted yearly speed profile and the magnitude of the latitudinally invariant SW energy flux for the given year (Le Chat et al. 2012).
The input data are presented in Section 2.1. The data filtering method is presented in Section 2.2. There, we estimate the mean bias as a function of latitude and describe its properties. Subsequently, we present the results of filtering the data in Section 2.3. The approximation model is presented in Section 3.1. The final SW speed binning and model fitting is discussed in Section 3.2. Then, the SW speed model is adjusted to the OMNI2 speeds in the ecliptic plane (Section 3.3). The results are discussed in Section 4, and the density model is calculated using the SW energy invariant measured in situ, as presented in Section 5. The results are discussed in Section 6, and we finish with summary and conclusions (Section 7).
2. DATA 2.1. General description of the data As input to the analysis, we use the most recent IPS Carrington speed maps publicly available from the Institute for Space-Earth Environmental Research (ISEE) at Nagoya University, Japan (Tokumaru et al. 2010, Tokumaru et al. 2021. These maps were derived using the Computer Assisted Tomography (CAT) method applied to the ISEE IPS observations (Tokumaru 2013), covering the entire interval of available observations from 1985 through 2020, with a gap in 2010 because of antenna system modification. These maps are available on a regular 1 • by 1 • mesh in heliographic coordinates. The maps are provided in time intervals corresponding to individual Carrington rotations of the Sun.
Because of snow loading on the two movable antennas in winter, the IPS velocity data collection is stopped from December to February. In effect, regular gaps in the IPS observations are present in each yearly Carrington speed maps for the CRs during this period, and the available number of CRs in the speed map for a given year is about 11; the boundaries of Carrington periods do not coincide with the beginning of the year, and thus the gap in data coverage over the winter months is variable.
In addition to these regular gaps, some additional gaps may appear due to observing conditions and the sparsity of the available radio sources in the sky. The ISEE IPS data are collected using a UHF radiotelescope array, which consists of antennas operating at a frequency 327 MHz. The data are collected from observations of 15-20 compact radio-sources (mostly quasars, observed daily at small angular elongations from the Sun). In winter months, the Sun is low on the horizon, and very few sources are able to be viewed south of the Sun at this time. Thus, sources in the regions southward of the Sun are not always sufficient to provide a uniform coverage of the Carrington speed map in the southern hemisphere. This factor also affects the mean speed profiles, because a lower number of data points may cause a larger scatter of binned speed points near the south pole of the Sun.
The irregular distribution of radio-sources in the sky is one of the reasons why the accuracy of SW speed reproduction varies within the Carrington maps. In particular, when very few radio sources are available in a large region in the sky, the tomography analysis sometimes returns outlying very low or very high speeds.
In the original analysis by Tokumaru et al. (2021), the boundaries of allowable speed values were set to 200 km s −1 at the lower end and 850 km s −1 at the high end of the speed range.
Based on comparison of IPS-derived and Ulysses in situ observations, Sokół et al. (2013) noticed that generally, individual monthly Carrington maps with less than 30 000 data points have a lower accuracy than the maps with a larger number of points, and to prevent a bias, they rejected these sparsely-populated maps from further analysis. But these sparsely-populated maps may include well-populated regions that can safely be used in the analysis. Therefore, we devised a data filtering method based on statistical search for outliers and employed it as the first step in our data processing. This method and its results are presented in the next section. To avoid bias arising due to the map background, we treat the Carrington maps as raw data, and we introduce additional processing of the maps before averaging. The need for data filtering is clearly seen in a density histogram of SW speeds, taken from a yearly Carrington map concatenated over heliolongitude, presented as an example in Figure 2. In the density histogram of speeds the majority of speed points generally follow a latitudinal V-shaped profile, which is a characteristic solar wind feature during solar minima, like in the example year 1996. But some of the speed points fall far outside the main V-shape area. Some of them reach the low-speed limit at 200 km s −1 , even though the majority of speed values at high heliolatitudes are about 800 km s −1 . The low speed values at higher heliolatitudes are not likely to be realistic during solar minimum, and will be referred to as the map background. It must be stressed here that some of the spread in the yearly-accumulated data is due to physical variation in SW speeds.

The need for data filtering
To identify and reject outlying data points, we used a novel statistical method based on the Generalized Extreme Studentized Deviate (ESD) Many-Outliers detection procedure by Rosner (1983). In this method, it is assumed that the data sample is a mixture of a certain number of outlier points in addition to the "correct data" population conforming with the normal statistical distribution. This method hypothesizes that the sample includes up to n outliers and assumes that they deviate greatly from others in the sample. We start by selecting n points the most distant from the median of the sample. We calculate the test statistics and critical values for each of the suspect-outlier points from 1 to n and reject the point for which the test statistic value is higher than the critical value, obtained from the Student T distribution. We reiterate this procedure: after each step, the i-th observation is rejected, and the statistics and the critical value are recalculated. The procedure returns the number of outliers, which is the largest i ∈ n at which the i-th test statistics is larger than i-th critical value.

Filtering Carrington speed maps against background
The analysis begins with application of initial map filtering procedure, which identifies and removes speed points belonging to the map background. The procedure is performed separately within individual heliolatitude bins, separately for each year, but the speed points from the whole year are analyzed simultaneously. Also, each heliolatitude bin is analyzed independently.
The filtering procedure is based on the fact that the map background points feature a large deviation from the median value of the distribution in the respective heliolatitude bin. We use the median value because the speed distributions are strongly asymmetric, as shown in the Figure 3 example. Not all speed distributions in heliolatitude bins are asymmetric but the use of median is not a problem in these cases. We assume that the background points are outliers of the speed distribution around the median value, so the map filtering relies on identification and removing of the outliers in each bin.
During the filtering, we operate on 1 • heliolatitude bins, coincident with the maximum resolution of the original Carrington maps. In each 1 • heliolatitude bin, we use all speed values in a given calendar year, without dividing them into individual CRs. Although we use narrow bins, the number of points in each 1 • heliolatitude bin is large. This is because the speed points in the Carrington maps do not represent the measured speed values, but they are a result of CAT reconstruction on a dense grid of heliolongitude and heliolatitude. Therefore, we treat the points in each heliolatitude bin as oversampled data, which in general follow the Student T-like distribution.
To identify and remove the outliers, we transform the speed distributions by calculating absolute values of speed differences between individual speed points and the median speed within a given bin. Next, we search for the maximum number of outliers in the distribution of absolute speed differences, using the one sided extreme Studentized distribution (ESD) test (Rosner 1983). After identification of the outliers, they are removed from the Carrington maps before averaging.
Even though for a given year the speed points from all CRs are analyzed at once, the information about CR location for each point is preserved and used subsequently during averaging. Since we operate on the original-resolution Carrington speed maps, we need to process 180 1 • equatorial bins per year. The central positions of the latitudinal 1°bins correspond to the original resolution of the Carrington speed map. This processing is time consuming because the maps consist of a large amount of speed points.
No initial averaging of any parts of the maps is performed. The sharp edges of the blank areas on the maps precludes this because averaging in the presence of such edges would result in reduction of speed values near the blank areas. Proceeding with the maximum resolution speed maps is done to preserve as many points as possible, since the data have extensive gaps. In this way, we can use speed points at the south pole belonging to CRs with extensive gaps of data. Preserving this information is important because of generally poor map coverage in the southern hemisphere is motivated by our intent to sustain the available statistics as much as possible. For example, for CR1908 in the yearly map presented in Figure 1, the southern area is poorly covered, but the speed points existing in this CR are used, while in previous models such CRs would be rejected entirely.
The filtering method is governed by only one parameter, i.e., the confidence level (CL). We use the same CL for all histograms in the analyzed data set, but the optimum CL value to use was subjected to further analysis. Different values of CL will effect in different intensity of map filtering. The optimal filtering intensity, driven by the CL value, is determined as follows. For each year, we calculate the mean SW speed profiles for the map with background (i.e., for not filtered). Subsequently, we repeat the calculations for the map with the background removed, using several typical values of the CL as a parameter of the one-tailed statistical test: 95%, 97.5%, 99.0%, 99.5%, 99.9%, 99.95%. Then we compare how the correction profile averaged over all years changes after filtering depending on the CL used. The magnitudes of such changes as a function of latitude averaged over all years can be regarded as a typical correction of the latitudinal profiles, obtained with respect to the profiles modeled without map filtering (see Figure 4). The behavior of the this correction for different levels of the filtering intensity allows us to select the optimum CL.  This is done by inspection of the mean correction profiles for different CLs. We found that the correction profiles have smooth shapes, and are similar for different CLs. An evolution of the correction profile with CL value is shown in Figure 4. Clearly, the modifications of a yearly profile due to filtering for various CLs has similar features. The interpretation of the typical correction is that on average, the map filtering increases the speed values by ∼ 25 km s −1 at both poles. This means that typically, the bias due to the map background causes an underestimation of speeds at the poles. At the equator, the filtering results in decreasing the modeled speed by ∼ 10 km s −1 . This is important also in the context of differences between the IPS and in-situ results in the ecliptic plane, discussed in Section 3.3. The typical correction is northsouth asymmetric. The asymmetry may arise because of non-uniform statistics and data coverage at the polar regions.
The filtering results in a systematic modification of speeds at all heliolatitudes. At the equator, probably the adopted lower speed limit of 200 km s −1 slightly biases the mean profile gathered from the raw maps. The latitudes at which the speed limits do not introduce any difference in the mean sense are around ±25 • , but this may vary between the individual years. It is also seen that for the CL 99.5% and less, a stabilization of the mean correction value occurs. This behavior suggests that, starting from a certain value of CL, the filtering method becomes stable for smaller CL values, otherwise the mean bias correction profile would increase with the CL decrease. This also confirms that the mean profile distortion, caused by the map background and the 200 km s −1 speed limit, is compensated for, since further elimination of outliers with a smaller CL brings negligible effect. However, the smaller CL, the lower deviation of points around the model fit, which leads to an underestimation of the accuracy of the model. Therefore, to avoid this effect, in further analysis we use CL = 99%, i.e., the largest CL at which the mean bias correction profile becomes stable.
With this, after identification and rejection of the map background, we obtain filtered Carrington speed maps, which are used for further steps in the model construction. The average yearly latitudinal profiles from these maps have a reduced or eliminated bias. The effect of filtration can be appreciated in Figure 5, where the histogram of the data from the raw map from 1996 is compared with the histogram of filtered speeds, and in Figure 6, where heliolatitude histograms of speeds before and after filtering are presented. The value of correction for individual years increases with heliolatitude, with extreme values from −35 km s −1 in the equatorial area to 78 km s −1 at higher latitudes.   Figure 6. Example density histograms of speed points taken from the 1996 Carrington speed map before (left) and after (right) filtering, plotted as a function of heliolatitude. Extensive elimination of the map background due to filtration is visible.

Preparation yearly profiles for fitting
Next, to prepare the yearly data for speed profile fitting, we change the domain of the Carrington speed maps from heliolatitude (φ) to cos(π/2 − φ), because we want to model the speed profiles using orthogonal functions at the −1, 1 interval. A division of the domain into bins with the same width yields equi-areal latitudinal bins, which means that the bin widths in degrees towards the solar poles are wider than those at the equator. The widening of the bins is a desirable feature, because the number of available points towards the poles decreases due to the lack of IPS data. With the equi-areal bins, the statistics in the bins is more homogeneous, and will result in reduction of fluctuations at the poles in the fitted models. An empirical study showed that the optimal latitudinal bin width is 0.05, i.e., we use 40 bins in the entire latitudinal −1, 1 range. We found that a change of the bin width by a factor of 1.5 does not substantially affect the results.
Subsequently, we average the speed points within the bins, separately for each CR. The speeds in each bin are averaged selectively, after being gathered from the CR to which they belong. So, for a given year, we obtain as many binned profiles, as there are CRs available that year. This is different compared with the processing performed during the filtering. During the filtering, all speed points from the entire year were processed together, but now data from individual CR are treated separately by using identical bins. This binning provides profiles for all CRs for which the data are available in the original data set. The data prepared in this way from all CRs for a given calendar year are used in the fitting of yearly model profiles. Some of the profiles, however, have gaps in the latitudinal coverage, and can be dealt with as demonstrated below.
An example visualization of binned points, ready to fit, is shown in Figure 7. A comparison of binned profiles for the raw and the filtered Carrington maps are presented in Figure 8. The data in Figure 8 are the same as in Figure 7; the binned profiles are shown with their CR numbers, to visualize the impact of map filtering. It is seen that after filtering the mean CR profiles are less noisy, but nevertheless a certain noise still persists. Clearly, the filtering removes many of the profile distortions and artifacts. The entire range of cos(π/2 − φ) is presented, as it is used for modeling of the mean yearly profile. The red points are used for fitting, the gray points have been filtered out. An example of the distortion removal can be appreciated in the profile from CR1917. This profile was obtained at the end of the IPS data collection period in 1996. The lack of sufficiently abundant IPS data in this CR caused that the profile is shifted down in the northern hemisphere, and is very distorted in the southern hemisphere. The filtering procedure removes the distortion, and moves the profile into an adequate place from the statistical point of view for the assumed CL. That in the case of extreme distortions the profiles are just cut out means that no extraction of information about the actual shape of the profile was possible because of insufficient data amounts.
To prepare data for fitting yearly model speed profiles, we do not average the profiles from individual CRs over a calendar year. This is to preserve the SW speed variations between the CRs during a given year. We do so as a result of an empirical study that we performed, which showed that the point distribution between individual CRs make the fitting procedure more efficient and reliable. The yearly models are fit to binned heliolatitude profiles for individual Carrington rotations for a given calendar year.
The final step of profile preparation for fitting is averaging the speed within the 40 bins.
3. MODELLING OF THE SOLAR WIND SPEED

Approximation model
The modeling is done using averaged speed profiles like that showed as an example in Figure 7. We model the SW speed profiles as a sum Legendre polynomials. These polynomials have the necessary flexibility needed to reproduce the various shapes of the profiles we need to approximate. However, tests showed that an additional condition of nulling the derivatives of the model at the poles is necessary to provide a general stability of the model and to avoid the propagation of the effects of low coverage of the polar regions into lower latitudes. Without these conditions, an unrealistic behavior of the fits at the poles sometimes occurs. Also, the use of binned speed values contributes to instability of the unconditioned model at the poles because of the lack of binned values exactly at the poles, i.e., at cos(π/2 − φ) = ±1. For this reason, we set the first derivatives of the Legendre polynomials at the poles to zero.
Denoting z = cos(π/2 − φ), a general form of the model is the following: with the conditions at the poles (i.e., at the points where z = ±1): where P i is Legendre polynomial of i-th order, Q i is the coefficient of the i-th polynomial, and N is the approximation order. Applying these conditions eliminates two of the free parameters Q i . The formulae for V(z) used in the fitting are algebraically rearranged by deriving a linear combination of the Legendre polynomials for a given N, in which the two coefficients standing at the highest orders of the Legendre polynomials are eliminated. The remaining free coefficients are obtained from the fit to the data. Derivation of the formulae was done using algebraic manipulation and differentiation in Wolfram Research Mathematica. Legendre polynomials were obtained from Mathematica function LegendreP. The analytic formulae are too large to show them here, especially since they change from one order to another with the change of the order used because of the additional condition specified in Equation 2.

Fitting yearly model speed profiles
The least-squares fitting was performed using Mathematica function LinearModelFit. Wolfram Research Mathematica version 12 was used. The mean SW speed changes from year to year. Depending on the year and the solar cycle phase, the complexity of the SW speed profile shapes varies. Since the Legendre polynomials are very agile, using too high an order in the fitting procedure may result in reproducing the remaining small data fluctuations in the model. The goal here is to avoid overfitting, i.e., in reproducing some remnant random variations in the data. In the case of the former model by Sokół et al. (2020), this issue was absent because of the simplicity of the approximating function (low-order polynomials). However, the same factor may have prevented that model from reproducing possible legitimate variations of solar wind speed with heliolatitude, like those proposed by Katushkina et al. (2013), Katushkina et al. (2019), who reported that the solar wind flux maximizes at mid-latitudes in the northern and southern hemispheres based on analysis of the helioglow distribution in the sky.
Our analysis showed that models of different orders are optimal in different years. The choice of the optimal polynomial order is done for each year separately, based on statistical analysis of the fit residuals. We set the lower order of the Legendre polynomials to 8, and the upper one to 20, and we assumed that the optimal order of the Legendre polynomial should return residuals with the distribution characterized by the highest p-value of the Pearson χ 2 test of normality.
However, using the ordinary χ 2 measure was abandoned because of the lack of the uncertainties for individual speed points in the input maps. Using the dispersion in individual bins as an estimator of the uncertainty measure, which might be used in the ordinary χ 2 measure, led to wrong results. This can be understood when one realizes that the input data points in high-resolution Carrington maps are results of mapping the results of CAT analysis on the coordinate grid. Therefore, the data points are not independent and may feature correlations within a certain spatial range in the maps. An example of such correlation can be seen in Figure 7, where among the gray points (filtered out from the sample), series of data points forming continuous gradients are visible. This illustrates that the CAT method used to derive the map can introduce correlations between nearby data points.
An example best fit for 1996 is shown in Figure 9. Having fit the SW speed model for all years, we performed a reverse transformation of latitude from cos(π/2 − φ) to φ, and we apply the OMNI2 adjustment (see Section 3.3).

The OMNI2 adjustment
A good agreement between IPS CAT speeds and Ulysses measurements was reported by Tokumaru (2013). However, in certain time intervals, systematic differences between IPS-derived speeds and those measured in situ and compiled within the OMNI2 collection (King & Papitashvili 2005) have been noticed (Sokół et al. 2020, Tokumaru et al. 2021.  Figure 10. A comparison of the mean yearly speeds obtained from filtered Carrington maps at the solar equator with the yearly averaged OMNI2 speeds. The magnitudes of the uncertainty bars are equal to 1σ taken from the standard deviation of the CR-averaged speeds at the equator. Differences between the SW speeds obtained from the ISEE IPS measurements and Ulysses data were also studied in detail by Sokół et al. (2013), where the overall usefulness of the IPS measurements at all latitudes was demonstrated for the years of the Ulysses mission. However, Sokół et al. (2020) reported that the mean yearly IPS speeds at the equator taken directly from the Carrington maps showed a divergence ∆v IPS −OMNI = v IPS − v OMNI after 2011, i.e., after the year of the ISEE IPS antenna system upgrade. For the revised up-to-date IPS speeds, a comparison of the IPS speeds with OMNI2 time series is presented in Figure 10, where the difference, previously reported by Sokół et al. (2020) for the previous version of IPS speed maps, is generally reduced, but is still present. Tokumaru et al. (2021) present a detailed study of the differences between the SW speed derived from IPS tomography and in situ measurements performed by Parker Solar Probe and those available in the OMNI2 collection on a much finer time scale than we have in this paper.
We studied in detail the ∆v IPS −OMNI time series. Although a look at the time evolution of this quantity may suggest it features a simple linear trend (see Figure 11), we performed a statistical analysis of the goodness of fit for several different functions. In addition to a linear trend, parabolic and step functions for the fit to the time series of ∆v IPS −OMNI (t) were studied. The boundaries in the step function were used to quantify a hypothetical coincidence between the antenna setup and changed the level of agreement of the IPS data with the OMNI2 measurements. The better agreement between the IPS and the OMNI2 data that exists during the years of the Ulysses mission, when in situ measurements were used for the IPS calibration of raw data, is another reason to consider the use of a step function to derive the adjustment for the IPS yearly profiles.
The three alternative models have different numbers of free parameters. We performed the fits to the three functions using Wolfram Mathematica and checked if the use of a model with larger parameter numbers yielded a statistically significant improvement of the model. We found that the values of the Akaike Information Criterion (AIC) for the three models are similar, and none of the models seems to be significantly favored from the current statistics. A possible explanation for the behavior of the systematic differences in time was pointed out by Tokumaru et al. (2021) who suggested that it might be due to the effect of time variation of the index α, as mentioned in Section 1. The physical property of SW microturbulence is likely to evolve differently from the past cycles in accompany with the significant weakening of the solar activity for Cycle 24. According to the earlier study (Tokumaru et al. 2018), a pronounced growth of the low-∆n e stream occurred for the slow SW in Cycle 24. This may result in a drastic modification of the ∆n e -V relation used in IPS SW speed reconstructions, particularly for the slow SW. The IPS-OMNI comparison is considered to be seriously affected by this change, because the OMNI data were collected in the ecliptic plane, where the slow SW dominates. The change in the ∆n e -V relation revealed here may suggest that the fractional density fluctuations ∆n e /n e became more dependent on V in this cycle than before, providing an important insight into the solar wind acceleration mechanism. While this approach resulted in a reduction of ∆v IPS −OMNI in comparison with these reported by Sokół et al. (2020), the present version of the IPS Carrington speed maps does not eliminate them altogether.
A detailed look at individual years on the ∆v IPS −OMNI plot shows a large difference appeared for year 1994, which was not present before. Additionally, a small systematic difference is visible for the years 1985-1992. In this period, all values of the IPS speeds are slightly lower than the OMNI2 speeds, which was not present in the previous version of the IPS data. However, taking into account the uncertainty, the differences seem to be reduced when compared with the previous version of the IPS speeds, but the divergence after 2011 is still noticeable. Almost no divergence is observed during the years of the Ulysses mission.
Since the OMNI2 time series is compiled from many different experiments and measurement methods, which are carefully curated (King & Papitashvili 2005, see also on-line description of the data scaling and calibration at https://omniweb.gsfc.nasa.gov/html/ow data.html), we decided to adjust the model derived from the IPS Carrington maps so that it agrees with the OMNI2 collection in situ time series.
After a detailed study of the ∆v IPS −OMNI , taking into account that the considered differences coincide with the above-mentioned variations of ∆v IPS −OMNI in different periods of data set, and having considered several potential methods of the adjustment derivations, we concluded that the optimal strategy for adjusting the IPS data was the method proposed by Sokół et al. (2020). Thus, we calculate the difference ∆v IPS −OMNI separately for each year and subtract it from all latitudinal bins. We found that only 9% of the adjustments are statistically significant, i.e., larger than the joint uncertainties of the yearly averages of the OMNI2 and IPS speeds. For the entire data set, this adjustment provides very small differences for the years that show a good agreement with the OMNI2 time series, while it reduces large values of ∆v IPS −OMNI for the years when these differences are large, like 1994, 2014, and 2020. The magnitudes of the adjustment values for individual years are listed in Table 1.

SW SPEED RESULTS
At this stage, with the OMNI adjustment applied, the model can be used directly to calculate the SW speed at an arbitrary latitude during a given year. The fitted coefficients in the formula defined in Equations 1 and 2 for all years for which IPS measurements are available are listed in Table 2. Note that application of the OMNI2 adjustment results in a decrease of the mean value of the speed over the full range of heliolatitudes.
The parameter of the model that is equal to the mean value is the parameter number 1. Hence, with the adjustment magnitudes available in Table 1 and the profile parameters in Table 2, the user can decide which version of the model to use: adjusted or not. The final step in the construction of the solar wind speed product in the WawHelIon (Warsaw Heliospheric Ionization) model of the factors relevant for modeling ionization losses of neutral species in the heliosphere is a projection of the speed model on the common time and heliolatitude grid. This grid is identical to that used in the models of the heliospheric ionization factors by Sokół et al. (2019aSokół et al. ( , 2020. This facilitates changing between ionization models in codes using the heliospheric ionization factors, such as the Warsaw Test Particle Model (WTPM) of interstellar neutral species distribution in the heliosphere (Tarnopolski & Bzowski 2009, Sokół et al. 2015a, Sokół et al. 2019b and in the WawHelioGlow model of the heliospheric Lyman-α backscatter glow (Kubiak et al. 2021a,b), as well as in the calculation of survival probabilities of heliospheric ENAs between the termination shock of the solar wind and the IBEX detectors (Bzowski 2008, McComas et al. 2012, McComas et al. 2017.
The time-heliolatitude grid in the time domain is based on a Carrington period mesh, with the nodes at the halves of the Carrington periods. In the heliolatitude domain, the grid nodes are located at multiples of 10 • , uniformly distributed between the solar north and the south poles. The SW speed is projected on this grid. The yearly SW profiles in heliolatitude are linearly interpolated to provide speed magnitudes at halves of Carrington rotation periods and on the heliolatitude grid.
We assume that the SW speed thus obtained is independent of the heliocentric distance. This assumption is well fulfilled between a few tenths of au, i.e., outside the region where the SW is being accelerated, and inside ∼ 5 − 10 au. Beyond these distances, a gradual, direction-dependent slowdown of SW speed due to mass-loading from the charge-exchange between the SW and interstellar neutral H begins (Isenberg 1986, Richardson et al. 2008, Nakanotani et al. 2020).
The solar wind model we are constructing is an element of a broader Warsaw Heliospheric Ionization (WawHelIon) model, devised as a source of ionization rates of heliospheric species for interpretation of measurements. Therefore, to accommodate the much more detailed information about the solar wind parameters available at the ecliptic plane, we replace the 0 • heliolatitude speeds with Carrington-averaged OMNI2 speeds, interpolated to halves of Carrington periods. The ±10 • nodes of the grid are filled with mean values between the equator and the latitudes −20 • and +20 • , respectively. This is the last step in our SW speed modeling. The final result of the heliolatitudinal SW speed evolution is presented in Figure 15. The newly-constructed model will be referred to as the WawHelIon 3DSW. The fitted models of the mean yearly SW speed profiles are presented in Figures 12-14. The model of the SW speed, averaged in three equatorial bands identical to those used for presentation by Sokół et al. (2020), is shown in Figure 23. The results from this latter paper are shown along with our present results for comparison, and a discussion is provided in Section 6.3.

DENSITY
The SW is composed of fast and slow wind, and features a strong anticorrelation between the flow speed and the density, as clearly demonstrated by the Ulysses mission (McComas et al. 2000. Both fast and slow solar wind seem to carry the same energy flux. As shown by Le Chat et al. (2012), the SW energy flux is invariable in heliolatitude within 10% and varies weakly with time. Hence, it can be regarded as a quasi-invariant in heliolatitude and used to infer the density of the SW at an arbitrary heliolatitude provided that the speed at this latitude is known and the magnitude of the SW energy flux is known at a different latitude, e.g., in this case in the ecliptic plane. This observation has been used in the modeling of the SW structure for several years (McComas et al. 2014, Sokół et al. 2019a and we use it in this paper. The speed-density relation, which is linked by the SW energy flux, allows us to determine the 3D SW density model. We use the relation given by Le Chat et al. (2012), also used in the previous model by Sokół et al. (2020), which has the following form: with where t is the time, m p the proton mass, m α the alpha particle mass, W the solar wind energy flux in Wm −2 , ν p the SW speed in km s −1 , n p the solar wind density in cm −3 , and the constant C = GM R −1 is the Sun's gravity potential at the photosphere (i.e., at R ). The SW speeds used to estimate the 3D SW density model are taken from the 3D SW speed model described above. The SW energy flux is obtained from in situ measurements of the SW density, speed, and the alpha-to-proton abundance n α /n p , taken from the OMNI2 time series. To find the most useful way to obtain the yearly SW energy flux magnitudes to be used for speed calculations, we first calculated Carrington period averages of the SW energy flux in the ecliptic plane. Based on these values, we computed mean values of the energy flux for individual years using the same Carrington rotations that were used for the given year to obtain the yearly speed profiles. The yearly mean values of the energy flux thus obtained are shown in Figure 16. It turned out that these values are very close to the simple yearly averages for all years in the IPS data interval.  The yearly mean values of the SW energy flux were then inserted into Equation 3 for all years and all heliolatitudes in the grid. This produced the yearly density profiles. Subsequently, we linearly interpolated the yearly density profiles to all CR halves similarly as we did for the speed profiles. Finally, like in the 3D wind speed model, the density values at the equator were replaced with the density values from OMNI2, while the ±10 • bins were taken as average values from 0 • and ±20 • , respectively. The resulting model of evolution of the SW density structure is presented in Figure 17. The density model in its full resolution, organized in five heliolatitudinal bands, is presented in Figure 25, along with the previous model from Sokół et al. (2020) for comparison. A discussion is provided in Section 6.4.
The density is assumed to drop with the square of heliocentric distance. This assumption is well fulfilled within a similar distance range to that for the SW speed. The density of the core SW decreases steeper than this because of the losses due to charge exchange with interstellar hydrogen atoms. However, outside ∼ 5 − 10 au, in addition to the core SW population, a pickup ion population appears. It comprises former neutral H atoms that were ionized and subsequently picked up by the magnetic field frozen in the SW plasma. Pickup ions form a separate proton population in the SW. When treating the core and the pickup ion SW populations together, the total SW proton density change becomes, in fact, a small density gain relative to the 1/r 2 relation. This is because pickup ions are created due to photoionization and ionization by electron impact. While charge exchange effectively shifts a proton from the core SW to the pickup ion population without an overall change in the proton density, the other two reactions effectively inject new protons to the pickup ion population. These "excess" pickup ions are former interstellar H atoms that were ionized by photoionization or electron ionization. This excess is small because the photoionization and electron impact rates are small for H atoms in comparison with the charge exchange rate (cf., e.g., Figure  B1 in Sokół et al. 2020), which shows that the charge exchange rate is responsible for 0.7-0.8 of the total ionization rate). Hence, the assumption of a 1/r 2 drop of the total proton density in the SW is relatively well fulfilled. However, an additional increase in the excess over 1/r 2 relation is due to the slowdown of the SW speed due to mass-loading, discussed in Section 4. Detailed analysis of these subtle effects requires accounting for interaction between SW and interstellar H gas and is outside the scope of this paper. 6. DISCUSSION 6.1. Is Carrington speed map filtering advantageous?
While filtering the data, we used the same CL for each year, to maintain the same intensity of map filtering. This is done to provide the same level of noise reduction for the entire data set, and to avoid scaling of the scatter of the residuals. Such scaling would influence the estimated uncertainty of the model. Having the same filtering level, we are able to compare the accuracy of the model with and without filtering of the map background for all individual years and to compare the effectiveness of filtering between the years. This comparison is presented in Figure 18.  Figure 18. Comparison of the SW speed uncertainty for the models obtained from the raw and the filtered Carrington speed maps.
Despite the use of a homogeneous filtering intensity, the SW speed model features different accuracies in different years. It is probably due to different scatter of the IPS data in different years, combined with different true SW speed variations during individual years. But an improvement of the SW speed model accuracy obtained owing to map filtering is clearly visible. While the uncertainty of the model without filtering is on average ∼ 70 km s −1 , after filtering it is reduced to ∼ 40 km s −1 . An improvement of the accuracy of the model due to filtering is present for each year. For some years, the map background reduces the model inaccuracy substantially, e.g., for the year 1996, for which the uncertainty is reduced by a factor of ∼ 3.5.
Looking at the individual years, the years 1986, 1996, 2005, 2014 and 2020 are characterized by outstanding values of the model uncertainty, when compared to the typical uncertainty values. For these years, the uncertainty is about 100 km s −1 and it is significantly reduced after map filtering. A closer look at the binned speed values for these years shows the reason of such a large uncertainty for the models obtained using the raw maps (see Figure 19). In the Figure, the black points represent speed values obtained from a raw map, the red points those from a filtered map. Some of the binned profiles obtained from the raw map sharply drop to 200 km s −1 for some CRs. Such profiles are usually from the CRs at the edge of the yearly Carrington speed maps, when the coverage of data is not sufficient for a proper reconstruction of the Carrington speed maps. Our filtering procedure removed the biased points. Without filtering out the biased points, a Legendre fit to such data would be distorted, which is shown in Figure 20.
In this example figure, the fit for the raw map is shifted to lower speeds, and additionally a wavy behavior of the fitted model appears. The waving arises when the Legendre fitting procedure ties to model fluctuations of speed points caused by the biased speed values. The difference between the two fits reaches 80 km s −1 at higher heliolatitudes, which results in the biased 3D SW speed model. Also, the speeds averaged with the map background cause the profile to be wider in heliolatitude (notice that the red -filtered -profile is inside the black -unfiltered profile). The filtering method is effective even for very noisy years, providing filtered speed maps, from which SW speed models are characterized by a very homogeneous final accuracy.
Summarizing this part of the discussion we point out that the filtering of the Carrington speed map background is an indispensable element to provide a reliable model based on the Legendre functions, without extensive limitations of the fitting function order. An arbitrary Legendre order limitations would be necessary to prevent introducing quasi periodic fluctuations of the model as a function of heliolatitude, but then the order would be probably too low to reproduce the salient details of the SW speed profile. Modeling of the SW speed profiles using the Legendre polynomials without map background filtering provides biased results, additionally affected by unrealistic fluctuations. Thus, the use of the filtering is advantageous not only because it removes the above mentioned drawbacks, but also because one is able to obtain objectively accurate modeling results with a reasonably low uncertainty.
Another issue is that the homogeneity of the model accuracy does not seem to be entirely random, unlike what one might expect. This is shown in Figure 18. We may distinguish three intervals in the uncertainty plot, for which the uncertainty values have similar properties. The boundaries of the intervals seem to be identical with the periods of the ∆v IPS −OMNI . Namely, the intervals are before 1995, and after 2011. Before 1995, almost for all years the uncertainty value oscillates around 50 km s −1 , to drop in the next interval, which covers the Ulysses operation years. Since 2011, the uncertainty seems to increase with time, but these variations are too low to make a statistically significant statement. Further monitoring of this uncertainty evolution will be necessary.

The significance of OMNI2 adjustments
The OMNI2 adjustments are applied for each year individually. We check the significance of the ∆v IPS −OMNI corrections compared with the effective uncertainty of IPS and OMNI2 in a given year. First, we compare the relative uncertainties of the IPS and OMNI2 speeds at the equator to see how they differ. Such comparison is shown in Figure 21. The uncertainties of the IPS model are estimated from the model residuals for the equatorial area only. Since the IPS data are organized into CRs, the IPS uncertainties at the equator are mainly associated with the intrinsic variability of the average equatorial SW speeds between individual CRs in a given year. In the case of the OMNI2 data, their uncertainty is estimated from the scatter of the measurements. During estimation of the OMNI2 uncertainties, the OMNI2 data are also organized into time basis similar to the manner used for the IPS data, i.e., according to the individual CRs, which is necessary to provide the same time basis for comparison studies and to smooth variability in smaller timescales. The individual CR numbers selected for the calculation of the uncertainty in a given year are the same for both the IPS model and the OMNI2 data, to provide the best accordance with the time coverage of the IPS data collecting periods. Figure 21 shows that the uncertainties of the IPS model and the OMNI2 data at the equator are of similar order. For the whole period 1985-2020, the mean relative uncertainties for the OMNI series and the model are 8% and 7%, respectively, with a ∼ ±3% scatter. At the beginning of the IPS data series, the relative uncertainty of the IPS model is systematically lower. Next, a slow increase of the relative uncertainty to the level of OMNI2 appears, to finally exceed its value for the most recent years.
On the other hand, the OMNI2 uncertainty seems to drop slightly with time. It is noticeable that the relative uncertainty is the largest in 1994, when the ∆v IPS −OMNI magnitude is also the largest. For the central period of the IPS data, the IPS and OMNI2 uncertainties are at the same level. An abrupt increase of the relative uncertainty for the IPS model is visible since 2008. Despite secular changes of the uncertainties, their distributions for the IPS model and the OMNI2 data are similar according to Pearson's χ 2 test and are not correlated. So the uncertainties may be treated as independent. Now we study the significance of the −∆v IPS −OMNI adjustments relative to the effective uncertainty of the IPS and the OMNI2 speed components used to calculate the adjustments (see Fig. 22). The comparison shows that the significant −∆v IPS −OMNI adjustments are obtained for the following years: 2014, 2019 and 2020, i.e. for 9% of all processed years (see Figure 22). Summarizing this part, we see that the scatter in the OMNI2 data is similar to that of the SW modeled speed, and for both data sets the relative errors are not correlated. Also, if a time trend exists in the IPS data, which at the moment is in question, it lies within the statistical accuracy of the model, so it should not introduce any systematic errors to the results. Taking the above into account, and having in mind that for the equator the yearly averaged OMNI2 and IPS SW speed models have the same statistical properties, and the adjustments are not statistically significant for 91% of the years in comparison with the model accuracy, the shifting of the helioatitude profiles of the model by −∆v IPS −OMNI seems appropriate. Statistical properties of the OMNI2 and the SW speed model uncertainty confirm that no significant systematic errors will be introduced this way. An additional argument in favor of application of the adjustments is that the OMNI2 speeds are derived from solar wind speeds measured at various spacecraft and cross-calibrated, so systematic differences resulting from application of different instrument and data analysis techniques are largely eliminated, and application of the adjustments to the IPS-derived speed profiles provides a good correspondence of the IPS SW 3D speed model with the OMNI2 data in the ecliptic plane.
Adoption of a different approach, i.e. adopting the hypothesis that a time trend in ∆v IPS −OMNI indeed exists, would cause nasty practical problems. The IPS observations suitable for our analysis arrive at a yearly cadence. An analytic formula would need to be fitted to the ∆v IPS −OMNI (t) time series. With each new year of data, the entire set of corrections to the yearly SW speed profiles would need to be applied, which would result in a new version of the entire data set starting from 1985. Given the fact that none of the simple analytic models we have tested shows a clear advantage, we decided to apply individual yearly corrections. The overall characteristics of the modeled SW speed is very similar to the previous model, which is shown in Fig. 23, but some differences are evident. Generally, while the WawHelIon 3DSW model tends to provide higher speed values during the solar minima at higher heliolatitudes, during the solar maxima the models are very consistent. This is seen at the beginning of the data set, as well as around 1994, when an exceptionally high difference between the OMNI2 speeds and those from IPS observations are present. In 1994, a comparison between the new and the previous model shows a kind of new feature for this year in the form of a deficit in the previous model, distinguishable especially in the north hemisphere. This feature is undoubtedly linked with a large difference between the OMNI2 and the IPS data for 1994, which was not present before, and should be monitored in future releases of the IPS data. Also, around 2008 and 2018 the WawHelIon 3DSW model provides systematically higher speeds, but for these years the difference in speed between the two models is smaller than around year 1994.
Looking at the dynamics of the SW speed, for the beginning of Solar Cycle (SC) 24 the WawHelIon 3DSW model provides an earlier decline of the SW speed. It is noticeable that this appears around year 2010. The difference in the dynamics is much larger than at the beginning of SC23 (between 1995 and 2000), when the decline started almost simultaneously in the previous and current models. Despite the earlier decline, the SW speed becomes consistent between the two models during the solar maximum.
We also compare the 3D IPS SW speed model with the in situ measurements from fast Ulysses scans. The Ulysses data are averaged in 5 • equatorial bands. The uncertainties were derived from the scatter of the data within the bins. The comparison is presented in Figure 24. A general good agreement with in situ measurements is seen during the first and the third fast scans. During second fast scan the Ulysses data show a large scatter due to physical changes, but the overall characteristics seems to agree.

Comparison of the densities
We point out an important difference between our approach to the calculation of the SW density and that used by Sokół et al. (2020). While we derive the density model for a given year from the corresponding mean yearly speed profile and the latitudinally invariant mean yearly solar wind energy flux, Sokół et al. (2020) use running yearly averages separately for each Carrington rotation of the invariant to calculate the Carrington period density profiles. The Carrington averages of the SW energy flux feature a larger scatter than the yearly averages, and these fluctuations are propagated into all heliolatitudes. As a result, the speed model by Sokół et al. (2020) features smooth variations of speed at heliolatitudes outside the equatorial band, in the case of the density, these variations are much more rapid. We are not convinced these variations are realistic since we are not sure if the invariance of the SW energy flux is maintained on such a short time scale as one solar rotation period. Therefore, we decided to use the approach presented above. After all, the existence of this invariant was based on pole-to-pole observations from Ulysses, as pointed out by Le Chat et al. (2012), and the time of Ulysses fast latitude scan between the poles was approximately 6 Carrington rotations (see, e.g., Figure 1 in Bzowski et al. 2014). The invariance at shorter time scales may exist but to our knowledge, has not been confirmed. Nevertheless, a comparison between the new and old models shows general similarities. Typical density variations correlated with the solar maxima and minima are present. Also, the north-south asymmetry in the density, as in the previous model, is present. Not all two-peak structures in time, present in the previous model, are confirmed by the new model. This may be due to the different approach to the calculation of the density that we used. Our model confirms the presence of the two-peak structure at the northern hemisphere in the middle latitudes for SC 22 and 24, and at high latitudes only for SC24. At the southern hemisphere, only in middle-high latitudes the two-peak structure appears during SC23.
A comparison of the absolute values of the densities between the two models shows that systematically higher values are obtained only for the highest latitudes in both hemispheres in SC24. The density values at the peak during SC24 at the highest latitudes in the northern hemisphere, returned by the WTPM 3DSW model, are the highest in the entire interval, reaching above 8 cm −3 . Also, in the southern hemisphere the WTPM 3DSW model provides high and more stable density values, while the previous model provides a lower density between the two-peak structures, i.e., around 1992 and 2014. It is likely that the lower value in the two-peak structures may be linked with the different approach to averaging and interpolating of the models.
A portion of the differences between the densities from the two models is due to the differences in the speed. They are propagated to the densities by the relation in Equation 3. Another portion of the differences comes up because of the different approach to the application of the invariant.

Concluding discussion
The agreement between WTPM 3DSW and the model by Sokół et al. (2020) in both speed and density are quite good, even though the original IPS data had been modified between these two releases of the model, as well as data filtering method and the analytic form of the fitting functions were changed. This illustrates the maturity and robustness of the line of modeling the 3D structure of the solar wind based on IPS CAT analysis and the approximations used in the papers by Sokół et al. (2013), Bzowski et al. (2013), Sokół et al. (2015b), Sokół et al. (2019a), and the present paper.
Analysis of the helioglow aimed at retrieval of the latitudinal structure of the solar wind, performed based on observations of the helioglow from the SWAN experiment onboard the SOHO mission (Bertaux et al. 1995) and presented in a series of papers by Lallement et al. (2010), Katushkina et al. (2013), Katushkina et al. (2019), Koutroumpa et al. (2019) suggested that the solar wind flux may feature maxima at midlatitudes. Our analysis of the SW speed structure, presented in figures 12-14, shows that while indeed, local minima and maxima of speed exist in different years, usually they are not pronounced, they are not regular and they occur at different latitudes. This may be an important finding because, unlike the previous model, WawHelIon 3DSW is able to reproduce any profile that is supported by the data. Hence, should the features suggested by these authors be real, they would be due to the evolution of the density with heliolatitude, but then very likely the reality of the SW energy flux invariance with heliolatitude would be challenged. This question certainly requires further study.

SUMMARY AND CONCLUSION
We present a new model of evolution of the latitudinal structure of the SW speed and density, calculated from a new IPS-derived data set of SW speed: the WawHelIon 3DSW model. The input data cover a newlyreprocessed entire set of available IPS observations subjected to CAT analysis, with an important addition of variable coefficient α, connecting the magnitude of the fluctuations of the SW electron density ∆n e with the SW speed V: ∆n e ∼ V α S W . The model covers a time interval from 1985 until the end of 2020. We developed and applied a new method of filtering of the IPS-derived Carrington maps of solar wind speed against outlying values. This method acknowledges that some of the data points in the Carrington SW speed maps may be biased. It identifies the suspect outlier points based on the ESD method without arbitrary assumptions and returns a filtered data set, suitable for approximation by an analytic function. We calculated yearly-averaged heliolatitudinal profiles of SW speed by fitting the filtered yearly data with an analytic formula based on Legendre polynomials of an order that can vary between the years, and with an additional condition of a null first derivative over the poles. The order of the formula for individual years was determined based on statistical analysis of the residuals. The formula is flexible enough to be able to reproduce even small-scale variations, should they be identified, and does not overfit the data.
We performed an extensive comparison of the predictions of the model with in-situ SW time series of the SW speed from the OMNI2 collection and requested that the two data sets agree at the solar equatorial plane. To that end, we studied differences between the OMNI2 and our model predictions and derived necessary corrections based on a meticulous statistical analysis.
We calculated yearly heliolatitudinal profiles of the SW density using yearly averaged corrected SW speed profiles assuming that the SW energy flux is latitudinally invariant on a time scale of one year. Finally, we projected the SW speed and density models on a regular heliolatitudinal grid, with time nodes at halves of individual Carrington rotation periods. The equatorial band was replaced with appropriate Carrington period averages of the SW speed and density from the OMNI2 time series. The newly developed model can be used in its tabular form, as a grid of densities and speeds tabulated in heliolatitude and time, or in its analytic form for individual years.
The results of the WTPM 3DSW model and that of Sokół et al. (2020) generally agree, even though the latter one used the former version of the input Carrington maps of SW speed, a different method of data filtering, and a different analytic approximating function. Also different was the assumptions made on the invariance of the solar wind energy flux. The WTPM 3DSW model suggests higher SW speed values during solar minima and a faster decline of SW speed at the beginning of SC24. The higher speed values at solar minima may be caused by the elimination of a bias in Carrington maps caused by outlier data in all CRs during all years, which was not fully addressed in Sokół et al. (2020).
The SW densities provided by the WTPM 3DSW model are also in general agreement with that by Sokół et al. (2020), but a part of fine features seen in the model by Sokół et al. (2020) as two-peak structures disappeared in our model. This is likely due to the different use of the solar energy flux invariant. Since it is currently unknown whether or not the invariance of the solar wind energy flux holds on a time scale of one solar rotation, as assumed by Sokół et al. (2020), or a time scale of 1 year, as assumed in the present paper, this discrepancy provides an estimate of the uncertainty of the density of SW at high latitudes. The general agreement between the two models suggest the robustness of data analysis and methodology adopted at all phases of research in the line of papers using the IPS observations and analytic function fitting that gives the time evolution of the SW structure.
The WawHelIon 3DSW model of the evolution of SW parameters is well suitable for use in global heliospheric studies, as well as analysis of heliospheric in-situ measurements. In particular, it can be used to compute attenuation of energetic neutral atoms between their origin in the outer heliosheath and detectors located deep in the inner heliosphere, like IBEX and the planned IMAP missions at 1 au.