A&A 432, 687-698 (2005)
DOI: 10.1051/0004-6361:20041606

Nonradial and nonpolytropic astrophysical outflows

VII. Fitting ULYSSES solar wind data during minimum

C. Sauty ¹ - J. J. G. Lima ^2,3 - N. Iro ⁴ - K. Tsinganos ⁵

1 - Observatoire de Paris, LUTH, 92190 Meudon Cedex, France
2 - Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal
3 - Departamento de Matemática Aplicada da Faculdade de Ciências, Universidade Porto, Rua do Campo Alegre, 657, 4150-007 Porto, Portugal
4 - Observatoire de Paris, LESIA, 92190 Meudon Cedex, France
5 - IASA and Section of Astrophysics, Astronomy & Mechanics, Department of Physics, University of Athens, Panepistimiopolis 157 84, Zografos, Greece

Received 7 July 2004 / Accepted 20 October 2004

Abstract
Exact axisymmetric analytical solutions of the governing MHD equations for magnetized and rotating outflows are applied to the solar wind during solar minimum as observed by ULYSSES. Using the spacecraft data, the latitudinal dependences of physical quantities such as the density, velocity, magnetic field and temperature are analytically described. The self-similar solutions are then compared to the global structure of the wind from one solar radius to 5 AU and beyond, including consistently the rotation of the outflow. The model makes it possible to describe the initial flaring of the magnetic dipolar structure, reproducing in a satisfactory way the observed profiles of the velocity, density and temperature with heliocentric distance. Finally, this model is in agreement with the conjecture that the solar wind should not be collimated at large distances, even close to its rotational axis.

Key words: magnetohydrodynamics (MHD) - solar wind - Sun: corona - Sun: magnetic fields - stars: winds, outflows

1 Introduction

The modelling of the solar wind is usually performed either by basically following the fluid approach (Parker 1963) or by emphasizing the kinetic aspects of the solar wind plasma as they emerge from relatively recent ULYSSES observations at several AU (e.g., Balogh et al. 2001). The kinetic approach has been successful in reproducing, for example, the profile of the electron temperature at large distances (e.g., Zouganelis et al. 2004), but it fails to explain proton behaviour closer to the Sun. In this study we shall largely follow the fluid approach, confining our attention to the dynamics of the ions in the solar wind. In that context, various numerical attempts at modelling the rotating steady solar wind (e.g. Sakurai 1985; Washimi 1990; Tsinganos & Bogovalov 2000; Keppens & Goedbloed 2000) have been successful at reproducing various of its observed features, by using a polytropic equation of state. In another study Usmanov et al. (2000) presented a sophisticated numerical simulation which reproduces many related solar wind observations by using a realistic energy equation which includes heating by Alfvén waves in the WKB approximation, albeit by excluding rotation. This model has been further refined (Usmanov & Goldstein 2003) by incorporating solar rotation and a tilted dipole magnetic field up to 10 solar radii. The full 3D MHD equations are solved in an inner computational domain while in the outer one only the dynamics of the flow along radial fieldlines is solved numerically. Although the study reproduces with some good agreement the latidudinal dependence of the solar wind, there is still a problem in fitting simultaneously the observational data for the density and temperature, both in the inner solar corona and in the extended wind. This is also the case with other MHD simulations of the solar wind which use a polytropic equation of state, unless they include a specific extra heating (Riley et al. 2001).

However, such numerical simulations are rather time consuming for exploring the full range of boundary conditions of the problem and solving at once the wind structure at all distances. Semi-empirical models have also been developed to model coronal holes following Munro & Jackson (1977) (e.g. Guhathakurta et al. 1999; Zangrilli et al. 2002; Zouganelis et al. 2004). Such semi-empirical models usually adopt a global 3D structure of the wind which is to zeroth order consistent with the observations, and then solve for the energy budget along each fieldline. In this way, such models are able to explore more sophisticated energy equations or kinetic effects but cannot solve consistently the structure of the flow itself. Conversely, here we propose to use exact MHD solutions to describe the average dynamical quantities and their latitudinal and radial variations in the wind. As a byproduct, the use of the ULYSSES data for the solar wind is a way to explore the limitations of such analytical models when they are used to model analogous outflows in jets from young stars where the observational constraints are not so well known.

Similarly to the pioneering work of Nerney & Suess (1975a,b,c) (see also Nerney & Suess 1985, and references therein), the approach we use is based on two different sets of meridionally self-similar solutions (see Sect. 2) for the non-polytropic MHD equations presented in earlier connections (Lima et al. 2001, hereafter LPT01; Sauty et al. 1999, hereafter STT99; Sauty et al. 2002, hereafter STT02; Sauty et al. 2004, hereafter STT04). These two approaches allow us to parametrize the latitudinal dependence of the solar wind density, velocity and magnetic fields by solving consistently the momentum equation, including not only solar rotation but also the rotation of the wind itself and the resulting deviations from Parker's spiral.

The first set (following STT99, hereafter called model A) allows only for a rather smooth variation of the density and wind speed with heliolatitude. Although it is an exact solution of the MHD equations, it is derived by keeping only the lowest harmonics in an expansion with the latitude of the various variables of the wind. This model may describe various MHD structures, including a fieldline flaring as in helmet streamers around the equatorial plane or the collimation of the fieldlines towards the axis due to magnetic stresses. It is this last property that enables the model to be suitable for describing jets from young stars by combining plasma ejection from both the central star and the surrounding disk (see STT02, STT04).

The second set (using LPT01, hereafter called model B) assumes that the wind is purely radial in the poloidal plane. Conversely to the previous set, it may allow for a sharper variation of the wind variables in the latitudinal direction. This is well adapted to the radial solar wind of the outer corona and to the observed sharp transition between the slow and the fast solar wind. Note however that it is not equivalent to a pure 1D modelling as it consistently solves the transfield equation in all space and not solely the conservation of energy along the lines.

