Evolution of the Solar Lyα Line Profile during the Solar Cycle

Our model is based on observations of disk-integrated profiles of the solar Lyα line measured by L15 using SUMER on board SOHO.¹ These profiles are shown collectively in Figure 1. There are 43 profiles observed between 1996 May and 2009 April. The times of observations are illustrated in Figure 2, where the daily values of the wavelength-integrated irradiance for the observation days are superimposed on daily irradiance values from the composite Lyα series from LASP (Woods et al. 2005). As is evident, approximately half of the observations were taken during low solar activity, while a little less than half were taken during maximum solar activity in 2000–2003.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Details of the observations and processing applied to obtain the profiles scaled in the absolute units are given in L15 and in the earlier papers from those authors; therefore, here we only point out some important aspects. The SUMER instrument had not been designed to provide disk-averaged profiles and therefore a specially devised observation scheme suggested by J. L. Bertaux had to be used. Effectively, the instrument measured the profile of the solar Lyα radiation scattered inside the instrument optics when the boresight was directed ∼2000'' off the disk center. To safeguard a more or less uniform contribution from all parts of the solar disk, observations from several pointings around the solar disk were averaged. As the angular distances of the boresight during a series of observations aimed at obtaining an individual profile for a given day were sometimes different, the spectral intensities from those individual pointings had to be rescaled to a common offset distance. The measurement process for one composite profile took more than four hours.

The resulting profile is composed of 281 points at wavelengths equally spaced with a 0.001 nm pitch. During the observations, the profiles had been sampled at several wavelengths and the observations were subsequently deconvolved to eliminate the effect of the instrument point-spread function. The final data product has a much higher resolution than the spectral resolution of the instrument and does not show any significant point-to-point scatter. The absolute scaling of the profiles was performed by the observers by requiring that the integral of the observed profile within the boundaries ±0.15 nm around the nominal wavelength of the line center is equal to the daily value of the composite series of the Lyα irradiance for a given day of observation (see Table A.1 from L15). The uncertainties in the data set are complex. The uncertainty of individual profiles are given by the authors of the observations at 15%. It seems, however, that an additional source of uncertainty is the systematic uncertainty of the LASP composite irradiance, which is also ±15% (Woods et al. 2005). This uncertainty results in systematic upward or downward offset of the entire data product, and it will result in a systematic uncertainty of the radiation pressure for H and D by the identical percentage. Additionally, because the integration time used to obtain the measurement of the LASP daily datum and that used during the spectral observations of the profiles were different, and because of the inevitable measurement scatter in the LASP observations, the absolute scaling of each of the individual profiles is biased relative to the other profiles from the sample. This results in random upward or downward shifts of the profiles relative to each other. Effectively, this uncertainty introduces correlations to the entire data set. The data points pertaining to individual profiles are more strongly correlated than data points from different profiles, and absolute accuracies of different profiles are different. Therefore, it is not reasonable to expect that the final model of the evolution of the profile during the solar cycle will fit all of the observed profiles with the same accuracy.

The composite time series of the solar Lyα irradiance was recapitulated by Bzowski et al. (2013). The total irradiance, which we denote ${I}_{\mathrm{tot}}$, is based on direct daily observations, with occasional gaps filled with values derived from proxies. The absolute calibration is based on UARS/SOLSTICE² (Woods et al. 1996, 2000). Observations from later experiments, TIMED and SORCE, are scaled to this calibration. The irradiance is integrated over the 1 nm interval from 121 to 122 nm. The series cover the time interval since 1947 until the present. The ${I}_{\mathrm{tot}}$ values for times before 1979 are obtained from scaling the solar F10.7 radio flux (Tapping 2013). More recently, gaps have been filled mostly using the Mg ii core to wing index (Viereck & Puga 1999). The use of proxies adds additional intermittent increase in the uncertainty of the daily irradiance values.

