A&A 427, 157-167 (2004)
DOI: 10.1051/0004-6361:20035873
P. Lesaffre ^{1,2,3} - J.-P. Chièze ^{3} - S. Cabrit ^{4} - G. Pineau des Forêts ^{5,6}
1 - Institute of Astronomy, Madingley Road, Cambridge
CB3 0HA, UK
2 - University of Oxford, Department of Astrophysics,
Oxford OX1 3RH, UK
3 - CEA/DAPNIA/SAp, Orme des Merisiers, 91191
Gif-sur-Yvette Cedex, France
4 - LERMA, UMR 8112 du CNRS, Observatoire de
Paris, 61 Av. de l'Observatoire, 75014 Paris, France
5 - IAS, UMR-8617 du CNRS, Université Paris-Sud, Bât. 121, 91405 Orsay Cedex, France
6 - LUTH,
UMR-8102 du CNRS, Observatoire de Paris, 92190 Meudon Cedex, France
Received 16 December 2003 / Accepted 4 June 2004
Abstract
In the first paper of this series (Paper I)
we computed time dependent
simulations of multifluid shocks with chemistry and a transverse
magnetic field frozen in the ions, using an adaptive moving grid.
In this paper, we present new analytical results on steady-state molecular shocks. Relationships between density and pressure in the neutral fluid are derived for the cold magnetic precursor, hot magnetic precursor, adiabatic shock front, and the following cooling layer. The compression ratio and temperature behind a fully dissociative adiabatic shock is also derived.
To prove that these results may even hold for intermediate ages, we design a test to locally characterise the validity of the steady state equations in a time-dependent shock simulation. Applying this tool to the results of Paper I, we show that most of these shocks (all the stable ones) are indeed in a quasi-steady state at all times, i.e.: a given snapshot is composed of one or more truncated steady shock. Finally, we use this property to produce a construction method of any intermediate time of low velocity shocks (u < 20 km s^{-1}) with only a steady-state code. In particular, this method allows one to predict the occurrence of steady CJ-type shocks more accurately than previously proposed criteria.
Key words: magnetohydrodynamics (MHD) - shock waves - ISM: general - ISM: molecules - time
In a previous paper (Lesaffre et al. 2004, hereafter Paper I), we presented a numerical method to compute the time-dependent evolution of molecular shocks with a realistic cooling and chemistry in the presence of a transverse magnetic field. The use of a moving grid algorithm allowed us to reduce the number of zones and the computational time. However, the computation of the evolution of a stable shock from formation to steady-state still involves one or two days of CPU time on a 500 MHz workstation. Oscillating shocks require one week or two. This prevents the use of this code to fit shock parameters to observations.
On the other hand, steady-state codes give a fast answer (in one or two minutes) and include much richer physics. They have therefore been extensively used to interpret observed spectra. Nevertheless, steady-state codes have their limits. Observed magnetic molecular shocks may not be fully in steady-state yet, in which case they show a combination of C-type and J-type features. They might also be mildly or strongly unstable (Lim et al. 2002; Smith & Rosen 2003, and Paper I of this series), in which case a steady state model will have limited success.
However, different attempts have been made to circumvent these problems. In the field of shocks in supernovae remnants, Raymond et al. (1988) among others were successful in interpreting spectra with truncated steady J-shocks. Those models allowed them to account for incomplete recombination zones in a filament of the Cygnus loop. Chièze et al. (1998) discovered that the nascent magnetic precursor in a C-type shock was identical to a truncation of the steady state shock (for sufficiently late ages). Flower & Pineau des Forêts (1999) used this result to reproduce H_{2} excitation diagrams in Cepheus A West.
The aim of this paper is to rigorously test the ideas of Raymond et al. (1988) and Flower & Pineau des Forêts (1999), by clarifying the relationship between time-dependent models and steady-state models. Indeed, we will show that time-dependence is within reach of steady state codes, as long as the shock is not subject to strong instabilities.
In Sect. 2 we study the stationary equations of magnetic shocks, and derive new analytical relations. Section 2 is independent of the other sections. Then, in Sect. 3, we assess the validity of the stationary approach by using the results of our fully hydrodynamical code (Paper I). In Sect. 4, we explain how to build time-dependent models of low velocity shocks with the only help of a steady-state code. We discuss our results in Sect. 5 and sum up our conclusions in Sect. 6.
When the flow is in a steady state, time derivatives can be skipped in the steady frame. We derive here a few analytical relations valid along such steady flows.
We recall here the time-dependent monodimensional equations of
multifluid hydrodynamics with a frozen transverse magnetic
field. We put them in their conservative form:
(2) |
(3) |
(4) |
(5) |
(6) |
(7) |
Mass flux:
(9) |
(10) |
= | |||
= | (11) |
= | |||
= | (12) |
= | |||
= | (13) |
= | |||
= | (14) |
= | |||
= | (15) |
Magnetic flux:
(16) |
In each sector of a steady J or C shock, we will now get algebraic relations between the dominant conserved quantities. We tackle successively the following features, in the order in which a parcel of gas entering a magnetised shock would meet them:
We assume here and in all the following that the ionisation fraction
is very low:
(17) |
(18) |
(19) |
(20) |
(21) |
In this part of the magnetic precursor, ram pressure is directly
transferred into radiation via friction. There is no more increase
in the thermal pressure. Therefore, the neutral thermal pressure
becomes very quickly negligible against magnetic
pressure. It also remains negligible against ram pressure in the rest
of the magnetic precursor.
This leads to the following reduced set of equations:
(23) |
In the shock front, the viscous pressure is one additional unknown. We will therefore assume that we know the magnetic field at the end of the magnetic precursor, as well as the amount of energy radiated away in the precursor . Since the shock front is very tenuous, it is fair to assume that neither the magnetic field nor the integrated radiative losses will vary across it.
We can then define the new conserved fluxes for this region:
(25) |
(27) |
In addition, the post-shock velocity
can be calculated by setting
in Eq. (26), which gives
the following quadratic equation:
(29) |
(30) |
(31) |
In cases of high shock velocities (greater than 30 km s^{-1} for
the models in Paper I), a recoupling of the velocities
of charges and neutrals may happen in the relaxation layer, which
builds up an additional magnetic pressure.
After the recoupling zone, we can assume that
:
(32) |
(33) |
In an adiabatic dissociative shock, energy losses can still be
accounted for by conserved quantities since thermal energy of the gas
is used to dissociate H_{2}, i.e. converted into internal energy.
We then have:
(34) |
(35) |
For more clarity, we neglect here the effects of the charged fluid,
but it would be straightforward to include the magnetic pressure and
radiative effects using primed shock parameters like in the previous
subsections. We then obtain a modified version of
Eq. (26), including the H_{2} dissociation energy:
(36) |
(37) |
The compression factor through such a shock can now be computed:
(39) |
We can get a simple expression for the temperature of
the atomic plateau (f^{*}=0) if we neglect the post-shock ram pressure
(so that
). This is only 20% accurate for the
compression factor obtained (we use
):
The
knowledge of the compression factor, along with the assumption of
steadiness of the adiabatic front provides us with the velocity of the
front relative to the piston in the adiabatic phase of a dissociative
shock front:
We verified the analytical relations derived above by comparing with numerical magnetic shock simulations from Paper I.
First, we find that the temperature and density of the atomic plateau are well predicted by the formulae (39) and (40) in most dissociative shocks. This was expected, since the width of the adiabatic shock is so small that it is very likely to be in a steady state: the sound crossing time of this feature is much shorter than the time of variation of the entrance shock speed. In addition, in strongly oscillating shocks (weakly dissociative case of Paper I), the maximum expansion of the front, corresponding to the adiabatic, fully dissociative phase, coincides with the velocity given by Eq. (41) (see Fig. 4d of Paper I).
We also verified the relations predicted for non-dissociative shocks.
As an example, in Fig. 1, we check the relations (22) and (24) against the final steady-state of a
C-type shock (diamonds). The agreement is very good, confirming that
magnetic compression is correctly treated in the code.
Figure 1: Steady state analytic relations between the velocities of neutrals and charges are compared to an overlaid steady C-shock (diamonds). The parameters of the shock are u=20 km s^{-1}, n=10^{4} cm^{-3} and b=0.1, time is t=10^{5} yr. The velocities in the shock frame are computed using a velocity of the shock front of 0.13 km s^{-1}, inferred from Fig. 3d of Paper I. | |
Open with DEXTER |
Finally, in Fig. 2 we plot in diamonds the state of the
gas for the same shock in a snapshot at t = 100 yrs, i.e. well
prior to steady (C-type) state. We overlay the algebraic relations
previously derived for the various shock regions. The agreement turns
out to be very good. This comforts us in the ability of the code to
reproduce the conservation equations. It also points out the fact that
steady equations may well be valid to describe a shock at early times,
even though it has not reached its final steady-state. It is to
address the domain of validity of this "quasi-steady''
approximation that we set up the technique and
tests described in the next section.
Figure 2: Steady state analytic relations between pressure and mass density are compared to an overlaid future C-shock (diamonds). The parameters of the shock are u=20 km s^{-1}, n=10^{4} cm^{-3} and b=0.1, time is t=10^{2} yr. The additional necessary parameters and are read in the shock model at the end of the precursor (they are not fitted). | |
Open with DEXTER |
Consider y, one of the N+6 state variables that enter the set of
Eqs. (1)-(8). Its evolution equation in the frame of the
piston can be cast in the following form:
(42) |
We define the local steady velocity v_{y} for variable y as: the
velocity of a reference frame in which the time derivative in the
evolution equation vanishes. Hence:
(43) |
Finally, expressions (44) and (45) yield a way of characterising the "steadiness'' of the flow. At a given position x, if v_{y} does not depend on y, then there is indeed a frame moving at a velocity v(x) in which none of the variables is changing in time. Furthermore, if this velocity is constant throughout an extended region, then this whole region is moving "en-bloc'' at velocity v and can be modelled with a truncated steady-state model. We say that this region is in a quasi-steady state.
We evaluated the numerical noise in the following way: we computed
the change
in the steady velocity v_{y} when each
variable y' was changed by 10^{-4} in relative value
(corresponding to our guess for the numerical precision). We then
estimated the numerical noise
on
variable y by:
(46) |
(47) |
(48) |
(49) |
(50) |
If the numerical noise is well estimated, a value of indicates that the dispersion of invidual v_{y} values is consistent with local numerical noise, i.e. that there may exist a common steady velocity for all the variables of the subset S at that position. Regions where this is fulfilled and v is constant are the quasi-steady regions for the set S.
A less strict criterion can be chosen for the local steadiness if we think of the ratio as a "signal to noise'': even if is high, it may be possible that the ratio is high. In this case, the steady velocities corresponding to different variables are not equal, but they are close to one another: therefore, a common steady velocity v is a good approximation.
It turns out that we compute very low values of in quite a few zones. This indicates that the numerical precision we have in these zones is far better than our estimate of 10^{-4}. We hence rather use the criterion based on the ratio and define the quasi-steady regions as the regions with a constant velocity v and a good "signal to noise'' ratio.
The subset of variables S on which the averages are computed should be the whole set of variables. However, the last section of this paper needs only that the dynamically important variables be in a quasi-steady state. Therefore we present here the results for a subset S including the temperatures of the three fluids, the four velocities (roots of Eq. (45) for and ), He, H_{2}, H_{2}O, CO and OH densities. We do not include the magnetic field in magnetised shocks, because its steady velocity differs from the other variables as we will show in the next section. We also did the calculation for S including all the variables but the magnetic field, and found that it did not change the general conclusions: the signal to noise ratio is slightly less good, and a few zones are not quasi-steady anymore because of marginal chemical species having a different steady velocity than the bulk of the variables.
We summarise the results of our investigation in the next two subsections, devoted respectively to non-magnetised and magnetised shocks.
Non-dissociative J-type shocks
Figure 3 shows the
as well as v with error bars
in each zone of a typical snapshot of a non-dissociative
J-type shock. The adiabatic front and relaxation layer show a good
"signal to noise'' ratio and a flat steady velocity. On the
contrary, the pre-shock zones show a huge dispersion around the steady
velocity, expected from the homogeneity of the medium there. The fact
that we retained both velocities from expression (45) does
not alter the results, because the numerical errors on velocities are
actually much larger than for the other variables.
Figure 3: We plot the values and steady velocities for each variables of S in each zone of snapshot t=200 years for the shock with parameters b=0, n=10^{4} cm^{-3}, and u=20 km s^{-1}. Each dot represents the steady velocity of one variable computed thanks to expressions (44) and (45). The error bars are evaluated zone per zone on these values. We indicate the computational domains associated to the adiabatic front and the relaxation layer of the shock. | |
Open with DEXTER |
We find that all of our non-dissociative J-type shocks are in a quasi-steady state from the adiabatic phase to the steady phase (only the initial formation of the shock front is not quasi-steady). Therefore, a snapshot of such a shock will always coincide with the truncated structure of a J-type steady state.
In Paper I, we plotted the trajectory of the point of maximum ratio of viscous over thermal pressure (see Fig. 1d). We compute here the average steady velocity v over the whole structure of the shock (trimmed from the pre-shock values) at various times and overlay it over this trajectory in Fig. 4. Error bars show the good consistency of the test, and the correspondence of v with the velocity at which the viscous shock front moves away from the piston.
Note that the entrance velocity in the shock is not the upstream
velocity u of the fluid towards the piston, but rather u^{0} = u+v.
Hence, the entrance shock speed for the truncated steady shock is
evolving in time. In Sect. 4, we present a way to
reconstruct this evolution.
Figure 4: Trajectory and velocity away from the piston of the J-shock with parameters b=0, n=10^{4} cm^{-3}, and u=20 km s^{-1} (from Paper I). Overlaid diamonds are the steady velocities v averaged over all variables and all zones, for each snapshot analysed. | |
Open with DEXTER |
Dissociative J-type shocks
On the contrary, dissociative shocks are almost never in a
quasi-steady state. For weakly dissociative velocities, we could not
come up with a coherent picture. This was expected, since
these shocks are highly unstable with large bouncing oscillations
between a fully dissociative expansion phase and a non-dissociative
recoil phase (Paper I). Partly ionised shocks are less
unstable. We illustrate their behaviour with a typical example
shown in Fig. 5. The adiabatic front is generally in a
quasi-steady state with velocity roughly equal to the velocity of
the viscous maximum, but the dispersion is much higher
than for non-dissociative shocks. The
first plateau that follows (with H ionisation and Lyman cooling) is not
at all in a quasi-steady state: the dispersion is huge and the mean
velocity is not even constant. The second plateau (H recombination) seems rather quasi-steady, but with a velocity much
lower than the adiabatic shock front.
This warns us that steady-state diagnostics may be
hopeless for weakly non-dissociative shocks, and that we have to be
cautious for partly ionising shocks.
Figure 5: Same as Fig. 3, but for a partly ionising shock of parameters b=0, n=10^{4} cm^{-3}, and u=40 km s^{-1} at time t=220 years. Here, the scale of the plot is linear, so that the dispersion is in fact much greater than for the non-dissociative shocks, even in the second plateau. | |
Open with DEXTER |
Figure 6: Same as Fig. 3, but for time t=100 years of the future C-shock with b=0.1, n=10^{4} cm^{-3}, and u=20 km s^{-1}. We also show the steady velocity for the magnetic field. | |
Open with DEXTER |
Figure 7 shows that the two steady velocities (relaxation layer and precursor) correspond well to the velocities of the viscous maxima of neutrals and charges determined by time-derivation of their trajectory. As time evolves, the velocity of the magnetic precursor and that of the relaxation layer get closer to one another, and finally coincide after the J-front has disappeared. The C-type structure is then in a quasi-steady state as a whole.
In principle, one should then be able to model an early age of a low velocity C shock by combining a truncated C-type model with a truncated J-type model (in which the magnetic field is treated appropriately). The problem is a bit more complex than in the J-type case because we need now to determine two different truncation distances and two sets of entrance parameters, but we will show in Sect. 4.2 that it is possible to solve.
This picture holds for all shocks with magnetic
field, as long as there is no dissociation or ionisation plateau. In
the case of dissociative velocities, the same problems described in
the previous subsection arise in the corresponding features of the
relaxation layer.
Figure 7: Same as Fig. 4, but for a C-shock of parameters b=0.1, n=10^{4} cm^{-3} and u=20 km s^{-1}. Curves plot the trajectory and velocity of the neutral (solid) and charged (dashed) viscous fronts. Diamonds are the quasi-steady velocities of the relaxation layer and adiabatic front. Triangles are the quasi-steady velocities of the magnetic precursor. | |
Open with DEXTER |
This picture is the same for all magnetic shocks. The only difference between CJ-type and C-type is whether or not the J-front has disappeared when the steady velocities of the relaxation layer and the magnetic precursor converge. From this remark, we will obtain in Sect. 4.2.1 a way to assess if a low velocity shock will eventually become a steady CJ-type shock.
Section 3.2.1 has shown that the whole structure of non-dissociative J-type shocks is at all times quasi-steady. One may then safely fit truncated steady models to observations. The fitted parameters would be the entrance parameters in the shock frame and the truncation distance. But one would then like a method to relate these parameters to the parameters of the shock in the piston frame, and to the age of the shock. Conversely, one would like to build at will a snapshot of a shock of given age and parameters (in the piston frame), with the only help of a steady state code. We come up with such a procedure in the following.
A steady-state code provides us with the steady profile of any variable in the frame of the shock front for a given set of entrance parameters (u^{0},n). Say, the steady velocity , where x is the distance from the shock front. If we are given an inflow speed and density (u,n) in the frame of the piston and a time t, the problem is to find what is the entrance velocity u^{0} at the same time t in the frame of the shock front as well as the corresponding distance r between the shock front and the piston.
To help set up the notations in both frames of the shock
and of the piston, we sketched in Fig. 8 the different
lengths and velocities involved.
Figure 8: Schematical view of a J-type shock in the piston frame and in the shock frame. | |
Open with DEXTER |
At any given time, velocities in the shock frame are found by
adding
to velocities in the piston frame. The entrance
shock speed is then:
(51) |
(52) |
Conversely, if the problem is to recover the time from a steady model truncated at distance R, we only have to integrate Eq. (53) backward in time up to the point where r=0 and compute .
High compression factor approximation:
For high compression factors, the final
is small
enough to be neglected with respect to u, which makes the
computation even easier. The integral yielding the age is dominated
by the very low velocities, which are also the most recent
ones. Therefore, we only need one steady state model
.
The age of such a truncated shock is
simplified in the following way.
(54) |
Figure 9: Schematical view of an early magnetised shock in the piston frame. | |
Open with DEXTER |
The analysis of Sect. 3 showed that low-velocity magnetised shocks are composed of two quasi-steady regions: a magnetic precursor, and a non dissociative J-type feature. In faster shocks, where the entrance velocity in the J-type feature is dissociative, the J-type structure is not in a quasi-steady state, although the magnetic precursor is. We thus restrict our analysis to the non-dissociative cases.
In principle, one should then be able to model a snapshot of such a shock by gluing together two truncated steady C and J models. The problem is to determine the entrance parameters and lengths of each of the two shock features, for a given time t and a given set of parameters (u,n,B) in the piston frame. A rigorous construction method is outlined below.
In the following, and denote respectively the distance of the fronts of the C-type and the J-type features with respect to the piston (see Fig. 9). Variables computed in the reference frame of the C-type feature are denoted with C subscripts. Variables in the J-type feature frame are specified with J subscripts. Entrance parameters in the shocks have a 0 superscript.
Now, let us assume that we know the positions and at current time t. To solve for their evolution, we need to find equations that will determine and .
The entrance velocities in the C-type feature are simply:
As in the non-magnetic case, the derivatives
and
are then determined by stating that both charges and
neutrals have to be at rest near the piston, i.e., in the frame of
the J-front:
However, one could also think about solving Eqs. (56) and (58) for given arbitrary values of the final front velocity . The result would be a series of physically consistent steady CJ-type states, each characterised by a different distance . Only one of these is selected by the time evolution but, if the entrance parameters (u,n,B) are allowed to evolve in time, it might be possible that several (or even all) of these final states can be realised. The final state would then depend on the evolution history of the entrance parameters.
In fact, a very easy way to exhibit one of those CJ-type steady-states would be to use a multifluid steady-state code, and trigger the viscous dissipation in the neutral at a given position where the neutral velocity is still supersonic.
For all the low velocity cases encountered in our simulations,
we noted that after one year of time, the velocity of the charges was already
almost brought to rest at the end of the magnetic precursor. This
approximation yields the following equation:
Just like for the J-type shocks, high magnetic compression factors
will lead to
negligible before u, and will facilitate
the integration of Eq. (59). In this case, and if in addition
,
the age of the shock is given by:
(60) |
The quasi-steady state analysis of shocks opens a new field of possibilities for the steady-state codes. We compare hereafter our method to previously used algorithms, and sketch possible extensions of our method.
Our quasi-steady state analysis of J-shocks justifies the use of truncated steady-state J-shocks by Raymond et al. (1988). We provide more theoretical basis to link the true age of the shock to the truncation distance used.
Flower & Pineau des Forêts (1999) and Le Bourlot et al. (2002) use simple algorithms to produce mixed C-type and J-type features to mimic time-dependent magnetised shocks. Le Bourlot et al. (2002) greatly improved the method used by Flower & Pineau des Forêts (1999) since they keep the multifluid treatment of the flow through the relaxation layer. They just switch on viscosity in the neutral fluid when they encounter a sonic point. The present analysis gives a less heuristic way to know at which point the viscosity should be switched on, and Paper I has already shown that it can be way upstream a sonic point. Furthermore, we specify that a change of velocity frame has to be done at the end of the magnetic precursor, except for the magnetic field equation. Finally, we state where the J-type structure has to be truncated for a given time t.
Our new method should therefore lead to more accurate results, and will allow the construction of much earlier phases of magnetised shocks. It shows as well that the criterion used by Le Bourlot et al. (2002) to assess whether steady-states will be of CJ-type (occurrence of a sonic point in the neutral fluid) has to be revised. CJ-type steady states may in fact occur at lower speeds, when velocity recoupling between neutral and charges enhances magnetic compression near the piston, and slows down the precursor to the expansion speed of the J-front.
Following the same idea, if non-dynamically important species happen to be non quasi-steady, they can be post-processed in parallel to the quasi-steady time-evolution with a Lagrangian code.
Here we present an algorithm for two kinematic flows (charges and neutrals), but the same method may be implemented for more flows. Especially, the treatment of charged grains could be envisaged, in relation to the questions raised by Ciolek & Roberge (2002) and Flower & Pineau des Forêts (2003). The only caveat is that we do not yet have a consistent check for the validity of the quasi-steady assumption.
Finally, our algorithm is straightforward to apply with slowly changing input conditions (u,n,B) in the piston frame. One has only to bear in mind that if these parameters change over time-scales much shorter than the crossing time scale of the shock, then the quasi-steady assumption is very likely to be violated.
Our algorithm is based upon the quasi-steady state assumption. However, it will give results with any shock. One problem is that we still have no other way to assess the validity of the steady-state assumption than computing the time-dependent evolution with a fully hydrodynamical code.
We encountered several cases where this assumption was not realised. Strongly unstable shocks like the weakly dissociative ones violate strongly this assumption. Fortunately, they seem to happen for a very restricted range of parameters. Partly ionising shocks are very slightly unstable, and are closer to meet the quasi-steady state assumption. They might therefore be accounted for by our algorithm. Finally, quite a few magnetised shocks have unstable entrance velocities for the J-shock only at early times, and are afterwards quasi-steady at all times. These shocks may be as well within reach of our algorithm if one is ready to skip the early evolution. However, one shouldalways be cautious when a plateau with dissociated molecules appears in a steady-state computation.
An other situation where the quasi-stationary assumption may be strongly violated is the case where a dynamically important chemical specie is not in a quasi-steady state. This might happen when a dominating cooling agent varies on very short time scales. Furthermore, diffusion effects, if they turn out to be important, will destroy the quasi-steady state as well.
In a companion paper (Paper I), we produced fully time-dependent numerical simulations of molecular shocks.
In the present paper, we derived new analytical relations valid at quasi-steady state, and successfully checked them on our simulations. These relations provide useful benchmarks to test existing and future multifluid codes.
In light of the simulations run in Paper I, we investigated carefully the validity of the quasi-steady state approximation. It was found that at all times stable shocks could be accounted for by truncated steady models. We point out as well that caution has to be kept regarding the use of steady-state models for dissociative velocities.
Finally, we produced a new algorithm based on the quasi-steady state assumption. With only a steady-state code, this method is able to compute time-dependent snapshots of shocks in the presence or not of a magnetic field. Therefore, it brings time-dependence within the reach of steady models, and should greatly improve the diagnostics of observed molecular shocks. Furthermore, it provides a way of assessing the CJ nature of a magnetised shock. Finally, this algorithm can be extended to many shocks other than molecular, provided that the quasi-steady state approximation is validated.
Acknowledgements
We thank the referee (Pr. T.W. Hartquist) for having kindly accepted to refer both paper I and paper II at the same time. This work was in part supported by a European Research & Training Network (HPRN-CT-20002-00303).