ULYSSES data for the first polar passages around the period of the last solar minimum offer a complete set of observations of the dynamical variables of the wind with latitude between 1.3 and 2.3 AU. Moreover, solar minimum is a period where the solar wind is exceptionally close to axisymmetry with two main polar coronal holes surrounded by equatorial streamers. In addition, the wind may well be approximated as being rather close to a steady state during that period. For all these reasons, we first constrain the various parameters of the two models from ULYSSES data, as explained in Sect. 3. Then, we use model A (Sect. 4) to predict the inner structure of the solar wind dynamics and geometry consistent with our analysis of the spacecraft data. In the inner region of the corona, the lowest orders in the multipole development with colatitude are more relevant, which justifies the use of model A. We postpone the complete use of model B and its comparison with the present results to a forthcoming paper. We present different ways of constructing global solutions that describe the solar wind from the stellar surface up to the heliopause, and discuss the properties of the solutions thus obtained. In particular, we compare the computed physical quantities close to the solar surface with observations, and check how consistent they are.

   
2 Governing equations for meridional self-similar outflows

2.1 Model A with nonradial streamlines

Model A (see also Tsinganos & Sauty 1992; Sauty & Tsinganos 1994) is an axisymmetric wind model obtained by a self-consistent solution of the full system of the ideal MHD equations. All three components of the velocity and magnetic fields are considered in this particular description and so the model is able to describe the usual flaring of the fieldlines observed in coronal holes, especially during the phase of solar minimum, something which plays an important role in the initial acceleration of the solar wind. In tackling the difficult problem of dealing with both variations, radial distance and latitude, an expansion of the MHD integrals (angular momentum, corotation frequency, etc.) in latitude is made using harmonics. Such an expansion makes the whole system of MHD equations tractable from an analytical point of view and also represents a heliolatitudinal variation of the wind variables which is in accordance with observations. In the limit of a hydrodynamical wind, application of this model to the data analysis of Munro & Jackson (1977) is given in Tsinganos & Sauty (1992). A more detailed parametric study of the properties of this particular solution has been presented in STT99, STT02 and STT04. Conversely to Nerney & Suess (1975a,b,c) where an expansion is made around the equatorial plane we use, in the present model, exact solutions of the whole set of equations restricting the latitudinal variations to their lowest order and using polar values as reference.

In the following we use spherical coordinates $[r,\theta,\varphi]$. As is already known from axisymmetric wind theory, the physically relevant solution passes through various critical points (Weber & Davis, 1967). One is at the Alfvén radius where the poloidal speed $V_{\rm p}$attains the poloidal Alfvén speed such that $B_{\rm p}=\sqrt{4\pi\rho_{\rm p}}V_{\rm p}$. For numerical purposes all quantities have been normalized with respect to values at the Alfvén radius r_* along the polar axis $\theta=0$. At this point the velocity is V_*, the density $\rho_*$, and from the definition of the Alfvén point the magnetic field is $B_*=\sqrt{4\pi\rho_*}V_*$.

The dimensionless spherical radius is denoted by R=r/r_*. The velocity, magnetic, density, and pressure fields of model A are given below. They are functions of the spherical dimensionless radius R and of the co-latitude $\theta$ and can be written in the northern hemisphere as:

$$V_r = V_{*} {f M^2\over R^2} { \cos\theta \over \sqrt{1+\delta f \sin^2\theta} } ~,$$ (1)

$$V_\theta =-V_{*} {M^2\over 2R}{{\rm d}f\over {\rm d}R} { \sin\theta \over \sqrt{1+\delta f \sin^2\theta}} ,$$ (2)

$$V_\varphi = {\lambda}V_{*} \left(\frac{1 - f M^2/R^2}{1- M^2}\right) {R\sin\theta \over \sqrt{1+ \delta f \sin^2\theta} } ,$$ (3)

$$B_r =B_{*} {f\over R^2}\cos\theta ,$$ (4)

$$B_\theta =-B_{*} {1\over 2R}{ {\rm d}f\over {\rm d}R}\sin\theta ,$$ (5)

$$B_\varphi = \lambda B_{*} \left(\frac{1-f/R^2}{1-M^2}\right)R\sin\theta ,$$ (6)

$$\rho=\frac{\rho_*}{M^2}\left(1+\delta f \sin^2\theta\right),$$ (7)

$$P={1\over2}\rho_* V_{*}^2\Pi\left(1+\kappa f\sin^{2}\theta\right) ,$$ (8)

where $\delta$, $\lambda$ and $\kappa$ are parameters while M(R), $\Pi(R)$ and f(R) are three radial functions. The parameter $\delta$ governs the latitudinal variation of density, $\kappa$the latidudinal variation of the pressure and $\lambda$ of the rotation. M(R) is the poloidal Alfvén Mach number, $\Pi(R)$ measures the radial variation of the pressure and f(R) the geometry of the flow.

2.2 Model B with radial streamlines

In the case of model B (Lima & Priest 1993; and LPT01, where the magnetic field is implemented) the main idea is to keep the latitudinal dependences of the various quantities as general as possible. This meant sacrificing one of the components of the velocity and magnetic field (the latitudinal $\theta$ component) which was taken to be zero. The two remaining non-zero components (the radial r and azimuthal $\varphi$ ones) yield stream and fieldlines that, when projected on the meridional plane, are always radial. Although this might constitute a drastic approximation especially in the accelerating part of the solar wind, it is in good agreement with observations after the wind has passed the critical points. The advantage of this approach is that the latitudinal dependences are versatile enough to reproduce sharp variations of the physical quantities, as was indeed observed during the first ULYSSES fly-by pass from the southern to the northern ecliptic hemispheres. A detailed parametric study of this solution can be found in LPT01.

For consistency between the two models and conversely to the original articles, we use below the same normalization (as in model A) at the Alfvén point along the polar axis, where r=r_*, by means of the dimensionless radius R=r/r_*. For this second model B the velocity, magnetic, density, and pressure fields in the northern hemisphere are:

$$V_r= V_{*}{M^2\over R^2} \sqrt{\frac{1+\mu\sin^{2\epsilon}\theta}{1+\delta\sin^{2\epsilon}\theta}},$$ (9)

$$V_\varphi={\lambda}V_{*} \left(\frac{1 - M^2/R^2}{1- M^2}\ri... ...sin^{\epsilon}\theta}{\sqrt{1+\delta\sin^{2\epsilon}\theta}} ,$$ (10)

$$B_r=\frac{B_*}{R^2}\sqrt{1+\mu\sin^{2\epsilon}\theta},$$ (11)

$$B_{\varphi}=\lambda B_* \left(\frac{1-1/R^2}{1-M^2}\right)R\sin^{\epsilon}\theta,$$ (12)

$$\rho=\frac{\rho_*}{M^2}(1+\delta\sin^{2\epsilon}\theta),$$ (13)