The values of the composite solar Lyα irradiance from LASP are listed at a daily cadence but they are not daily averages. Actually, the measurements are taken during a fraction of the day. Also, the data used for proxies are not daily average values, even though they are listed at a daily cadence. Therefore, the daily ${I}_{\mathrm{tot}}$ values used by L15 to scale the Lyα profiles do not precisely correspond to the times of profile observations. The daily values of ${I}_{\mathrm{tot}}$ directly from the LASP are shown as gray dots in Figure 2, and the Carrington averages by the blue line. The daily values of ${I}_{\mathrm{tot}}$ for the dates corresponding to profiles observed by L15 are shown with black dots, and the subset from Lemaire et al. (2002) used by TB09 are taken in black circles.

For modeling radiation pressure for neutral H and D atoms in the heliosphere at an arbitrary longitude and latitude in the heliosphere and an arbitrary time t, we recommend to use the Lyα total irradiance averaged over Carrington rotations, linearly interpolated to t, and not the corresponding daily value of ${I}_{\mathrm{tot}}$ extracted from the daily time series from LASP. The reason for this is the following. For a given date, the value of ${I}_{\mathrm{tot}}$ from the LASP time series is only valid for the geometric location of the Earth. The instantaneous irradiance for other longitudes is unknown. Using the daily values directly from the LASP time series would imply that the entire Sun is fluctuating in its Lyα output precisely as measured by LASP. However, from the fact that the measured daily ${I}_{\mathrm{tot}}$ values feature quasi-periodic variation during the Carrington period (see the right panel of Figure 3), it follows this is not the case: most of those variations are due to active regions, which are visible from the observer's hemisphere, but are invisible from (and thus cannot affect) the opposite hemisphere. These active regions change in strength and gradually fade out, while other ones appear in the solar disk at various heliolongitudes and heliolatitudes. Therefore, averaging over solar rotation period seems to provide a good balance between the true variations in the solar output and our ignorance of the solar radiation field in the heliosphere away from the Earth.

Zoom In Zoom Out Reset image size

The relation of the short-time variations in ${I}_{\mathrm{tot}}$ to the distribution of active regions, to their locations structured in heliolatitude, one may suspect that radiation pressure will be a function of heliolatitude. This aspect is still poorly investigated, as discussed by Bzowski et al. (2013). However, observations of the solar corona in EUV lines (Auchère 2005) and of the heliospheric backscatter glow (Pryor et al. 1992), as well as some theoretical considerations (Cook et al. 1981), suggest that ${I}_{\mathrm{tot}}$ is a mild function of heliolatitude ϕ. Therefore, we follow the recommendation by Bzowski (2008) and adopt the monthly averaged ${I}_{\mathrm{tot}}$ value as the following function of heliolatitude ϕ:

$\begin{eqnarray}&&{I}_{\mathrm{tot}}(\phi )={I}_{\mathrm{tot}}(0)\sqrt{{a}_{\mathrm{Lya}}{\sin }^{2}\phi +{\cos }^{2}\phi }\end{eqnarray} \tag{ 2 }$

with ${a}_{\mathrm{Lya}}=0.8$.

In this paper, we propose a more physical approach to the problem of defining an analytic approximation of the line profile. The shape of the line is determined by the velocity distribution of H atoms that are excited to the second energy level. Recent studies (Jeffrey et al. 2017) have shown that in the absence of thermodynamical equilibrium, the atom velocity distribution can be described by the kappa distribution and the resulting line profiles are also kappa-like. Hence, in our model we used a kappa function to reproduce the main shape of the Lyα emission line

$\begin{eqnarray}&&{F}_{K}={A}_{K}{\left[1+\displaystyle \frac{{({v}_{r}-{\mu }_{K})}^{2}}{2{\sigma }_{K}^{2}\kappa }\right]}^{-\kappa -1}.\end{eqnarray} \tag{ 8 }$

Here, ${\mu }_{K}$ represents the radial velocity corresponding to the central wavelength of the line, which we allow to differ from the nominal wavelength of 121.567 nm, κ is the κ parameter of the kappa-like line profile, and ${\sigma }_{K}$ represents the overall width of the profile.

The core of the Lyα line is affected by the absorption in the solar transition region. To include this effect in our model, we added a Gaussian component

$\begin{eqnarray}&&{F}_{R}=\displaystyle \frac{{A}_{R}}{{\sigma }_{R}\sqrt{2\pi }}\exp \left[-\displaystyle \frac{{({v}_{r}-({\mu }_{K}+d\mu ))}^{2}}{2{\sigma }_{R}^{2}}\right]\end{eqnarray} \tag{ 9 }$

