A&A 400, 769-778 (2003)
DOI: 10.1051/0004-6361:20021822
S. Landi^{1} - F. Pantellini^{2}
1 - Dipartimento di Astronomia e Scienza dello Spazio,
Largo Enrico Fermi 2, 50125 Firenze, Italy
2 -
Observatoire de Paris, LESIA,
5 Place Jules Janssen, 92195 Meudon, France
Received 30 August 2001 / Accepted 9 December 2002
Abstract
We present and discuss a completely self-consistent kinetic simulation of
a steady state transonic solar type wind. The equations of
motion of an equal number of protons and electrons plunged
in a central gravitational field and a self-consistent
electric field are integrated numerically.
Particles are allowed to make binary collisions with
a Coulombian scattering cross-section. The velocity
distributions of the particles injected
at the boundaries of the simulation domain
are taken to be Maxwellian.
As anticipated by previous authors we
find that the transonic solution implies the existence
of a peak in the proton equivalent potential at some
distance above the sonic critical point. Collisions appear to be
the fundamental ingredient in the process of accelerating the
wind to supersonic velocities. For a given temperature at the
base of the simulation domain the acceleration efficiency
decreases with decreasing density. The reason is that
the plasma has to be sufficiently collisional for the
heat flux to be converted efficiently into plasma bulk velocity.
Concerning the heat flux we find that
even when in the vicinity of the sonic point
the collisional mean free path of a thermal
particle is significantly smaller than
the typical scales of variation of
the density or the temperature, the electron heat flux cannot be described
conveniently by the classical Spitzer-Härm conduction law;
not even in most of the subsonic region. Indeed, in the simulations where
a transonic wind forms the heat flux has been found to
strongly exceed the Spitzer-Härm flux, in opposition
to recently published results from multi-moment models.
We emphasize that given the high coronal temperatures we use in our simulations
(3 times the typical solar values) we do not expect the results
presented in this report to be uncritically transposable to the case of
the "real'' solar wind. In particular, the quantitative aspects of our results
must be handled with some care.
Key words: Sun: solar wind - stars: winds, outflows - plasmas - conduction - methods: numerical
At all heights, from the bottom of the corona up into the interplanetary space, the solar atmosphere is a permanently expanding, out of thermodynamic equilibrium and fully ionized plasma. During the 1950's the recognition of the weak collisionality of the solar wind conveyed some doubts concerning the ability of fluid models to describe the solar atmosphere conveniently. As a consequence Chamberlain (1960), largely influenced by the theories on gas evaporation from planetary atmospheres, published the first kinetic model of the solar wind. In Chamberlain's evaporation model the wind is subsonic at the Earth's orbit in clear opposition with the supersonic solution of the fluid equations proposed a few years earlier by Parker (1958). During the early 1960's in situ measurements confirmed the supersonic nature of the solar wind and kinetic models just fell in disuse for some time. In the early 1970's it became clear that Chamberlain's erroneous prediction of a subsonic wind was the consequence of having mistakenly assumed that the charge neutralizing electric field was the Pannekoek-Rosseland field (e.g., Rosseland 1924). The latter is based on the assumption of a static solar atmosphere which has been a privileged working hypothesis since Laplace's Traité de mécanique céleste, published in the early years of the nineteenth century, but has been shown to be completely at odds with observations. After the definitive relaxation of the static approximation for the electric field, kinetic models of the solar wind were back on stage again (Jockers 1970).
The simplest, and most widely used kinetic models are the so called exospheric models (e.g., Lemaire & Sherer 1971). In these models the solar atmosphere is assumed to change from fully collisional to collisionless at a sharply defined level called the exobase. Above the exobase, conventional exospheric models assume a monotonically decreasing equivalent proton potential (cf. Eq. (4)) and no protons coming into the system from infinity, so that, by construction, all protons have anti sunward directed velocities. For non pathological distribution functions this means that the plasma bulk velocity at the exobase is of the order of the proton thermal velocity, i.e. of the order of the sound velocity. For example, is the mean velocity of a proton population with Maxwellian velocity distribution truncated for (subscripts and refer to the radial direction with respect to the center of the Sun). Since the typical velocity of a proton at the exobase is by construction of the order of the radial bulk velocity, it follows that the exobase is located at a heliocentric distance comparable to the distance r_{*} where the subsonic-supersonic transition is located (the sonic critical point). However, as demonstrated graphically by Jockers (1970), the proton potential cannot be monotonic from deep inside the corona, where the bulk velocity is supposed to be small compared to sound speed and where static approximation may apply, out to infinity, where the wind is supersonic. Jockers anticipated that on its way from the corona to infinity a transonic wind must overcome a maximum in the proton potential such that , where is the location of the maximum. In two recent papers Scudder (Scudder 1996a,b) pushes a step farther by identifying the critical point of Parker's fluid model with the location of the maximum of the proton potential energy . Based on that assumption he derives a number of constraints on the possible radial variations of both the proton and the electron temperatures near the sonic point. However, even though the existence of a transonic wind seems to be intimately related to the existence of a peak of the potential , there is no reason for to coincide with the sonic point of fluid theories unless very special, and therefore unlikely, conditions are met there. For example, in the simulations presented in this paper we do always find . It can be shown analytically that this is indeed the normal case for radially decreasing temperature profiles provided T decreases more slowly than r^{-1} (Meyer-Vernet et al. 2002).
In this paper we present self-consistent kinetic simulations of a stationary solar type wind, where we concentrate on those aspects which cannot be addressed by fluid theories such as the electric field, the heat flux and the collisionality of the plasma. We deliberately treat only the most simple case of Maxwellian boundary conditions for the particles' velocity distribution function. The effects of plasma instabilities and waves are also not included in the model. Such additional "complications'' may hide part of the fundamental physics of the acceleration process and shall be discussed elsewhere. In this respect we do merely mention that the effect of resonant waves on a hybrid (fluid + kinetic) solar wind model has been discussed by Tam & Chang (1999) who conclude that ions may well be accelerated more efficiently by resonant waves rather than by the radial electric field. A similar model has been used by Lie-Svendsen & Leer (2000) to show that the two temperature electron velocity distribution functions often observed in the solar wind can be generated by Coulomb collisions without the need of assuming the presence of non-Maxwellain distribution in the corona. Olsen & Leer (1999) and Li (1999) use a closed system of transport equations based on an anisotropic bi-Maxwellian approximation for the velocity distribution functions to simulate the solar wind from the lower corona outward. Lie-Svendsen et al. (2001) extend the model down to the chromosphere based on the argument that chromosphere, transition region, corona and solar wind constitute a coupled system. The system of equations used by these authors is known as the 16-moment approximation (e.g., Demars & Schunk 1979) is a fluid-type model including transport effects such as heat flow and viscosity and even Coulomb collisions between interacting species. The main purpose of the above authors was to reproduce as good as possible the characteristics of the coronal plasma by including ad-hoc heating functions supposed to mimic the local deposit of energy due to plasma waves. If suitably chosen the heating functions can reproduce the temperature profile and temperature anisotropies which observations suggest to prevail in the solar corona (e.g., Esser et al. 1999). In many respects, our model is much more limited than the above multi-moment models which include most of the ingredients (e.g. radiation, plasma heating through waves, collisions, etc.). However, these models are fundamentally fluid models and many ingredients are not self-consistent. Our model is kinetic and fully self-consistent, but neither waves nor radiation and not even the lower layers of the corona are taken into account. In addition, because of computational limitations, we use an artificially low proton to electron mass ratio, and an exceedingly high coronal temperature, so that transposition of our results to the case of the real Sun must be done critically. One substantial difference between our results and the results from the mentioned multi-moment models is that we find that the electron heat flux in the corona is one order of magnitude larger than the classical value predicted by the Spitzer-Härm formula (Spitzer & Härm 1953), whereas the multi-moment models find it to be of the same order. Of course, both models are subject to their own limitations so that the question of whether the heat flux in the solar wind is classical or not appears to remain an open question.
Even though the physical parameters characterizing the wind simulated in the following section do not correspond exactly to those observed for the Sun, we shall refer to the simulated wind as the solar wind and to the central star as the Sun.
Details of the simulation model have been given in two previous papers
(Pantellini 2000; Landi & Pantellini 2001)
and shall not be repeated here to full extent.
The model is spatially one dimensional, i.e. all fields
depend on the heliocentric distance r only.
An equal number of protons and
electrons are allowed to move freely in the domain
,
where r_{0} is the solar
radius and
is the outer boundary of the system
located several solar radii beyond the sonic
point. The equations of motion are those of
a particle of mass m and charge q in a central gravitational field
produced by a star of mass M and a radial,
charge neutralizing electric field,
,
i.e.
The physical state of the solar corona at heliocentric distance
r_{0} (the lower boundary of our simulation domain)
is characterized by the dimensionless parameter
defined as
(
is the Boltzmann constant)
(3) |
The equations of motion Eqs. (1) and (2) are integrated for N protons, and an N electrons in the radial distance range between r=r_{0} and . The number N is determined by the requirement of the collision frequency of an electron in the system (near r=r_{0}) being roughly equal to the Fokker-Planck collision frequency of a plasma with an electron number density , which is a typical figure in the solar corona. The so calculated number N turns out to be of the order 10^{3} (Landi & Pantellini 2001). Each time a proton or an electron hits the boundary at r=r_{0} it is injected back into the system according to a non drifting isotropic Maxwellian velocity distributions with a temperature T_{0}. Given that the protons reaching the top boundary at are generally either supersonic or nearly supersonic most of it must be re-injected into the system at r=r_{0}. On the other hand, electrons reaching the boundary are injected back into the system either at r=r_{0} or at depending on what is needed to make the electron flux to be equal to the proton flux (zero charge current condition). The injection method ensures that there are always N protons and N electrons in the system. The velocity distribution of the electrons injected at the top boundary is chosen to be a drifting bi-Maxwellian with radial and perpendicular temperatures equal to the radial and perpendicular temperature of the outgoing electrons. Finally, the drift velocity of the electrons at is taken to be equal to the drift velocity of the protons measured at . In case of a supersonic wind this is just the average velocity of the protons escaping from the top of the system. In a subsonic wind some protons have to be re-injected into the system from the top. The temperature and the number of the re-injected protons is then adjusted iteratively until a coherent solution is obtained, in a manner similar to the one used by Landi & Pantellini (2001) to simulate the static corona. The electric field profile is adjusted iteratively, during an initialization phase, until zero charge flux and local charge neutrality is achieved in all points of the system.
The number density n and the collisionality
of the simulated plasma is dependent on the number
of particles N used in the model.
Figure 1 shows the results of 4 simulations which only differ
in the number N of simulated particles, i.e.
N=400, 784, 1600, 6400, the
N=400 run being the one with the most tenuous (i.e. less collisional)
atmosphere. The corresponding number densities at the base of the
system are
and 13.4, respectively.
The differences between the 4 simulations are substantial
in many respects. The most evident difference is that the wind acceleration is
much more efficient in the high density case, even though the thermal
Knudsen number
(where
is the electron-proton collisional mean free path
based on the electron-proton rate of momentum exchange
(Landi & Pantellini 2001, Appendix B), and where
is the radial
thermal electron velocity)
is much smaller than unity for all runs, ranging from
10^{-3}, for the densest case, to 10^{-2} for the most tenuous case.
From the figure it appears that the two more tenuous cases do not even
become supersonic with respect to radial proton thermal velocity
(which coincides with the fluid isothermal sound speed when
). The curves
on the bottom panel of Fig. 1
illustrate the effect of collisions on the proton potential energy
.
The latter results from the sum of the gravitational potential and a
charge neutralizing electrostatic potential :
Figure 1: Flow velocity, Mach number and proton potential energy profiles for four different values of the number of particles N in the system. From bottom to top the velocity profiles correspond to N=400, 784, 1600, 6400, respectively. In all runs and . The normalizing velocity v_{0} is the proton thermal velocity . TheMach number is defined as the ratio of the radial bulk velocity ofthe plasma divided by the proton thermal speed . The normalizing energy is given by . Thus, as a reference, if the charge and current neutralizing electric field was the Pannekoek-Rosseland potential for . | |
Open with DEXTER |
In summary, the main consequence of increasing the plasma
density beyond some threshold appears to be a reduction
of the potential barrier the protons
have to overcome in order to escape to infinity, accompanied by
the formation of a local maximum in the
profile.
As we shall see below the formation of the maximum in the
proton potential energy is intimately related to the
existence of both an outward directed, and radially
decreasing heat flux, and a radially decreasing temperature profile.
Both contribute in strengthening the outward directed electric field
,
thus facilitating the extraction of the protons.
In this context we shall remember that if the plasma was static
(impermeable boundaries) with equal electron and proton temperatures,
the charge neutralizing potential
would be the
celebrated Pannekoek-Rosseland potential (e.g., Rosseland 1924)
Let us now address the question of the wind acceleration mechanism.
In order to do so we write the energy equation for the species j
for the case of a steady state and purely radial wind.
Indeed, the second moment of Boltzmann's equation (e.g., Endeve & Leer 2001)
leads to the following expression
Figure 2: Proton and electron energy profiles obtained by evaluating Eq. (6) for the N=6400 simulation. Note how both, and are separately constant over most of the simulation domain. The electrostatic profile is plotted as a reference. | |
Open with DEXTER |
(11) |
Figure 3: Relative importance of each term in Eq. (10) for the most dense case N=6400 (top panel) and the most tenuous case N=400 (lower panel). The dense case supports a transonic wind (the vertical line in the top panel gives the position of the sonic point r_{*}) while the tenuous case does not. From the comparison of the two figures it appears that the wind acceleration is primarily driven by the heat flux term q/(nv). | |
Open with DEXTER |
Figure 4: Radial and perpendicular temperature profiles for 4 different values of the N. Note the log-log axis of the top two panels. | |
Open with DEXTER |
Figure 4 shows that the radial electron temperature profile is essentially insensitive to the density while is not. Indeed, in the collisionless limit the parallel and perpendicular temperatures are independent of each other and one should observe in order to satisfy to the angular momentum conservation law of individual particles (cf. Eq. (2)). As a consequence, in the rigorously collisionless case, should decrease by a factor 50^{2} from bottom to top of the simulation domain. Collisions do significantly contribute in limiting this bottom to top perpendicular electron temperature gradient which ranges from 10 to 30 depending on the value of N. On the other hand, for all four cases the parallel temperature only drops by a factor 3 from bottom to top, leading to strong temperature anisotropies (bottom panel in Fig. 4). This is not particularly surprising as in the collisionless limit the parallel temperature of a plasma plunged in a potential field is constant as long as the velocity distribution function is close to Maxwellian.
We can now summarize the wind acceleration scenario from a
kinetic point of view. The natural decrease of the temperature
with distance (essentially due to the
fact that in the collisionless limit
)
implies the existence of a radial heat flux
predominantly conveyed by the electrons.
The heat flux is transfered from the electrons
to the protons which become accelerated
in the outward direction.
Since the momentum of the
wind is mainly carried by the heavy protons (rather than by
the light electrons) the plasma as a whole becomes accelerated
in this way.
This mechanism must be particularly efficient in the
region located inside the spherical shell
(location of the maximum of
)
where the protons have to climb uphill in order to escape
from the potential energy well (cf. Fig. 1).
As already stated above, collisions contribute in increasing
the electric field strength. Since the electric field is directed
outward, increasing the electric field favors the extraction of the
protons from the gravitational well by reducing the height of the
maximum in the proton potential
.
The fluid estimate of the macroscopic electric field
for a
spherically symmetric electro-proton plasma can be obtained
by differentiating Eq. (6) for the electrons.
Neglecting the small terms proportional to
and taking advantage of the fact that
is approximately constant (in particular if one
compares it to the
profile shown
in Fig. 2) we then obtain
It is instructive to apply Eq. (12) to a
a weakly collisional system. In such a case the parallel temperature
is roughly constant and the perpendicular
temperature
.
This leads to
the collisionless approximation
Figure 5: Parallel electron velocity distribution functions (solid lines) in arbitrary units at two different positions in the system, near the base and in the supersonic region. Shown are results of the N=6400 simulation. The dashed lines represent Maxwellian distributions for which the first three moments (density, mean velocity and temperature) are those of the measured distributions. Velocities are normalized to the electron thermal velocity . | |
Open with DEXTER |
For the N=6400 simulation the parallel electron velocity distribution functions at the base of the system at r=r_{0}, and in the supersonic region at r=45 r_{0} are shown in Fig. 5. Since the collisional mean free path is proportional to v^{4} high energy electrons are nearly unaffected by collisions on their journey through the system. Thus, the velocity distribution of the high energy electrons flowing downward is the imprint of the upper boundary condition whereas the velocity distribution of the high energy electrons flowing upward is the imprint of the lower boundary at r=r_{0}. On the other hand, the low energy electrons, which populate the core of the velocity distribution function, are strongly affected by collisions. As a consequence, at low velocities, the electron velocity distributions are approximately isotropic Maxwellians which do not carry any heat flux. Instead, the heat flux is carried by the high energy electrons which are responsible for the asymmetry of the distribution function. As illustrated by the profiles in Fig. 5 the heat flux is due to a deficiency of downward flowing high energy electrons near the lower boundary and to an excess of upward flowing particles in the supersonic region.
Figure 6: Electron heat flux profile (solid line) and thermal Knudsen number (dashed line) for the most strongly collisional case N=6400 normalized to , where is the electron thermal velocity at the base of the system at r=r_{0}. The triangle on the heat flux axis indicates the heat flux expected near the base of the system using the Spitzer-Härm formula (Spitzer & Härm 1953). The dot-dash line shows the law normalized to the measured flux at r=r_{0}. | |
Open with DEXTER |
Figure 6 shows the heat flux and the thermal Knudsen number K_{T} measured in the N=6400 run.
One observes that while
it remains approximately
constant in the supersonic region above the
level,
K_{T} grows steeply in the subsonic region,
where it increases from 10^{-3} to 10^{-2}.
Since
in whole simulation domain, it is not
particularly surprising that the heat flux closely follows a
dependence (dot-dash curve) as in the case of the
classical Spitzer-Härm heat conduction formula (Spitzer & Härm 1953).
This conclusion is misleading, since, despite the
smallness of the Knudsen number, the heat flux
is strongly non classical.
As already pointed out by several authors in the past
the classical heat conduction
formulation breaks down either because the heat flux intensity exceeds
a value of the order
10^{-2}q_{0} (Gray & Kilkenny 1980),
or because the Knudsen
number is larger than 10^{-3} (Shoub 1983),
or because the flow velocity is a significant fraction of the sound
speed (Hollweg 1974; Alexander 1993), or
even because the electric field in the system
is of the order of the Dreicer field
(Scudder 1996b).
Indeed, in the N=6400 simulation the electric field at the sonic point is
,
reaching
an intensity
in the N=1600 case. Concerning the
heat flux intensity, Fig. 6 shows that
it is small enough for the
the low heat flux intensity condition established by
Gray & Kilkenny (1980) to
be satisfied. One can therefore conclude that
the heat flux intensity is low enough for the plasma
to be capable to support a Spitzer-Härm flux.
Let us now address the problem of the heat flux
in a flowing, and weakly collisional plasma. As discussed by
Hollweg (1974) and Alexander (1993)
a non negligible fraction of the of the energy is carried
by a collisionless term of the form
where
is a positive numerical factor of order unity
(note that the electron temperature has been supposed to be isotropic
by these authors).
Given that collisions are still relatively important in our simulations
we make the ansatz that the observed electron heat flux
is made of the sum of a classical (collisional) term
(e.g., Braginskii 1965)
and a collisionless term
Figure 7: Electron heat flux calculated using the collisionless approximation (top panel), the classical Spitzer-Härm approximation (middle panel). The lower panel shows the heat flux profiles effectively measured in the simulations. Fluxes are in arbitrary units, but the same normalization has been used for all simulations. | |
Open with DEXTER |
Given the artificially low mass ratio
in our simulations there is a concern about the
sensitivity of the results on the value of
.
In order to address this question we show the same simulation for two different
values of
in Fig. 8.
The other parameters are
identical for both simulations , i.e.
and N=6400.
In both cases the formation of a transonic wind occurs,
in association with the formation of a
maximum in the proton potential. However, the maximum's amplitude is
substantially higher, and less peaked, in the low mass ratio simulation.
The discrepancy is likely due to the fact that in the high mass
ratio case the scattering of the electrons in velocity space
by the protons is more efficient than in the low mass ratio case.
Indeed, the temperature ratio
reaches a value of 3 at the upper boundary in the
case (cf. Fig. 4) and a value of 4 in the
case. As a result the absolute value of second term on the right hand side
of Eq. (9) is significantly
smaller in the high mass ratio case than in the low mass ratio case.
Since the sign of this term is negative it contributes in reducing the
the strength of the overall positive electric field.
From Fig. 4 one may argue that
a similar argument applies to the
observation that the electric field strength increases with increasing
plasma density (i.e. with increasing collisionality) as does effectively
show Fig. 1.
Figure 8: Dependence of the Mach number and the proton potential energy on the proton to electron mass ratio for the simulation with N=6400. Even though the qualitative behavior is similar it appears that a high mass ratio is in favor of a stronger acceleration of the wind. The definitions for the Mach number and the normalizing energy are the same as in Fig. 1. | |
Open with DEXTER |
Extrapolating these observations to and N=6400 one therefore expects the maximum of the proton potential to drop to an even lower level. The peak is expected to be at least as marked as for the case.
From self-consistent kinetic simulation of a solar type wind we find that, unless an efficient isotropization mechanism for the electron velocity distribution (e.g. wave-particle interaction), or some source of suprathermal electron distributions are invoked (e.g. shock produced), the formation of a transonic wind is only compatible with a sufficiently high collisionality in the vicinity of the sonic point r=r_{*}. In oder words, the coronal density must exceed a threshold density for the wind acceleration to be sufficiently strong to become supersonic. Given the admittedly oversimplifications in our model, combined with the fact that the parameters we use are rather unrealistic (excessively high coronal temperature and low ratio) makes it impossible for us to specify an upper limit for the thermal Knudsen number at the sonic point r_{*}. We do merely show that the number cannot be arbitrarily large, the limiting value most likely being of the order unity, or less. As already stated, non Maxwellian boundary conditions, feeding an excess of suprathermal particles into the system (e.g. kappa distributions) may help overcoming the low Knudsen number condition. However, the existence of non thermal particle distributions rises the question of their origin. We chose not to address this question and assume Maxwellian boundary conditions which have the advantage of not requiring the introduction of additional ad hoc parameters into the model. Thus, in the absence of any electron scattering mechanism (apart from collisions), and unless special boundary conditions are imposed at the base of the simulation domain, collisions appear to be the essential ingredient for the wind to become accelerated to supersonic velocities, mainly because collisions are necessary to convert the electron heat flux into bulk energy of the plasma. The enthalpy gradient does also contribute to the acceleration of the wind. However, the acceleration associated with the radially decreasing enthalpy is found to be weakly dependent on the plasma collisionality and doesn't seem to be the discriminating factor in the acceleration of the wind to supersonic velocities.
In simulations where a transonic wind forms we find that the proton potential has a maximum near (but above) the sonic point. Typical values of the electric field near the sonic point r_{*} are found to be of the order of Dreicer's field, or larger. The presence of such strong electric field intensities may contribute in making the electron heat flux depart from the classical Spitzer-Härm formula (which requires the electric field being much weaker than Dreicer) but the main reason for the observed heat flux to depart from the Spitzer-Härm prediction is due to the presence of a strong "non collisional'' heat flux . The latter appears to contribute significantly to the total electron heat flux, even in the region where the wind velocity is much smaller than the sound speed. We are aware of the fact that extrapolating the above results to the "real'' Sun is a perilous exercise. However, we do not expect the qualitative behavior of a system with real coronal temperature and real proton to electron mass-ratio to behave in a substantially different way from the high density case discussed in this paper. In particular, increasing the proton to electron mass ratio from 400 to 1836 implies a factor two increase in the electron thermal velocity, only. Given that the electron thermal velocity at the base of the system is already one order of magnitude larger than the escape velocity for the case, not much difference is expected in a system with twice this thermal velocity. The skeptical reader may also argue that using a realistic mass ratio would substantially modify the transport properties of the plasma. This is certainly true, but we expect the modifications to be small, essentially because neither the classical electron heat flux nor the collisionless heat flux depend on the mass ratio, at least as long as . We expect the excessively high coronal temperature used in our simulations to be a more crucial limitation in the process of transposing the above results to the solar case just because the proton thermal velocity and the escape velocity are of the same order, i.e. . Indeed, the coronal temperature of the real Sun is roughly 3 times smaller that the value we use here. This implies a factor difference for the order unity quantity . The impact is certainly non negligible from a quantitative point of view but there aren't any reasons for us to believe that the qualitative aspects of our results do not apply to the solar case. For instance, whether or not the collisionless heat flux near the solar sonic point is really one order of magnitude stronger that the classical Spitzer-Härm flux remains an open question since this finding discords with recent results from multi-moment models.
Acknowledgements
We thank the referee for the detailed comments which very much helped us during the revision of the paper.