$$P={1\over2}\rho_* V_{*}^2\left(\Pi_0+\Pi_1\sin^{2\epsilon}\theta\right) .$$ (14)

Where $\delta$, $\mu$, $\epsilon$ and $\lambda$ are parameters while M(R), $\Pi_0(R)$ and $\Pi_1(R)$ are three radial functions. The parameters $\delta$ and $\lambda$ and M(R) have the same role as in model A. The parameter $\mu$ controls the latitudinal variation of the momentum flux density and $\epsilon$ the sharpness of the density latidudinal variation. $\Pi_0(R)$ and $\Pi_1(R)$ control the radial variation of the pressure, along the polar axis and out of the polar axis, respectively.

2.3 On the assumptions of meridionally self-similar models

Substituting the forms of the physical quantities into the mass and momentum conservation equations we get three ordinary differential equations in radius for the three functions of the radial distance defined in Sects. 2.1 and 2.2. Details of the numerical and analytical techniques can be found in LPT01 and STT02.

The function f(R) simply gives the geometry of the wind. It is a constant if the poloidal fieldlines are radial. In model B, we just need to consider that f=1 throughout the whole wind region, which produces radial poloidal fieldlines. The function f(R) is, in fact, the inverse of the well known expansion factor used in solar wind theory (Kopp & Holzer 1976).

The function M(R) corresponds to the poloidal Alfvén Mach number which is unity at the Alfvén radius, where corotation more or less ceases. In the present approach it is assumed that Alfvénic surfaces of constant M are spherical and do not depend on colatitude. This is a priori a very restrictive assumption but crucial for using self-similarity (Vlahakis & Tsinganos 1998). To verify that it is consistent with ULYSSES observations, we will start by assuming that between 1 and 5 AU M varies like R. This comes from the fact that in this region the poloidal field lines are radial, velocity is constant with distance and density drops as the inverse of the square of the distance (cf. the following section). Then, plotting the variations of M with latitude scaled down at 1 AU, as shown in Fig. 1, we may conclude that this assumption is a good approximation even far from the sun and despite the obvious large scattering in the slow wind region due to the change in magnetic polarity. We get as an average $M(1\;$AU $)=M_0\approx 19$.

   
3 Using ULYSSES data to constrain the parameters

3.1 General procedure

Using the least mean squares method we have fitted the ULYSSES hourly averaged data for the protons at various latitudes with models A and B. In the following sections we will take hourly averaged data from the Swoops Ions experiment onboard ULYSSES, try to fit this data and, from those fits, estimate the corresponding model parameters.

To plot the variation with latitude of the physical quantities to be fitted at the orbit of the spacecraft, we assume that the variation with distance is known such that we can normalize all quantities to their values at 1 AU. For instance, we assume that the velocity at a given latitude is constant with distance and that the poloidal streamlines and fieldlines are radial. From mass and magnetic flux conservation we deduce as well the variation with distance of the mass and particle density, assuming that the radial magnetic field decreases like 1/r². This implies that the toroidal magnetic field varies like 1/r (Tsinganos et al. 2003) such that the magnetic fieldlines trace out Parker's spiral. To deduce the pressure variation with latitude we need an extra assumption on the temperature profile, which we shall explain later.

All quantities evaluated at 1 AU along the polar axis are noted with the subscript 0, except the toroidal field $B_{{\rm T}0}=B_{\varphi}(r=1\;{\rm AU},\theta=\pi/2)$ which is evaluated in the equatorial plane.

	Figure 1: Scaled Alfvén Mach number at 1 AU. The points represent daily averaged data from the Swoops Ions experiment on ULYSSES.
Open with DEXTER

3.2 Particle density

Following our previous argument, the density at 1 AU can be written for model A as

$$n(1\;{\rm AU})=n_0\left(1+\delta f_0\sin^2\theta\right),$$ (15)

while for model B

$$n(1\;{\rm AU})=n_0\left(1+\delta\sin^{2\epsilon}\theta\right) .$$ (16)

n₀ is the particle density on the polar axis and f₀ is the value of f, both at 1 AU. The function f is a constant if the streamlines are radial. To fit these functions to the data on the density we simply scaled them up by plotting (solid curves in Fig. 2),

$$n_{\rm scaled}=r_{({\rm AU})}^2\;n_{\rm measured} .$$ (17)

If we want to go back to the mass density used in the model, we can simply assume a fully ionized gas and use $\rho_0=\mu\;m_{\rm proton}\;n_0$, where $\mu$ is the mean molecular weight.

	Figure 2: Plot of the scaled proton density versus latitude. The points represent hourly averaged data from the Swoops Ions experiment on ULYSSES. The solid curve corresponds to the fit using model A in panel a) and model B in panel b).
Open with DEXTER

Fitting the data points with the curves resulting from model A (Fig. 2a) and model B (Fig. 2b) we arrive at n₀, $\delta f_0$ and, for the case of model B, $\epsilon$ (see Table 1).

3.3 Momentum flux density

With the streamlines assumed radial, the second step is to fit the momentum flux density

$$\Phi = n V_r^2 r^2 .$$ (18)

We choose this function instead of the radial velocity because in the models addressed here this quantity is independent of both the distance from the sun and the value of $\delta$ previously determined. For model A, at 1 AU we have

$$\Phi(1\;{\rm AU}) = n_0 V_0^2 \cos^2\theta ,$$ (19)

while for model B

$$\Phi(1\;{\rm AU}) = n_0 V_0^2 \left(1+\mu\sin^{2\epsilon}\theta\right) .$$ (20)

Thus, from the two fits (Fig. 3) we deduce the value of V₀. In addition, for model B, we get the value of $\mu$, as displayed in Table 1, assuming the value of n₀ deduced from the fit of the density data. In Table 1 (in parenthesis with *) we also show the two values of n₀ and V₀ for model A deduced by fitting the data on the momentum flux density with two parameters (n₀ and V₀) rather than one (V₀). Although the discrepancy between the two values of V₀ obtained by fitting the momentum flux density as given by models A and B, using a single parameter, is rather high, note that the values of n₀ and V₀ obtained by fitting both parameters using model A are close to the ones obtained with model B. In the following, we have chosen to restrict ourselves to the first set of parameters obtained by fitting both expressions with one single free parameter.

3.4 Radial magnetic field

A similar procedure can be used to deduce parameters referring to the magnetic field. In particular, we plot in Fig. 4