centered at the wavelength corresponding to radial velocity ${\mu }_{K}+d\mu$, i.e., we allow the central wavelength of the Lyα line and that of the central reversal to vary by $d\mu$. The width of the central reversal is given by ${\sigma }_{R}$.

Finally, we took into account the fact that the solar spectrum is not flat and we included a linear spectral background

$\begin{eqnarray}&&{F}_{\mathrm{bkg}}={a}_{\mathrm{bkg}}{v}_{r}+{b}_{\mathrm{bkg}}.\end{eqnarray} \tag{ 10 }$

The final form of the function is given by

$\begin{eqnarray}&&{F}_{\mathrm{model}}={F}_{K}-{F}_{R}+{F}_{\mathrm{bkg}}.\end{eqnarray} \tag{ 11 }$

It is schematically presented in Figure 4 along with its components. The relation between the wavelength in the absolute units and the radial velocity is given by

$\begin{eqnarray}&&{v}_{r}(\lambda )=-\displaystyle \frac{{\rm{\Delta }}\lambda }{{\lambda }_{0}}c,\end{eqnarray} \tag{ 12 }$

where ${\rm{\Delta }}\lambda =\lambda -{\lambda }_{0}$ is the value of the shift from the center of line (in the original observations), ${\lambda }_{0}=121.567\,\mathrm{nm}$ is the wavelength of the Lyα line, and c is the speed of the light. Note that negative radial velocities correspond to longer wavelengths. In consequence, the blue horn, which was located on the blue side (left part of Figure 1), is on the right-hand side of Figure 4, where we are using radial velocity units.

Zoom In Zoom Out Reset image size

With the model given by Equations (8) through (11), we were able to reproduce the characteristic shape of the line profile, with two peaks and the asymmetries that are visible in the observed profiles, especially the asymmetry between the heights of the two horns. Further refinement of the model is discussed in the following section.

We found that parameter values obtained in fits to individual profiles showed a systematic dependence on ${I}_{\mathrm{tot}}$, which had to be recognized and accounted for in the final model. Thus, we constructed the final model in two steps. First, we fitted the parameters of the model defined in Equations (8)–(11) individually to each of the 43 profiles. Then, we analyzed the relations between the fit parameters and ${I}_{\mathrm{tot}}$ and inscribed them into the model.

4.2.1. Fitting Individual Profiles

We performed the fitting after converting the original wavelength scale given in nanometers to the Doppler scale in km s⁻¹ using Equation (12), and the spectral irradiance originally given cm⁻² s⁻¹ nm⁻¹ to the dimensionless μ parameter, which represents the fraction of solar gravity compensated by radiation pressure, by dividing the spectral irradiance by the scaling factor $\mathrm{sclFctH}=3.34\times {10}^{12}$. The μ units are very convenient from the viewpoint of calculating the trajectories of H atoms in the heliosphere. For all profiles, we used for fitting the portions of the data corresponding to ±200 km s⁻¹ around 0. This was done to avoid biasing the results by the far wings of the profiles.

Each profile was fitted with the nine-parameter function described by Equations from (8) to (10) using the nonlinear least-squares method from Wolfram Mathematica (Wolfram-Research 2016). We obtained a very good match between the data and the fitted functions in most of the cases. Within ±200 km s⁻¹, only 5 out of 43 profiles have residuals (defined here as (data—model)/data) larger than 10%, and this happens only for v_r close to ±200 km s⁻¹. Generally, most of the residuals are less than 5%, as shown in Figure 5. Example profiles for the minimum and maximum of solar activity and the corresponding residuals are shown in Figure 6. Near the maximum of solar activity (2001 October 28), the magnitude of the profile was much higher than during solar minimum (1996 December 5) and consequently, the magnitude of radiation pressure was also larger. The residuals in both cases are safely within the ±5% margin for radial velocities ±100 km s⁻¹, i.e., within an interval the most important for ISN H in the heliosphere. These values are within the data uncertainty level given by L15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result of fitting individual profiles, we obtained 43 sets of parameters, each describing a different profile observed by the SUMER detector. Their magnitudes and uncertainties are shown in Figure 7 as a function of ${I}_{\mathrm{tot}}$.