$$B_{r,\rm scaled}=r_{(\rm AU)}^2 B_{r,\rm measured}$$ (21)

and compare it to

$$B_r(1\;{\rm AU})=B_{0}\cos\theta ,$$ (22)

for model A and

 

$$B_r(1\;{\rm AU})=B_{0}\sqrt{1+\mu\sin^{2\epsilon}\theta} \quad {\rm if~ }\theta<\pi/2{\rm ~(north)},$$ (23)

$$B_r(1\;{\rm AU})=-B_{0}\sqrt{1+\mu\sin^{2\epsilon}\theta} \quad {\rm if~ }\theta>\pi/2{\rm ~(south)},$$

for model B. The fit gives us the value of the magnetic field on the polar axis at 1 AU, B₀ in Table 1. Because of the relative flip of the magnetic axis and the rotation axis, northern and southern polar coronal holes tend to mix at low latitudes around the equatorial plane. Besides, the system is not exactly axisymmetric and part of the fast solar wind coming from the northern polar coronal hole can reach southern latitudes and vice versa. To avoid bias to the mixing at low latitudes, we can redo the same procedure excluding data points in the range of latitudes between -20 and 20 degrees. This yields the values in Table 1 indicated by **. As far as the radial magnetic field is concerned this does not change drastically the values obtained.

Table 1: Parameters obtained from models A and B at 1 AU deduced from the fitting of ULYSSES hourly averaged data. All quantities evaluated at 1 AU along the polar axis (except the toroidal field evaluated on the equatorial plane) are noted with the subscript 0. (*) correspond to those values determined by fitting the data on the momentum flux density using least-mean squares and determining the two parameters simultaneously. (**) correspond to using the data on the magnetic field without the points between $-20\hbox {$^\circ $ }$ and  $20\hbox {$^\circ $ }$ of latitude. Values in parentheses were not taken into consideration in further calculations. Note that for model B f₀=1 and $\kappa$ has been defined using Eq. (29).

3.5 Toroidal magnetic field

The toroidal or azimuthal component of the magnetic field $B_\varphi$, denoted by $B_{\rm T}$, gives us a way to determine the value of the remaining free parameter $\lambda$. We proceed similarly to what we did for the radial magnetic field, plotting in Fig. 5

$$B_{\varphi , \rm scaled}= r_{(\rm AU)} B_{\varphi,{\rm measured}} ~$$ (24)

and comparing it to

 

$$B_\varphi(1\;{\rm AU})=-B_{{\rm T}0}\sin\theta~, \quad {\rm if~ }\theta<\pi/2{\rm ~(north)},$$ (25)

$$B_\varphi(1\;{\rm AU})=B_{{\rm T}0}\sin\theta~, \quad {\rm if~ }\theta>\pi/2{\rm ~(south)},$$

for model A and

 

$$B_\varphi(1\;{\rm AU})=-B_{{\rm T}0}\sin^{\epsilon}\theta , \quad {\rm if~ }\theta<\pi/2{\rm ~(north)},$$ (26)

$$B_\varphi(1\;{\rm AU})=B_{{\rm T}0}\sin^{\epsilon}\theta , \quad {\rm if~ }\theta>\pi/2{\rm ~(south)},$$

for model B.

The fit gives us the value of the toroidal magnetic field in the equatorial plane at 1 AU, $B_{{\rm T}0}$, to which we refer in Table 1. Because of the relative flip of the magnetic axis and the rotation axis, northern and southern coronal holes of opposite polarities tend to mix, as we have already pointed out in the previous section. As it is almost impossible to disentangle them, we fit after excluding the data points in the range of latitudes between -20 and 20 degrees (solid lines in Fig. 5). This yields the parameter $B_{{\rm T}0}$ in Table 1 indicated by **. Conversely to the radial magnetic field, the reversal around the equator is a rather drastic effect and tends to yield an underestimate of the real value of  $B_{{\rm T}0}$. For this reason we prefer to use those values to determine the parameter $\lambda$. Unfortunately, $\lambda$depends on the values of the magnetic field, M₀ and R₀, the last one being

$$R_0=\frac{r(1\;{\rm AU})}{r_*} \cdot$$ (27)

So, $\lambda$ cannot be determined without a knowledge of the location of the Alfvén radius in real space.

In principle, we could constrain the value of $\lambda$ by using both the toroidal magnetic and velocity fields. Unfortunately, as shown in Fig. 6, the data on the toroidal velocity cannot be used directly. Firstly, we would expect the velocity to be of the same sign on both hemispheres as on the surface of the Sun, being zero on the two polar axes. Secondly, the rather small values of the measured velocities may be considered to be lower than the capability of the instrument to measure them. Strangely, what Fig. 6 really shows is a kind of reversal of the toroidal velocity close to the equator, similarly to what happens with B_rand $B_\varphi$. In fact, in these models if the poloidal velocity is constant after the Alfvén point and the lines radial, then V_r=V_* and fM² /R²=1 (Eqs. (1) and (9); of course in model B, f=1) then the toroidal velocity is exactly 0 from Eqs. (3) and (10). Though we did not consider this assumption before starting the integration, the terminal toroidal velocity ended up being rather small compared to the radial one and the terminal radial velocity close to the radial velocity at the Alfvén point.

	Figure 3: Plot of the scaled momentum flux density versus latitude. As in Fig. 2, the points represent hourly averaged data from the Swoops Ions experiment, the solid curve corresponds to the fit using model A in panel a) and model B in panel b).
Open with DEXTER

	Figure 4: Plot of the scaled radial magnetic field versus latitude. The points represent hourly averaged data from the Swoops Ions experiment, the solid curve corresponds to the fit using model A in panel a) and model B in panel  b).
Open with DEXTER

Figure 5: Plot of the scaled toroidal magnetic field versus latitude. The points represent hourly averaged data from the Swoops Ions experiment. The curves correspond to the fit using model A in panel a) and model B in panel b). In the fit with the dashed line all points have been taken into account while in the fit with the solid line the points between -20 and 20 degrees latitude (i.e. between the two dotted vertical lines) have been excluded. Note that in panel a) the solid and dashed lines almost coincide.

Open with DEXTER

	Figure 6: Plot of the toroidal velocity. The lack of precision in the data prevents them from being used to constrain the parameters. As explained in the main text, this corresponds however to rather small toroidal velocities consistent with our calculations.
Open with DEXTER

3.6 Temperature

The last information from the Swoops Ion Experiment to constrain the two models is the proton temperature profile. With the single fluid approximation used, the average temperature corresponds mainly to that of the protons. The temperature can be used to determine the last parameter $\kappa$ of model A. However, as the energy equation is not properly defined, we may obtain at the end of our calculation an effective temperature profile resulting from considering all sorts of pressures, like wave pressure. This effective temperature might be much larger than the kinetic one. The adopted procedure is used to constrain the latitudinal variation of the pressure from the temperature profile but we cannot control the temperature variation with distance. All we can do is to check a posteriori how consistent it is with the measured profile and predict the distribution of the heating in the corona needed to support the wind.

As the temperature varies, both with latitude and radial distance, we have assumed that the variation with latitude should be negligible inside the fast solar wind part. Thus, taking out latitudes between $-20^\circ$ and $20^\circ$, we fit the proton temperature at high latitudes with a power law function for the polar temperature of the form

$$T(R,\theta=0)=T_{0}\left(\frac{R_0}{R}\right)^\tau \cdot$$ (28)

This gives $\tau \approx 1.047$ (Fig. 7). This is somewhat similar to the fit for the electron temperature given in Issautier et al. (2001). As another test, we also tried to use only data points beyond 1.52 AU instead of 1.3 AU which yielded $\tau\approx 1.02$ which is very similar to the above value.

We then fitted a more general form of the temperature profile to the whole set of points. The latitudinal dependence was chosen to be consistent with Eqs. (7) and (8), assuming that the poloidal streamlines are radial such that $f(r\ge1\;{\rm AU})=f_0$,

$$T(R,\theta)=T_{0}\left(\frac{R_0}{R}\right)^\tau \frac{1+\kappa f_0 \sin^2\theta}{1+\delta f_0 \sin^2\theta} \cdot$$ (29)

As $\delta f_0$ is already known from the fit of the density, we deduce from the present fit (Fig. 8) the values of $\kappa f_0$ shown in Table 1.

	Figure 7: Plot of the scaled proton temperature versus distance, using only data for the fast wind for colatitudes $\varphi$ above $70^{\circ }$ and below ${-}70^{\circ }$. We fit the data with a simple power law $\sim$$1/R^\tau$ which gives $\tau \approx 1.047$.
Open with DEXTER

Of course we cannot use directly the equations of model B to find the latitudinal dependence of the pressure, as this is known only a posteriori once the full solution is constructed. However, to constrain model A using parameters of set II we can use the same latitudinal dependence assuming of course $f(r\ge1\;{\rm AU})=f_0=1$. Then we get the value of $\kappa$ using Eq. (29), consistent with the value of $\delta$ deduced previously from model B.

	Figure 8: Plot of the scaled proton temperature versus latitude. As in Fig. 2, the points represent hourly averaged data from the Swoops Ions experiment, the solid curve corresponds to the fit using model A in panel a) and model B in panel b).
Open with DEXTER

To fit the temperature data we have been using the profile of the so called small temperature from the experiment, which is likely to be an underestimate of the real temperature. However, even using data on the high temperature leads to equivalent values for $\kappa f_0$.

3.7 Parameters at 1 AU

The parameters obtained are summarized in Table 1. From the previous sections note that the more general form of the latitudinal dependences of the various quantities shown in Eqs. (9)-(14) (model B) yields a more convincing fit to the data than the form of Eqs. (1)-(6) (model A). Model A tends to overestimate the polar speed V₀ and underestimate the polar density n₀. Note, however, that the values in parenthesis for SET I in Table 1, which were determined using other statistical methods, are very close to those of set II. In any case, considering the precision to which we worked out the parameters, the discrepancy between the two methods remains at a rather satisfactory level in the sense that the two models did not give drastically different outputs.

   
4 Construction of a complete solution

4.1 Initial conditions at the Alfvén surface

So far, the appropriate fitting of the ULYSSES data has made it possible to deduce all parameters of models A and B. Thus, we have analytical expressions for the heliolatitudinal dependence of the variables of the wind at 1 AU. The next step is to proceed towards the construction of a complete solution for all heliocentric radial distances. This may be accomplished by integrating over radial distance the ODEs of each of the two models. In order to do so, we need to estimate the parameters and physical quantities such as the polar velocity and density at the Alfvén radius. This cannot be done without further assumptions. A posteriori we will verify that our estimations are in agreement with the general numerical output. If not we will iterate and change consistently all our estimates until we find an agreement between the input parametric values and their output.

In the following we have chosen to use the two sets of parameters in Table 1 and integrate the equations of model A. This choice is related to the fact that model A is able to reproduce a more general geometry of the fieldlines from the extreme case of flaring towards the equator, to collimation towards the pole, passing by a configuration in which those lines are purely radial. Model B is more limited in this sense since it corresponds to this later geometry. However, we plan in future work to integrate the equations of model B with the same parameters and compare the obtained results with the present ones. These two sets of parameters will thus yield two sets of solutions: set I and set II, corresponding to the parameters previously obtained by models A and B, respectively. On one hand, model B includes higher order latitudinal variations than model A. Thus, as we have already mentioned, the former fits more accurately the rapid latitudinal variations of the observed data. On the other hand we can assume that the dynamics is mainly affected by the lowest order variations with latitude. Assuming this, we use set II and integrate the equations of model A, ignoring the values of $\mu$ and $\epsilon$.

If we know the physical quantities at 1 AU, on the polar axis, $\rho_0$, V₀ and B₀, the value of the Alfvén Mach number at that same distance follows from

$$M_0^2=\frac{4\pi\rho_0V_0^2}{B_{0}^2} \cdot$$ (30)

Substituting Eq. (30) into Eq. (7) we get the value of the mass density at the Alfvén surface,

$$\rho_*=\frac{4\pi\rho_0^2V_0^2}{B_{0}^2} \cdot$$ (31)

The location of the Alfvén surface r_* depends on the degree of streamline flaring towards the equator or collimation towards the polar axis and on the plasma acceleration beyond the Alfvén surface.

The expansion factor is given by f. At 1 AU, bearing in mind that f=1 at R=1, the relative variation in the streamline expansion is given by $f(1\;{\rm AU})=f_0$. From Table 1, we see that guessing an initial value for f₀ provides the initial values for $\delta$ and $\kappa$.

The total acceleration gained between the Alfvén surface and 1 AU can be parametrized by the ratio