Zoom In Zoom Out Reset image size

4.2.2. Parameter Correlation with ${I}_{\mathrm{tot}}$

The model parameter values, shown in Figure 7, clearly show a correlation with the line-integrated irradiance ${I}_{\mathrm{tot}}$, albeit with a large scatter. As the procedure being used to compute the density and higher moments of the distribution function of ISN H and D in the WTPM model (Tarnopolski & Bzowski 2008, 2009) must calculate thousands of trajectories of individual atoms, and each trajectory is calculated by numerically solving Equation (1) with the μ factor being a function of v_r and t, our goal was to create a model for $\mu ({v}_{r},{I}_{\mathrm{tot}}(t))$ that can be expressed by a relatively simple analytical formula. On the other hand, the model must be able to reproduce the true radiation pressure relatively accurately because the density and bulk velocity of ISN H and D are sensitive functions of radiation pressure.

TB09 used a simple proportionality scaling of the line profile by ${I}_{\mathrm{tot}}$, which was demonstrated by Schwadron et al. (2013) and Katushkina et al. (2015) to likely be too simple. With more than 40 sets of parameters for profiles observed during the entire solar cycle, we could search for relations between profile parameters corresponding to various levels of solar activity. We correlated every parameter of the model (A_K, ${\mu }_{K}$, ${\sigma }_{K}$, κ, A_R, $d\mu$, ${\sigma }_{R}$, a_bkg, and b_bkg) with the total irradiance in Lyα scaled by the solar cycle average value (computed over 23rd solar cycle for which observed profiles are available) ${I}_{\mathrm{tot}}/\langle {I}_{\mathrm{tot}}\rangle$. We found that the correlation of all of the parameters with ${I}_{\mathrm{tot}}$ can be approximated by a linear function, as shown in Figure 7.

We used this finding in the construction of the final model of $\mu ({v}_{r},{I}_{\mathrm{tot}}(t))$. We fitted each of the nine parameters of the model with the linear function using the least-squares method. A parameter P_i is expressed in following convenient form:

$\begin{eqnarray}&&{P}_{i}={\beta }_{i}\left(1+{\alpha }_{i}\displaystyle \frac{{I}_{\mathrm{tot}}}{\langle {I}_{\mathrm{tot}}\rangle }\right),\end{eqnarray} \tag{ 13 }$

where ${P}_{i}=({A}_{K},{\mu }_{K}$, ${\sigma }_{K},\kappa ,{A}_{R},d\mu$, $\sigma ,{b}_{\mathrm{bkg}},{a}_{\mathrm{bkg}})$, $\langle {I}_{\mathrm{tot}}\rangle \,=4.25\times {10}^{11}$ ph cm⁻² s⁻¹, and ${\alpha }_{i}$ and ${\beta }_{i}$ are the fit parameters, listed in Table 1. In each of those nine fits, the data were the 43 values of given parameters for the 43 Lyα profiles, and the independent variable was the corresponding daily ${I}_{\mathrm{tot}}/\langle {I}_{\mathrm{tot}}\rangle$ value. The data were weighted by inverse squares of the parameter errors, obtained from the fits presented in Section 4.2.1. Results of fitting the parameters from Equation (13) are listed in Table 1 and the fitted lines shown in Figure 7.

Table 1.  Coefficients of the Linear Correlations between the Model Parameters for H and Total Irradiance in Lyα, Defined in Equation (13)

Parameter (Pi)	${\alpha }_{i}$	${\beta }_{i}$
AK	4.873	1.10341
${\mu }_{K}$	4.349	−1.1064
${\sigma }_{K}$	38.251	0.0969713
κ	2.074	−0.273553
AR	446.14	0.606729
$d\mu$	−0.285	−0.804873
${\sigma }_{R}$	32.224	−0.0431647
bbkg	0.026	0.503516
abkg	$0.435\times {10}^{-4}$	−1.3148

Note.The model along with all parameters is available online: http://users.cbk.waw.pl/~ikowalska/index.php?content=lya.

Download table as:  ASCIITypeset image