$$\eta=\frac{V_0}{V_*} \cdot$$ (32)

By guessing an initial value for $\eta$ we immediately deduce some initial guess for V_* and $B_*=\sqrt{4\pi\rho_*}V_*$. Then, the value of r_* follows from,

$$R_0=\frac{1\;{\rm AU}}{r_{*({\rm AU})}} =\frac{215\;r_\odot}{r_{*({r_\odot})}}= \sqrt{f_0\frac{B_*}{B_{0}}} ,$$ (33)

i.e., from R₀ we get r_* in AU or solar radii.

The next step is to get the value of $\lambda$ from r_*, f₀, B₀ and  $B_{{\rm T}0}$ by combining Eqs. (4) and (6) or Eqs. (11) and (12). Assuming that the lines are almost radial and that we are at large distances from the Alfvén radius ($R_0\gg 1$), we get

$$B_{{\rm T}0}\approx\lambda B_{*}\frac{R_0}{M^2_0} \Leftrighta... ...mbda \approx \frac{B_{{\rm T}0}}{B_{*}}\frac{M^2_0}{R_0} \cdot$$ (34)

Finally, $\nu$ provides also the value of the gravitational potential. Since the mass of the Sun ${\cal M}$ and the gravitation constant ${\cal G}$ are known, we can write

$$\nu =\sqrt{{2\cal G M} \over r_* V^2_*}\cdot$$ (35)

Altogether, we can summarize our procedure as follows:

• First, determine M₀ and $\rho_*$ from the data at 1 AU.
• Second, guess f₀ and get as an output the parameters $\delta$ and $\kappa$.
• Third, guess $\eta=V_0/V_*$, which then gives V_* and B_*.
• Finally, by combining all of the above, r_*, $\lambda$, and $\nu$ will follow.

4.2 Iterating

The first set of guessed parameters allows us to compute the whole solution using our set of equations and a Runge-Kutta integration algorithm. At this point there is, of course, no reason why the final output should correspond to the initial input. In particular, the computed f₀, $\eta$ and r_* may not necessarily coincide with the initial values guessed for these parameters.

Then we either use these new values of f₀ and $\eta$ (or r_*) from our computation or some other estimates of those parameters to recalculate consistently all the other parameters. Then, after trial and error, we hope to get all parameters (on the input and output phases of our calculation) to converge to each other.

A priori there is no reason why f₀, $\eta$ and r_* should converge simultaneously in this process. We were lucky enough that such convergence was attained almost completely in the case of the second set of parameters (set II, fom model B). However, for the first set (set I, from model A), as $\eta$ is completely arbitrary we did not required that it remains similar on input and output. Instead, we have rescaled V_* such that everything remains consistent with observations at 1 AU.

A more serious problem was that we were not able to produce physical solutions with positive pressure everywhere considering the high values of the magnetic field at 1 AU. This may be due to the fact that the simple latitudinal dependence of model A does not allow consistent treatment of the role played by the equatorial sheet in the total force balance. Nevertheless, we were able to produce reasonable solutions with consistent magnetic field values at the surface of the sun by dividing the obtained values of B₀ by a factor of 7 in the case of set I and by 10 for set II. We note that Zangrilli et al. (2002) had similar problems in modelling physical quantities and their variation with latitude and distance both in the low and in the extended corona. A similar problem of negative pressure values at low latitudes has been encountered in Tsinganos & Trussoni (1991).

4.3 Results

In Table 2 we present the parameters of the final iteration with which it is possible to solve the consistency problem mentioned above.

Table 2: Parameter sets I and II integrating with model A, on last iteration.

Figure 9: Plot of the poloidal velocity as a function of distance for the parameters of Table 2. In panel a) we show the solution using parameters of set I and in panel b) we show the solution using parameters of set II. In each graph, the upper curve corresponds to $\alpha =0$ (at the pole) and $\alpha$ increases by 0.1 between consecutive curves towards the lower curve ( $\alpha =0.5$ which corresponds to the line crossing r* at $\theta =45^\circ$).

Open with DEXTER

	Figure 10: Plot of the effective poloidal temperature as a function of distance. As in Fig. 9, panel a) corresponds to parameter set I and panel b) to set II and $\alpha$ ranges from $\alpha =0$ to $\alpha =0.5$.
Open with DEXTER

	Figure 11: Plot of the fieldlines on a) for the solution using parameters of set I and on b) for the solution using parameters of set II. Distances are given in units of r* which is given in Table 2.
Open with DEXTER

	Figure 12: The same as Fig. 11 but zooming in on the central region, in a) for the solution using parameters of set I and in b) for the solution using parameters of set II. Distances are given in units of r*.
Open with DEXTER

Clearly, set II is better, both from the viewpoint of consistency between the parameters fed at the last iteration and obtained after computation and also because it gives more reasonable values of the velocity and effective temperature (compare Figs. 9 and 10). This is related to the fact we have already mentioned that parameters of model A give overestimates of the polar velocity and underestimates of the polar density.

The values for the effective temperature for set I are high with a maximum of around $35\times 10^6$ K on the polar axis. Even considering this temperature as an effective one this value may be regarded as unreasonably high. In fact, the plotted effective temperature is defined as

$$T=\frac{P_{\rm total}}{n_{\rm proton}k_{\rm B}} ,$$ (36)

where the total pressure may include various terms besides the usual kinetic pressure,

$$P_{\rm total}= P_{\rm kinetic} + P_{\rm waves} + P_{\rm radiative} + ... .$$ (37)

The values of the effective temperature for set II are still too high but they are in the usual range of most MHD models of the solar wind. For instance, in the kinetic simulations by Landi & Pantellini (2003), they obtain temperatures around 3 times higher than the values usually adopted for the solar corona. Similar problems appear in Zouganelis et al. (2004) for the electron temperature, though they seem to have a better treatment of the proton temperature. The high values obtained in the present study might be related to the fact that a full energy equation is not included in the models. Yet, with set II, the maximum effective temperature of $15\times 10^6$ K can be related to an important deposit of torsional Alfvén waves.

Let us consider that the kinetic temperature of the protons may be of the order of $5\times 10^6$ K at the maximum temperature. This may even be an underestimate (e.g. see the values given by Zangrilli et al. 2002). The maximum temperature in the plot (set II) is below $15\times 10^6$ K, corresponding to distances between 1.1 and $1.3~ r_\odot$. In this range the density decreases from around $n=6\times 10^5$ cm^-3 to $n=2.5\times 10^5$ cm^-3. The extra $10\times 10^6$ K, if assumed to be a result of torsional Alfvén waves, would then correspond to a wave pressure that is at most

$$P_{\rm Alfven}=\frac{\delta B_\varphi^2}{8\pi} = n_{\rm proton}k_{\rm B} T_{\max} \leq 8.3\times 10^{-4}~ {\rm dynes} ,$$ (38)

yielding a reasonable fluctuation in the magnetic field of the order

$$\frac{\delta B_\varphi}{B_r}\approx \frac{0.144}{0.427} \approx 0.34 .$$ (39)

By comparison, Usmanov et al. (2000) use a maximum of relative Alfvén wave amplitudes of the order $\delta B /B \approx 0.6$, although using a WKB approximation. So we may conclude that we get sensible results. In fact, we could use our constraints on the dynamic of the wind to predict the distribution of wave damping and heating. This is something we intend to do in the future in a more precise way.

We plot in Fig. 9 the poloidal velocity for both sets of parameters and for $\alpha=f\sin^2\theta$ ranging between 0 and 0.5, i.e., between the pole and somewhere at mid-latitude as the line $\alpha =0.5$ is at an angle of $45^\circ$ at r=r_*. In both cases the asymptotic speed has been attained by $4\;R_{\odot}$ along the polar axis, which we believe to be a reasonable value. It may sound as if the acceleration is rather fast compared to the usual measurements close to the Sun. However, we should bear in mind that this constitutes an upper limit on the velocity and acceleration, since at mid-latitudes the radial speed is a fraction of the polar speed as shown in Fig. 9. Observed velocity profiles in coronal holes are always integrated over the line of sight, mixing various latitude inside the coronal hole. This means that to make a detailed comparison with observations of our velocity profile in the low corona we need a more sophisticated way to plot the average speed inside the coronal hole by integrating over latitudes along the line of sight.

Figure 10 shows the effective poloidal temperature for both sets of parameters and for $\alpha$ ranging between 0 and 0.5. The values of the effective poloidal temperature, especially at mid-latitudes, are in quite good agreement with observations for set II, although we have already discussed that it is rather high close to the polar axis. Moreover, we get very good agreement of the temperature profile at large distances where ULYSSES measurements are made. This is a nice surprise as we cannot constrain the temperature profile from the model. It is a byproduct of our calculation.

However, another limitation of our modelling is that, although the initial values of the magnetic field at the solar surface in the coronal holes are quite reasonable, we have failed to reproduce these values observed by ULYSSES at 1 AU by a factor of 7 in set I to 10 in set II as seen from the values we have used as input in Table 2. Higher values of the magnetic field at large distances give too high a value of the magnetic field at the surface of the sun ($\approx$100 G) and negative pressure, which is related to the dead zones in the equatorial plane being too large. By taking only the lower dipolar order in the latitudinal dependence we do not reproduce correctly the effect of the current sheet on the equatorial plane at large distances. Usmanov et al. (2000) showed that the current sheet is a rather important ingredient of the force balance at large distances. We may solve this problem by using model B instead of A to model the solar wind in the outer corona.

On the opposite, we should stress that our solutions, despite being constrained from measurements at 1 AU, reproduce fairly well the dipolar structure of the magnetic field close to the sun up to 2 solar radii as shown in Fig. 12 with $f\sim 1/R$, $B_r\sim \cos\theta/R^3$, $B_\theta \sim \sin\theta/2R^3$. Moreover, one main result of these calculations is that we are not getting any collimation at large distances so far. This is in contradiction with our initial expectation (see Heyvaerts & Norman 1989; Nerney & Suess 1985; Usmanov et al. 2000) as the wind is carrying some current which closes into the equatorial current sheet. We have shown in STT02 that collimation is not a necessary condition as long as the force from the toroidal magnetic field can be balanced by either the curvature of the poloidal velocity field or by the pressure gradient. We also note that the efficiency of the magnetic rotator as defined in STT99 is rather low with a value of the order $\varepsilon/(2\lambda^2)\approx-50$. We refer the reader to STT99 Eq. (3.12a) for the definition, as it would be lengthy and out of order to reproduce it here. Thus, the absence of collimation was not really surprising considering the parametric analysis of STT02.

   
5 Conclusions

We have proposed here a method for modelling the rotating and magnetized solar wind constraining the parameters from ULYSSES data. We have used two models: model A from STT99 and STT02 and model B from LPT01. These axisymmetric solutions described in Sect. 2 include the effects of both the fast and the slow components on the global structure and the dynamics of the wind. Yet, they aim at describing the wind mainly close to its polar axis, especially model A. They also consistently include rotation, conversely to the usual axisymmetric models for the solar wind.

Model A tends to over- or underestimate some physical quantities, especially along the polar axis. Model B allows a better fit because it can account for faster variations with latitude. In Sect. 3 we have described in detail how to analyse the data in the frame of our modelling.

Then we have used model A and the two sets of parameters obtained through the latitudinal dependences appearing in the solution for models A and B (sets I and II, respectively) to construct the full solution from the solar surface up to ULYSSES orbit as described in Sect. 4. Our results show that there are still difficulties in trying both to fit the latitudinal profile of the solar wind using ULYSSES data and to reconstruct the solution all the way back from the spacecraft to the solar wind base. These inconsistencies might be connected to intrinsic limitations of the models used. They might also be associated with the contamination of pure fast wind by pure slow wind flows, which affects the variation with latitude (Usmanov & Goldstein 2003).

On the whole, we have been able to reproduce the density and velocity profiles from the solar surface to ULYSSES orbit reasonably well. We have also obtained a good value for the rotation frequency of the Sun, which was not in the initial input. The acceleration along one streamline may appear quite strong when compared to the usual velocity profile measurement. Yet one should bear in mind that the observed profiles in coronal holes are integrated over different latitudes. With variations with latitude being quite rapid, the averaged velocity along the line of sight increases much more slowly than along the polar axis.

One characteristic of this non-polytropic model is that flux tubes and coronal holes with large expansion factors tend to give lower asymptotic speeds than those with small expansion factors. This is in contradiction with the usual polytropic models (e.g. Kopp & Holzer 1976) but in good agreement with observations as shown by Wang (1995) and Wang & Sheeley (2003).