Instead of the composite total irradiance in Lyα, one could potentially use the solar radio flux F10.7. These two quantities are strongly correlated, and the first one is calculated from the latter when there are no direct observations (Woods et al. 1996). We performed the whole model fitting procedure using the radio flux F10.7 instead of ${I}_{\mathrm{tot}}$, but the results were a little less accurate than in the case of total irradiance in Lyα.

The final result for radiation pressure model for neutral H in the heliosphere is given by Equations (8)–(11) and (13), summarized in Equation (14), with the coefficients listed in Table 1. The free parameters are the total irradiance in the Lyα line ${I}_{\mathrm{tot}}$ in ph cm⁻² s⁻¹ and radial velocity v_r in km s⁻¹. The result is in the units of compensation of the solar gravity force (the μ-units).

Summarizing, the model of radiation pressure $\mu (r,{v}_{r},t,\phi )$ for radial velocity v_r (in km s⁻¹), time t, heliolatitude ϕ, heliocentric distance r is given by the following set of equations

$\begin{eqnarray}{I}_{\mathrm{tot}}(t,0) & = & {I}_{\mathrm{tot},n}+\displaystyle \frac{{I}_{\mathrm{tot},n+1}-{I}_{\mathrm{tot},n}}{{t}_{n+1}-{t}_{n}}(t-{t}_{n})\\ {I}_{t,\phi } & = & {I}_{\mathrm{tot}}(t,0)\sqrt{{a}_{\mathrm{Lya}}{\sin }^{2}\phi +{\cos }^{2}\phi }\\ {P}_{i} & = & {\beta }_{i}\left(1+{\alpha }_{i}\displaystyle \frac{{I}_{t,\phi }}{\langle {I}_{\mathrm{tot}}\rangle }\right)\\ {v}_{r} & = & -\displaystyle \frac{\lambda -{\lambda }_{0}}{{\lambda }_{0}}c\\ {F}_{K} & = & {A}_{K}{\left[1+\displaystyle \frac{{({v}_{r}-{\mu }_{K})}^{2}}{2{\sigma }_{K}^{2}\kappa }\right]}^{-\kappa -1}\\ {F}_{R} & = & \displaystyle \frac{{A}_{R}}{{\sigma }_{R}\sqrt{2\pi }}\exp \left[-\displaystyle \frac{{({v}_{r}-({\mu }_{K}+d\mu ))}^{2}}{2{\sigma }_{R}^{2}}\right]\\ {F}_{\mathrm{bkg}} & = & {a}_{\mathrm{bkg}}{v}_{r}+{b}_{\mathrm{bkg}}\\ \mu (r,{v}_{r},t,\phi ) & = & ({F}_{K}-{F}_{R}+{F}_{\mathrm{bkg}}){({r}_{E}/r)}^{2},\end{eqnarray} \tag{ 14 }$

where ${t}_{n},{t}_{n+1}$ are the times of the starts of the Carrington rotations immediately preceding and following the time of calculation t, ${a}_{\mathrm{Lya}}=0.8$, ${I}_{\mathrm{tot},n};$ ${I}_{\mathrm{tot},n+1}$ are the values of the LASP time series averaged over the nth and $n+1$ th Carrington periods, respectively; $\langle {I}_{\mathrm{tot}}\rangle =4.25\times {10}^{11}$ ph cm⁻² s⁻¹; and ${P}_{i}=({A}_{K},{\mu }_{K},{\sigma }_{K},\kappa ,{A}_{R},d\mu ,\sigma ,{b}_{\mathrm{bkg}},{a}_{\mathrm{bkg}})$ are the parameters of the profile, obtained from Table 1 for H and Table 2 for D.

Table 2.  Coefficients of the Linear Correlations between the Model Parameters for D and Total Irradiance in Lyα, Defined in Equation (13)

Parameter (Pi)	${\alpha }_{i}$	${\beta }_{i}$
AK	2.438	1.10341
${\mu }_{K}$	−77.031	0.196728
${\sigma }_{K}$	38.251	0.0969713
κ	2.074	−0.273553
AR	223.24	0.606729
$d\mu$	−0.285	−0.804873
${\sigma }_{R}$	32.224	−0.0431647
bbkg	0.013	0.503516
abkg	$0.218\times {10}^{-4}$	−1.3148

Note. The model along with all parameters is available online: http://users.cbk.waw.pl/~ikowalska/index.php?content=lya.

Download table as:  ASCIITypeset image

The same formulae can be used to model radiation pressure for neutral D. This is a result of a simple shift in the wavelength due to isotope effect (${v}_{r}\to {v}_{r}+81.3802$ km s⁻¹) and scaling down by the H/D mass ratio, equal to 0.5003. The coefficients thus obtained are given in Table 2 and example profile (2001 October 28) presented in Figure 8 along with a H profile for comparison.

Zoom In Zoom Out Reset image size

The data from L15 offer an opportunity to study the variation of the full-disk profiles within a solar rotation. We have identified two time intervals when at least three profiles are available within a time interval less than the Carrington rotation period: one during low solar activity (in 1997 December, with five profiles), and the other one during high solar activity (in 2001 August, three profiles). These profiles are shown in Figure 3 along with our model profile, calculated for ${I}_{\mathrm{tot}}$ equal to the mean value over the 27.4 day interval centered at the arithmetic mean of the dates from the data corresponding profile subsets.

The full-disk profiles vary considerably practically overnight, both during high and, surprisingly, low solar activity. These variations are reflected in the total irradiance ${I}_{\mathrm{tot}}$. This finding illustrates why it is not realistic to expect that a profile measured for any individual date will be representative for the entire heliolongitude range. It also illustrates that adopting a model where the value of ${I}_{\mathrm{tot}}$, which is the driver for the model, obtained from interpolation of Carrington period averages of the daily LASP series is a reasonable choice. The scatter between the profiles observed within a few days from each other may be adopted as a measure of the uncertainty of our model resulting from adopting average profiles characteristic for a the entire range of heliolongitudes.

The quality of the approximation we propose is further illustrated in Figure 11. We show that because, in our approach, the sole model parameter is the disk-integrated solar irradiance ${I}_{\mathrm{tot}}$, it is reasonable to expect that average profiles, obtained from measurements with similar values of ${I}_{\mathrm{tot}}$, will be well reproduced by the model evaluated with ${I}_{\mathrm{tot}}$ equal to the mean value of this parameter for all profiles within such a group. To reduce statistical scatter, we decided to split the observed data set into three classes according to ${I}_{\mathrm{tot}}$ values: low activity (${I}_{\mathrm{tot}}\,\leqslant 4\times {10}^{11}$ cm⁻² s⁻¹), medium activity ($4\lt \times {10}^{11}\lt {I}_{\mathrm{tot}}\leqslant 5\,\times {10}^{11}$ cm⁻² s⁻¹), and high activity (${I}_{\mathrm{tot}}\gt 5\times {10}^{11}$ cm⁻² s⁻¹), and compare the model profiles with the observed profiles averaged over the three subsets of ${I}_{\mathrm{tot}}$. We also show profiles that are obtained for the aforementioned three subsets from averaging the individually fit profiles, discussed in Section 4.2.1. The average individually fit profiles agree very well with the average observed profiles. The agreement between the averaged observed profiles and the profiles obtained from the final model is also good, but a little worse than for individually fit profiles. This is because in the final model, we use the linear dependence of the fit parameters on the total irradiance.

One of important aspects of the Lyα line profile is the depth of the central reversal and the height of the two horns. It can be represented as ratios of the spectral irradiance for one of the horns to the local minimum in the spectral irradiance, located between the two horns. Schwadron et al. (2013) and Katushkina et al. (2015) suggested that the ratio of the blue horn to the reversal depth needs to be larger than that in the model by TB09 to reproduce the spectrum of interstellar neutral H flux observed by IBEX-Lo.

The horn-to-minimum ratios for all observed profiles are shown in Figure 9, separately for the blue (positive velocities) and red (negative velocities) horns, as a function of the total irradiance ${I}_{\mathrm{tot}}$. Clearly, this ratio drops with the increase of solar activity approximately linearly, i.e., the depth of the central reversal approximately linearly decreases with the increase of solar activity. This can be understood as on one hand, Tian et al. (2009b) showed that the Lyα profiles observed in active regions do not feature considerable self-reversal, and on the other hand, the number of active regions visible on the solar disk increases with the increase of solar activity, so their contribution to the disk-averaged profiles should be larger than during low solar activity.