We also reproduce fairly well the temperature profile at large distances. In the low corona, the agreement is not too bad providing that this is interpreted as an effective temperature, especially around its maximum between 1 and $2\; R_\odot$. There, the total effective temperature can be divided roughly into 1/3 of kinetic temperature and 2/3 of temperature coming from wave pressure. If those waves are torsional Alfvén waves we find that the toroidal magnetic field fluctuations need to be more or less of the order of the radial magnetic field magnitude itself. Of course a more complete study of heating and wave damping is needed to support this conclusion. A more quantitative study is in order but this has been postponed to the future because we need a realistic energy equation to calculate it properly.

We could reproduce the dipolar structure of the magnetic field in the first few solar radii with approximately 0.5 G at the surface of the coronal hole. However, we were not able to model the correct magnitude of the field up to 1 AU but only the shape of Parker's spiral. This also needs to be improved. The ultimate goal would be to combine both this model (which uses STT99) at small distances and the LPT model at large distances where the lines are radial, to make a more consistent description of the whole wind. Assuming radial fieldlines at large distances may solve the magnetic field problem by reducing the way it decreases with distance and allowing the second model to give better fits to the observed magnetic fields.

Acknowledgements

The authors thank the anonymous referee for his/her valuable comments which helped to improve the presentation of this paper. They wish to acknowledge the SWOOPS instrument team and the ESA and NSSDC COHO Web archives for ULYSSES data. J. J. G. Lima wish to acknowledge the financial support of the Université Paris VII and the hospitality of everyone at the Observatoire de Paris, Meudon. This work was supported by grant POCTI/1999/FIS/34549 approved by FCT and POCTI, with funds from the European Community programme FEDER and by the EEC RT Network HPRN-CT-2000-00153.

References

• Balogh, A., Marsden, R. G., & Smith, E. J. 2001, The Heliosphere near solar minimum (Springer Praxis Books) (In the text)
• Guhathakurta, M., Fludra, A., Gibson, S. E., Biesecker, D., & Fisher, R. 1999, JGR, 104, 9801 [NASA ADS] (In the text)
• Heyvaerts, J., & Norman, N. 1989, ApJ, 347, 1055 [NASA ADS] [CrossRef] (In the text)
• Issautier, K., Meyer-Vernet, N., Pierrard, V., & Lemaire, J. 2001, Ap&SS, 277, 189 [NASA ADS]
• Keppens, R., & Goedbloed, J. P. 2000, ApJ, 530, 1036 [NASA ADS] [CrossRef] (In the text)
• Kopp, R. A., & Holzer, T. E. 1976, Sol. Phys., 49, 43 [NASA ADS] (In the text)
• Landi, S., & Pantellini, F. 2003, A&A, 400, 769 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Lima, J. J. G., & Priest, E. R. 1993, A&A, 268, 641 [NASA ADS] (In the text)
• Lima, J. J. G., Priest, E. R., & Tsinganos, K. 2001, A&A, 371, 240 [EDP Sciences] [NASA ADS] [CrossRef] (LPT01) (In the text)
• Munro, R. H., & Jackson, B. V. 1977, ApJ, 213, 874 [NASA ADS] [CrossRef] (In the text)
• Nerney, S. F., & Suess, S. T. 1975a, ApJ, 196, 837 [NASA ADS] [CrossRef] (In the text)
• Nerney, S. F., & Suess, S. T. 1975b, ApJ, 200, 503 [NASA ADS] [CrossRef]
• Nerney, S. F., & Suess, S. T. 1975c, Sol. Phys., 45, 255 [NASA ADS] [CrossRef]
• Nerney, S. F., & Suess, S. T. 1985, ApJ, 296, 259 [NASA ADS] [CrossRef] (In the text)
• Parker, E. N. 1963, Interplanetary Dynamical Processes (New York: Wiley Intersience)
• Riley, P., Linker, J. A., & Mikic, Z. 2001, JGR, 106 (A8), 15, 889 (In the text)
• Sakurai, T. 1985, A&A, 152, 121 [NASA ADS] (In the text)
• Sauty, C., & Tsinganos, K. 1994, A&A, 287, 893 [NASA ADS] (In the text)
• Sauty, C., Tsinganos, K., & Trussoni, E. 1999, A&A, 348, 327 [NASA ADS] (STT99) (In the text)
• Sauty, C., Trussoni, E., & Tsinganos, K. 2002, A&A, 348, 327 [NASA ADS] (STT02) (In the text)
• Sauty, C., Trussoni, E., & Tsinganos, K. 2004, A&A, 421, 797 [EDP Sciences] [NASA ADS] [CrossRef] (STT04) (In the text)
• Tsinganos, K., & Bogovalov, S. 2000, A&A, 356, 989 [NASA ADS] (In the text)
• Tsinganos, K., & Sauty, C. 1992, A&A, 255, 405 [NASA ADS] (In the text)
• Tsinganos, K., & Trussoni, E. 1991, A&A, 249, 156 [NASA ADS] (In the text)
• Tsinganos, K., Vlahakis, N., Bogovalov, S., et al. 2003, Ap&SS, 287, 103 [NASA ADS] (In the text)
• Usmanov, A. V., Goldstein, M. L., Besser, B. P., & Fritzer, J. M. 2000, JGR 105 (A6), 12675 (In the text)
• Usmanov, A. V., & Goldstein, M. L. 2003, JGR, 108 (A9), 1354 (In the text)
• Vlahakis, N., & Tsinganos, K. 1998, MNRAS, 298, 777 [NASA ADS] [CrossRef] (In the text)
• Weber, E. J., & Davis, L. J. 1967, ApJ, 148, 217 [NASA ADS] [CrossRef]
• Wang, Y.-M. 1995, ApJ, 449, L157 [NASA ADS] (In the text)
• Wang, Y.-M., & Sheeley, N. R. 2003, ApJ, 587, 818 [NASA ADS] [CrossRef] (In the text)
• Washimi, H. 1990, Geophys. Lett., 17, 33 [NASA ADS] (In the text)
• Zangrilli, L., Poletto, G., Nicolosi, P., Noci, G., & Romoli, M. 2002, ApJ, 574, 477 [NASA ADS] [CrossRef] (In the text)
• Zouganelis, I., Maksimovic, M., Meyer-Vernet, N., Lamy, H., & Issautier, K. 2004, ApJ, 606, 542 [NASA ADS] [CrossRef] (In the text)