Zoom In Zoom Out Reset image size

The present model does show this dependence of the horn-to-minimum ratios on the total irradiance. The magnitude of the model horn-to-minimum ratios are in agreement with those observed, for both the red and blue horns, and as well for the models fit to individual profiles as for the global model presented in Section 7, where the model parameters are linear functions of ${I}_{\mathrm{tot}}$. This also holds for the model and data profiles averaged within the three aforementioned data subsets.

Another aspect is the relation between ${I}_{\mathrm{tot}}$ and the level of radiation pressure at the line center: ${\mu }_{0}({I}_{\mathrm{tot}})=\mu ({v}_{r}=0,{I}_{\mathrm{tot}})/{I}_{\mathrm{tot}}$. Historically, this is an important parameter because it was used to assess the radiation pressure from the wavelength-integrated intensity in the models where the dependence of radiation pressure on v_r was neglected.

In the past, this ratio was modeled using a nonlinear relation. Vidal-Madjar & Phissamay (1980) suggested that the ratio between the spectral flux and line-integrated intensity is ${f}_{0}=0.54\times {I}_{\mathrm{tot}}^{1.54}$. Emerich et al. (2005) modified this relation to ${f}_{0}=0.64\times {I}_{\mathrm{tot}}^{1.21}$, and Lemaire et al. (2005) to ${f}_{0}=0.62\times {I}_{\mathrm{tot}}^{1.23}$. In all of these formulae, ${I}_{\mathrm{tot}}$ was to be used in the units 10¹¹ cm⁻² s⁻¹. While clearly nonlinear in general, for the values of ${I}_{\mathrm{tot}}$ characteristic for the Sun, these relations are linear in a very good approximation, which was pointed out by Bzowski et al. (2013) and Lemaire et al. (2015).

In our model, the quantity ${\mu }_{0}({I}_{\mathrm{tot}})$ can be calculated analytically

$\begin{eqnarray}{\mu }_{0} & = & {b}_{\mathrm{bkg}}+{A}_{K}{\left(1+\displaystyle \frac{{\mu }_{K}^{2}}{2\kappa {\sigma }_{K}^{2}}\right)}^{-1-\kappa }\\ & & -\displaystyle \frac{{A}_{R}}{\sqrt{2\pi }{\sigma }_{R}}\exp \left[-\displaystyle \frac{{({\mu }_{K}+d\mu )}^{2}}{2{\sigma }_{R}}\right],\end{eqnarray} \tag{ 15 }$

where the coefficients must be calculated for a given value of ${I}_{\mathrm{tot}}$ from Equation (13) using the appropriate parameter values from Table 1. A plot of this formula is shown in Figure 10, along with the ${\mu }_{0}$ ratios taken directly from the data and with a linear fit to Equation (15). Clearly, the linear approximation for the dependence of the radiation pressure at the line center on the total irradiance is a very good one, both for the data and for our model. Our model systematically underestimates the radiation pressure in comparison with the observed profiles by ∼9% for the lowest values of ${I}_{\mathrm{tot}}$ and by ∼3% for ${I}_{\mathrm{tot}}\gt 5.2\times {10}^{11}$ cm⁻² s⁻¹. Therefore, for the models neglecting the dependence of radiation pressure on radial speed of H atoms, we recommend using the appropriately scaled formula from L15. We point out, however, that ISN H atoms in the heliosphere spend very little time at ${v}_{r}\simeq 0$, so this systematic difference is not expected to significantly affect heliospheric models. Note that this systematic difference vanishes for larger radial speed, as can be inferred from the accuracy of the reproduction of the the horn/minimum ratio (see Figure 9).

Zoom In Zoom Out Reset image size

The present model differs from TB09 in some important aspects. The model function that we used is much more sophisticated than the simple three-Gaussian model used by TB09, and the quality of data reproduction that we have obtained is generally better. This is understandable, as the present model was obtained from many more profiles than were available to TB09. The differences between the two models are illustrated in Figure 11.

Zoom In Zoom Out Reset image size

Due to the definition of the model function, TB09 has the horn-to-minimum ratios almost independent of the total irradiance ${I}_{\mathrm{tot}}$. These ratios are illustrated in Figure 9. Clearly, TB09 underestimates the horn-to-minimum ratios for the solar minimum conditions, as suggested by Schwadron et al. (2013) and Katushkina et al. (2015). It is not clear at the moment if this profile change will be sufficient to understand the spectrum of ISN H observed by IBEX, but definitely, the model we have derived agrees in this aspect with the suggestions in those papers.

Effects of the new model of radiation pressure as a function of atom radial velocities are likely the most important for ISN H and D distributions, especially for the low solar activity conditions. The radiation pressure dependence on radial velocity is stronger than in TB09. This will likely affect estimates of the PUIs and solar wind ENA production, as well as the heliospheric backscatter Lyα glow. On the other hand, effects for ENA survival probabilities will likely be small, because ENAs due to their large velocity are mostly sensitive to the regions of the profile outside the self-reversal.

The effect for the distribution of ISN D in the inner heliosphere needs careful evaluation, but we believe it will likely be smaller than for ISN H, because ISN D is mostly sensitive to the blue horn region of the Lyα profile, and in this region, as can be seen in Figure 11, the magnitudes of the μ factor returned by the present model and that of TB09 are similar.

This discussion certainly must be regarded as speculative. Effects of the new model of radiation pressure on the models of heliospheric ISN H and D and some of their derivative populations (PUIs, solar wind ENAs, direct samples collected by IBEX-Lo) and the helioglow will be a subject of a future paper.

The uncertainties of our model are a result of superposition of many different errors and uncertainties. They include the accuracy of the absolute calibration of the total irradiance ${I}_{\mathrm{tot}}$, the accuracy of the measurements by the SUMER detector and of the deconvolution process used by L15, as well as fit errors obtained in our two fitting procedures. We have no knowledge of some of those uncertainties; therefore, it is almost impossible to estimate the accuracy of our model in a formal way.

L15 suggest the total uncertainties of their profiles are at a level ∼15%. The residuals in our individual profile fits are much less than that (see Figure 5). However, profiles obtained for specific days are not fully representative for the purpose of modeling the solar radiation pressure for the needs of heliospheric studies.

The uncertainties related to this aspect may be estimated from discrepancies between our mean profile and profiles measured by L15 within several days (less than one Carrington period), as illustrated in Figure 3. They seem to be at a level of ∼15% for locations close to solar equator. At larger heliolatitudes, there is additional uncertainty related to the heliolatitude dependence of the total irradiance ${I}_{\mathrm{tot}}$, which currently cannot be verified experimentally, but is believed to be another ∼15% (comparable to the magnitude of the departure of the pole-to-equator ratio from 1).

An uncertainty that can be quantified is the influence of the uncertainty in the absolute calibration of ${I}_{\mathrm{tot}}$. This uncertainty is on the order of 15% (Woods et al. 2005). Because the mathematical form of our model (Equation (14)) is strongly nonlinear in ${I}_{\mathrm{tot}}$, it follows that if the ${I}_{\mathrm{tot}}$ time series used to find the coefficients of our model is biased by a factor $1\pm \delta$, $0\leqslant \delta \ll 1$, then the entire model derivation procedure must be repeated with ${I}_{\mathrm{tot}}$ values rescaled by a $1\pm \delta$; however, the data would also need to be identically rescaled. Therefore, the shapes of the profiles (both data and model) will not change, only their magnitudes will be modified by a factor $1\pm \delta$.

We verified this by looking for differences between the radiation pressure models obtained as a result of the full fitting procedure with the bias factors $1\pm \delta ,\delta =0.15$ applied to both ${I}_{\mathrm{tot}}$ and data, and the model resulting from simple rescaling of the nominal model, evaluated with the nominal (non-scaled) value of ${I}_{\mathrm{tot}}$, by the factor $1\pm \delta$. We found no differences. Therefore, if a systematic rescaling of the ${I}_{\mathrm{tot}}$ by a factor of $1\pm \delta$ is needed in the future, then the only change of the radiation pressure in our model will be multiplication by $1\pm \delta$ and using the old (i.e., non-scaled) value of ${I}_{\mathrm{tot}}$.