text:
# A 3+1 covariant suite of Numerical Relativity Evolution Systems

## I Introduction

In a recent article  (which will be hereafter referred to as the KST paper), a wide class of hyperbolic first order formalisms has been studied with a view to Numerical Relativity applications. All these formalisms do use the well known $`3+1`$ decomposition of spacetime

$`ds^2`$ $`=`$ $`\alpha ^2dt^2`$ (1)
$`+`$ $`\gamma _{ij}(dx^i+\beta ^idt)(dx^j+\beta ^jdt)i,j=1,2,3`$ (2)
which allows one to use Einstein’s field equations as an evolution system for the metric (1), namely

$`(_t_\beta )\gamma _{ij}`$ $`=`$ $`2\alpha K_{ij}`$ (3)
$`(_t_\beta )K_{ij}`$ $`=`$ $`_i\alpha _j`$ (4)
$`+`$ $`\alpha [{}_{}{}^{(3)}R_{ij}^{}2K_{ij}^2+trKK_{ij}\widehat{R}_{ij}]`$ (5)
(ADM system). The KST paper describes many ways of obtaining a first order version of the second order ADM system by introducing a number of arbitrary parameters. Different parameter choices correspond to different ways of using the energy and momentum constraints

$`{}_{}{}^{(3)}Rtr(K^2)+(trK)^2`$ $`=`$ $`2\tau `$ (6)
$`_kK_{}^{k}{}_{i}{}^{}_i(trK)`$ $`=`$ $`S_i`$ (7)
to modify the structure of the resulting systems. The goal is to get well posed systems, so that existence, unicity and stability of the solutions can be ensured . In this sense, the KST paper can be understood as a generalization of previous works by other groups , which have in common the fact that the only independent quantities are the metric coefficients and their first derivatives.

A similar approach can be used to generalise the hyperbolic first order formalisms proposed by the Palma group , as suggested in Ref. . The main difference with the formalisms studied in the KST paper is that in this case there are three supplementary variables, independent of the metric derivatives, whose evolution equations are obtained by reversing the order of space and time derivatives in the momentum constraint (7). This use of the momentum constraint to evolve three supplementary variables can also be found in another context , where a second order system of the form (4) is considered (BSSN system).

To summarize, we have different formalisms in which the constraints (6, 7) are used to modify the structure of the ADM evolution system. The momentum constraint is used in some cases to evolve supplementary quantities, either keeping the system to be of second order  or getting first order versions . But the same momentum constraint is used instead in other cases to modify the first order versions without introducing extra quantities .

The purpose of the present article is to point out the strong relationship among all these formalisms from the theoretical point of view. With this aim, we propose a new covariant framework in which the momentum constraint is used to evolve a vector quantity, so that the properties of the resulting systems are independent of the choice of space coordinates.

The paper is organised as follows: In Section II, we introduce the vector quantity $`Z_i`$, which will replace the corresponding non-vector supplementary quantities in . The original ADM system is then generalised by making use of the vector $`Z_i`$ and the energy constraint (6). In Section III, the resulting second order system (which will be hereafter referred to as system A) is compared with the BSSN system .

In section IV, we obtain a first order version of system A (which will be hereafter referred to as system B), which is compared with the formalisms presented in . Finally, in section V, system B is further modified in order to be compared with the formalisms discussed in the KST paper.

## II Integrating the momentum constraint

Let us define the vector $`Z^i`$ as follows:

$$(_t_\beta )Z_i=\alpha [_j(K_{i}^{}{}_{}{}^{j}\delta _{i}^{}{}_{}{}^{j}trK)S_i]$$
(8)
so that it is clear that $`Z_i`$ vanishes for physical solutions, which can be defined as the ones verifying both the energy and momentum constraints (6, 7). The ”zero” vector $`Z_i`$ can then be considered as a good covariant measure of the deviation of any solution $`(\alpha ,\gamma _{ij},K_{ij})`$ of the ADM system from the submanifold of physical solutions, as far as the momentum constraint is concerned. This means that the cumulative effect of the differential constraint (7) can be allowed for by looking at the simple algebraic equation

$$Z_i=0$$
(9)
which is to be regarded as the integral version of (7).

But these $`Z_i`$, besides being useful auxiliary quantities, can be also considered as supplementary dynamical variables. This is the kind of approach proposed in Ref. , where a ”$`\lambda `$-system” was proposed containing as much as 40 supplementary quantities to be evolved by using the constraints. We will limit ourselves here to only three quantities $`Z_i`$ (corresponding to $`\lambda _i`$ in ), the ones related with the momentum constraint. We will use these quantities $`Z_i`$, together with the energy constraint (6) to implement covariant quasilinear modifications of the principal part of the ADM system, without affecting the submanifold of physical solutions (the ones verifying (6, 9)). We will call ”system A” the most general second order evolution system obtained in this way, namely

$`(_t_\beta )K_{ij}`$ $`=`$ $`_i\alpha _j`$ (10)
$`+`$ $`\alpha [{}_{}{}^{(3)}R_{ij}^{}2K_{ij}^2+trKK_{ij}\widehat{R}_{ij}]`$ (11)
$``$ $`{\displaystyle \frac{\alpha n}{4}}[Rtr(K^2)+tr^2(K)2\tau ]\gamma _{ij}`$ (12)
$`+`$ $`{\displaystyle \frac{\alpha \mu }{2}}(_iZ_j+_jZ_i){\displaystyle \frac{\alpha \nu }{2}}(Z)\gamma _{ij}`$ (13)
where $`\mu `$, $`\nu `$ and $`n`$ are arbitrary parameters.

The modified evolution system (8,10) (system A) is of a mixed type: first order in ($`Z_i,K_{ij}`$) but second order in ($`\alpha ,\gamma _{ij}`$). In this sense, it is very similar to the BSSN system . In order to compare these systems, let us state some definitions:

a) Two evolution systems will be said to be equivalent if one can be obtained from the other by a change of variables. It is clear that two equivalent systems have the same space of solutions.

b) Two evolution systems will be said to be quasiequivalent when their principal parts are equivalent. It follows that two quasiequivalent systems have the same causal structure.

For instance, if we introduce ”damping terms” to modify eq. (8) (as proposed in  with a view to Numerical Relativity applications), we are not affecting the principal part (damping terms are algebraic in $`Z_i`$), so that the resulting ”damped” system would be quasiequivalent to system A. Also, if we introduce additional terms of the form $`Z_iZ_j`$ or $`Z^2\gamma _{ij}`$ in (10), the resulting modified system would be quasiequivalent to system A as well.

On the contrary, different values of the parameters $`\mu `$, $`\nu `$, $`n`$ in (10) do correspond to different forms of the principal part, so that they would lead to different evolution systems which are not quasiequivalent, even if they have (as they do) the same submanifold of physical solutions. In this context, it is natural to use system A as a framework to compare different formalisms that have been proposed with a view to Numerical Relativity applications, where the causal structure of the evolution system is important to ensure the well posedness of the problem, namely the existence, unicity and stability of the solutions.

We can also take advantage of the covariance of system A under coordinate transformations of the generic form:

$`t^{}`$ $`=`$ $`f(t)`$ (14)
$`x_{}^{i}{}_{}{}^{}`$ $`=`$ $`g(x^i,t)`$ (15)
(3+1 covariance) so that we can choose $`\beta ^i=0`$ (normal coordinates) in what follows without any loss of generality.

## III Comparing with the BSSN system

Let us consider here the system introduced in by Shibata and Nakamura  and Baumgarte and Shapiro  (BSSN system). The metric coefficients $`\gamma _{ij}`$ are expressed in terms of a conformal metric:

$$\stackrel{~}{\gamma }_{ij}=e^{4\varphi }\gamma _{ij}$$
(16)
with unit determinant, so that

$$e^{4\varphi }=\gamma ^{1/3}=[det(\gamma _{ij})]^{1/3}$$
(17)
The second fundamental form $`K_{ij}`$ is then decomposed into its trace and trace-free components, namely

$`K`$ $`=`$ $`\gamma ^{ij}K_{ij}`$ (18)
$`\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\varphi }(K_{ij}{\displaystyle \frac{1}{3}}K\gamma _{ij})`$ (19)
so that the principal part of their evolution equations can be written as :

$`_tK`$ $`=`$ $`\gamma ^{ij}_i\alpha _j+\mathrm{}`$ (20)
$`_t\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\varphi }(_i\alpha _j+\alpha R_{ij})^{TF}+\mathrm{}`$ (21)
where the superscript TF denotes the trace-free part of a tensor.

The principal part of the Ricci tensor in (21) is decomposed as follows :

$`R_{ij}`$ $`=`$ $`2_{}^{2}{}_{ij}{}^{}\varphi {\displaystyle \frac{1}{2}}\stackrel{~}{\gamma }^{rs}\stackrel{~}{\gamma }_{ij,rs}+\stackrel{~}{\gamma }_{k(i}_{j)}\stackrel{~}{\mathrm{\Gamma }}^k+\mathrm{}`$ (22)
where the “conformal connection” quantities $`\stackrel{~}{\mathrm{\Gamma }}^i`$ can be defined as:

$$\stackrel{~}{\mathrm{\Gamma }}^i=\stackrel{~}{\gamma }_{}^{ij}{}_{,j}{}^{}$$
(23)
Up to now we have only performed a transformation of variables of the original ADM system. But notice that the system (20, 21) contains seven evolution equations, compared with the six equations in (4). The extra quantity is the trace of $`\stackrel{~}{A}_{ij}`$, which is supposed to vanish by construction. This opens the way to two different approaches:

a) to consider $`tr(\stackrel{~}{A})`$ as an independent variable, which can take non-zero values during numerical evolution. The equation

$$tr(\stackrel{~}{A})=0$$
(24)
is to be considered as an additional constraint which can be used to monitor the accuracy of the simulation (free evolution approach).

b) to enforce the vanishing of $`tr(\stackrel{~}{A})`$, by removing systematically any contribution that could arise from truncation errors in the discretization of equation (21), so that $`tr(\stackrel{~}{A})`$ is not considered as a new dynamical variable (constrained evolution approach)

In what follows, we will adopt the second approach, so that the system (20, 21) is then quasiequivalent to the ADM system.

The BSSN system does consider the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ as independent dynamical variables, so that the full set of independent BSSN variables is given by

$$u=\{\alpha ,\varphi ,\stackrel{~}{\gamma }_{ij},K,\stackrel{~}{A}_{ij},\stackrel{~}{\mathrm{\Gamma }}^i\},$$
(25)
where we must remember that both the trace of $`\stackrel{~}{A}_{ij}`$ and the determinant of $`\stackrel{~}{\gamma }_{ij}`$ do have fixed, non-dynamical, values. An evolution equation for the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ could be derived  by permuting a time derivative with the space derivative in (23)

$$_t\stackrel{~}{\mathrm{\Gamma }}^i=2\alpha \stackrel{~}{A}_{}^{ij}{}_{,j}{}^{}+\mathrm{}$$
(26)
But the BSSN system uses instead the momentum constraint (7) to transform eq. (26) into

$$_t\stackrel{~}{\mathrm{\Gamma }}^i=\frac{4}{3}\alpha \stackrel{~}{\gamma }^{ij}K_{,j}+\mathrm{}$$
(27)
Now we are in position to compare the BSSN system with the one introduced in the previous section (system A), with a set of variables given by

$$u=\{\alpha ,\gamma _{ij},K_{ij},Z_i\}$$
(28)
The conformal decomposition (16, 17, 18, 19) can be easily performed on system A as well. The key point is then the passage from (26) to the BSSN evolution equation (27). In the framework of system A, the momentum constraint is just the right-hand-side of the evolution equation for $`Z_i`$. It follows that the use of the momentum constraint to transform (the principal part of) the right-hand-side of (26) into (27) does correspond to the following transformation of $`\stackrel{~}{\mathrm{\Gamma }}^i`$:

$$\stackrel{~}{\mathrm{\Gamma }}_i=\stackrel{~}{\gamma }_{ik}\stackrel{~}{\gamma }_{}^{kj}{}_{,j}{}^{}+2Z_i$$
(29)
which is obviously consistent with the previous expression (23) if we remember that $`Z_i=0`$ for the physical solutions.

Now we can compare easily (21) with the trace-free part of (10): it follows that they coincide if and only if $`\mu =2`$. The same comparison can be performed again between (20) and the trace part of (10): it follows that they coincide if and only if $`\nu =n=\frac{4}{3}`$.

To summarize, we have shown then that the BSSN system is quasiequivalent to system A for the parameter choice

$$\mu =2,\nu =n=4/3$$
(30)
provided that the “conformal connection” quantities $`\stackrel{~}{\mathrm{\Gamma }}^i`$ in the BSSN system are related with the vector $`Z_i`$ by equation (29).

## IV A first order evolution system

A first order version of the system A can be obtained in the standard way by considering the first space derivatives

$`A_k`$ $`=`$ $`_k(\mathrm{ln}\alpha )`$ (31)
$`D_{kij}`$ $`=`$ $`{\displaystyle \frac{1}{2}}_k\gamma _{ij}`$ (32)
as independent dynamical quantities whose evolution equations are given by

$`_tA_k+_k(\alpha Q)`$ $`=`$ $`0`$ (33)
$`_tD_{kij}+_k(\alpha K_{ij})`$ $`=`$ $`0`$ (34)
where the quantity $`Q`$, defined as

$$_t\mathrm{ln}\alpha =\alpha Q$$
(35)
can be related with $`trK`$ in order to fix the time coordinate gauge, namely

$$Q=ftrK$$
(36)
where $`f`$ is an arbitrary function of the lapse $`\alpha `$.

This allows to transform system A into a first order evolution system for the set of variables

$$u=\{\alpha ,\gamma _{ij},K_{ij},A_i,D_{kij},Z_i\}$$
(37)
with a non-trivial principal part given by (8, 33, 34) and

$$_tK_{ij}+_k(\alpha \lambda _{}^{k}{}_{ij}{}^{})=\mathrm{}$$
(38)
where we have noted

$`\lambda _{}^{k}{}_{ij}{}^{}=D_{}^{k}{}_{ij}{}^{}`$ $`+`$ $`{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{i}{}^{}(A_j+D_j2D_{rj}^{}{}_{}{}^{r}\mu Z_j)`$ (39)
$`+`$ $`{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{j}{}^{}(A_i+D_i2D_{ri}^{}{}_{}{}^{r}\mu Z_i)`$ (40)
$``$ $`{\displaystyle \frac{n}{2}}(D^kD_{}^{rk}{}_{r}{}^{})\gamma _{ij}+{\displaystyle \frac{\nu }{2}}Z^k\gamma _{ij}`$ (41)
and

$$D_k=\gamma ^{ij}D_{kij}.$$
(42)
This first order system will be called system B in what follows.

System B is very similar to the Bona-Masso system  (we will follow here the notation of Ref.  to avoid confusion), which is expressed in terms of the set of variables

$$u=\{\alpha ,\gamma _{ij},K_{ij},A_i,D_{kij},V_i\}$$
(43)
where the quantities $`V_i`$ can be defined as

$$V_i=D_iD_{ji}^{}{}_{}{}^{j}$$
(44)
An evolution equation for the $`V_i`$ could be derived by combining (44) with the evolution equations (34) for the $`D_{kij}`$

$$_tV_i=\alpha _j(K_{i}^{}{}_{}{}^{j}trK\delta _{i}^{}{}_{}{}^{j})+\mathrm{}$$
(45)
but the Bona-Masso system uses instead the momentum constraint to transform (45) into an evolution equation with a trivial principal part, namely

$$_tV_i=\mathrm{}$$
(46)
The remaining evolution equations in the Bona-Masso system are identical to those in system B, with the replacement of (39) by

$`\lambda _{}^{k}{}_{ij}{}^{}`$ $`=`$ $`D_{}^{k}{}_{ij}{}^{}+{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{i}{}^{}(A_j+2V_jD_j)`$ (47)
$`+`$ $`{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{j}{}^{}(A_i+2V_iD_i){\displaystyle \frac{n}{2}}V^k\gamma _{ij}`$ (48)
In order to compare system B with the Bona-Masso system, the key point is again the passage from (45) to the Bona-Masso evolution equation (46). The use of the momentum constraint to transform (the principal part of) the right-hand-side of (45) into (46)) does correspond to the following transformation of $`V_i`$:

$$V_i=D_iD_{ji}^{}{}_{}{}^{j}Z_i$$
(49)
which is again consistent with (44) because $`Z_i`$ vanishes for physical solutions.

Now we can substitute (49) into (47) and compare with (39): it follows that they coincide if and only if

$$\mu =2,\nu =n$$
(50)
To summarize, we have shown that the Bona-Masso system is quasiequivalent to system B for the parameter choice (50), provided that the supplementary quantities $`V_i`$ in the Bona-Masso system are related with the vector $`Z_i`$ by equation (49). As a consequence of this fact, system B inherits (for the parameter choice (50)) the causal structure of the Bona-Masso system, which is known to be strongly hyperbolic (real characteristic speeds and a complete set of eigenfields ) if and only if $`f>0`$.

## V Reducing variables space

System B has two sets of non-trivial standing eigenmodes (the ones with zero characteristic speed in normal coordinates). The first one is related with gauge evolution, namely

$$A_ifD_i$$
(51)
while the second one is the combination (49), related with the momentum constraint. We will see how we can use this fact to reduce the set of independent dynamical variables (37) of the system B to the smaller one used in the KST paper, namely

$$u=\{\gamma _{ij},K_{ij},d_{kij}\}$$
(52)
where, in spite of the fact that

$$d_{kij}=2D_{kij}$$
(53)
we will keep using here the KST original notation to avoid confusion.

First of all, one can explicitly integrate the relationship (35, 36) between the lapse and the volume element in the case $`f=2\sigma =constant`$, namely

$$_t(\alpha \gamma ^\sigma )=0$$
(54)
so that the value of $`\alpha `$ can be defined in terms of $`\gamma `$ for every initial condition (densitized lapse). The same thing can be done with $`A_i`$ and $`D_i`$, so that

$$A_i2\sigma D_i+\mathrm{}$$
(55)
The list of dynamical variables is then reduced from (37) to

$$u=\{\gamma _{ij},K_{ij},D_{kij},Z_i\}$$
(56)
where the gauge-related standing modes (51) have just disappeared from the dynamical system.

In doing this, we have restricted ourselves to the $`f=constant`$ case, in spite of the fact that this case is known  to be prone to numerical gauge instabilities if $`f1`$. We do that with a view to mimic the KST system, where $`\sigma `$ is defined to be a constant and, as a consequence of it, only examples with $`\sigma =1/2`$ ($`f=1`$) are provided. In constrast, system B would allow more general gauge choices, like the ”1+log” one ($`f=2/\alpha `$) which is frequently used in Numerical Relativity.

Disposing of the supplementary variables $`Z_i`$ in (56) is not so easy. We will take as a guide the KST evolution equations for $`\gamma _{ij}`$ and $`K_{ij}`$, which are very similar to the corresponding ones in system B: we need just to replace (39) by

$`2\lambda _{}^{k}{}_{ij}{}^{}`$ $`=`$ $`d_{}^{k}{}_{ij}{}^{}+{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{i}{}^{}[(1+2\sigma )d_j2d_{rj}^{}{}_{}{}^{r}]`$ (57)
$`+`$ $`{\displaystyle \frac{1}{2}}\delta _{}^{k}{}_{j}{}^{}[(1+2\sigma )d_i2d_{ri}^{}{}_{}{}^{r}]{\displaystyle \frac{n}{2}}(d^kd_{}^{rk}{}_{r}{}^{})\gamma _{ij}`$ (58)
where we have taken $`\zeta =1,\gamma =n/2`$ in the original KST expression to reproduce (39) more closely.

The main difference between system B and the KST one comes instead from the evolution equations for $`d_{kij}`$, namely

$`_t`$ $`d_{kij}=2\alpha _kK_{ij}+\chi \alpha \gamma _{ij}(_rK_{}^{r}{}_{k}{}^{}_ktrK)`$ (59)
$`+`$ $`{\displaystyle \frac{\eta }{2}}\alpha [\gamma _{ki}(_rK_{}^{r}{}_{j}{}^{}_jtrK)+\gamma _{kj}(_rK_{}^{r}{}_{i}{}^{}_itrK)]+\mathrm{}`$ (60)
where $`\eta ,\chi `$ are arbitrary parameters. These equations are obtained by using the momentum constraint (7) to modify the standard equations (34) which were used in system B. The key point again is to realize that the use of the momentum constraint to transform (the principal part of) the right-hand-side of (34) into that of (59) amounts, allowing for (8), to modify the relationship (53) between $`D_{kij}`$ and $`d_{kij}`$ in the following way:

$$2D_{kij}=d_{kij}\eta \gamma _{k(i}Z_{j)}\chi Z_k\gamma _{ij}$$
(61)
which is consistent again with (32) because $`Z_i`$ vanishes for physical solutions. A straightforward substitution of (55) and (61) into (39) shows that the vector $`Z_i`$ disappears from the principal part if we just relate the free parameters in system B with those of the KST system as follows:

$`f`$ $`=`$ $`2\sigma `$ (62)
$`\nu `$ $`=`$ $`\chi n(\chi \eta /2)`$ (63)
$`\mu `$ $`=`$ $`\eta \chi /2\sigma (\eta +3\chi )`$ (64)
so that the Einstein-Christoffel system  (corresponding to $`\sigma =1/2`$, $`\eta =4`$, $`\gamma =\chi =0`$ in the KST paper) is recovered when

$$f=1,\mu =2,\nu =n=0$$
(65)
To summarize, we have shown here that the KST systems are quasiequivalent (in the case $`\zeta =1`$) to a first order system which can be obtained by reducing dynamical variables in system B to the minimal set (52). The reduced system will be called system B’ in what follows. The four parameters of system B’ can be expressed in terms of those of the KST system by (62), where we must remember that $`n=2\gamma `$.

## VI Concluding remarks

The evolution systems A, B and B’ presented in this paper are closely related in a transparent way, up to the point that we have been able to keep the same set of parameters ($`f,\mu ,\nu ,n`$) in the three cases. The fact that, for suitable parameter choices, they are quasiequivalent respectively to the BSSN , Bona-Masso  and Einstein-Christoffel  systems points out the strong relationship between these systems, proposed independently in different contexts by different authors, all of them extremely useful in Numerical Relativity applications.

There are at least two specific points in which the suite A, B, B’ can improve our understanding of the previous systems. The first one is that the supplementary quantities $`Z_i`$ in systems A and B are three-dimensional vectors, contrary to what happens with their counterparts $`\stackrel{~}{\mathrm{\Gamma }}^i,V_i`$ respectively, which do not have tensor behaviour. Then, systems A and B can be used to build up covariant counterparts of the second order BSSN  and first order Bona-Masso  systems, respectively.

The second point can be seen when analysing the list of non-trivial characteristic speeds in the KST systems, namely

$$\{\pm 1,\pm c_1,\pm c_2,\pm c_3\}$$
(66)
where in our case ($`\zeta =1`$) we have

$`c_{1}^{}{}_{}{}^{2}`$ $`=`$ $`2\sigma `$ (67)
$`c_{2}^{}{}_{}{}^{2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(2\eta 2\eta \sigma \chi 6\chi \sigma )`$ (68)
$`c_{3}^{}{}_{}{}^{2}`$ $`=`$ $`1+2\gamma \gamma \eta +\chi +2\gamma \chi `$ (69)
which can be translated, allowing for (62) into

$$c_{1}^{}{}_{}{}^{2}=f,c_{2}^{}{}_{}{}^{2}=\mu /2,c_{3}^{}{}_{}{}^{2}=1+\nu n$$
(70)
where we have used $`n=2\gamma `$.

Expression (70) allows to decouple the parameter set in a way such that every parameter is related with just one characteristic speed, making the physical interpretation more transparent: the gauge parameter f is associated only with $`c_1`$ (gauge speed), the parameter $`\mu `$ (which appears in the terms affecting to the longitudinal part of $`K_{ij}`$) is associated only with $`c_2`$, whereas the two parameters $`\nu `$ and $`n`$ appearing in the transverse trace terms are related with $`c_3`$ (notice that there is some degree of redundancy between $`\nu `$ an $`n`$).

It is not surprising then that in all the particular sets of parameters corresponding to previous works  (30, 50, 65) one has

$$\mu =2,\nu =n$$
(71)
which amounts to the “physical speed” requirement for the degrees of freedom not related with the gauge ($`c_{2}^{}{}_{}{}^{2}=c_{3}^{}{}_{}{}^{2}=1`$). We are aware that in the case of the BSSN system we can not speak properly about characteristic speeds, but the close relationship between systems A and B can bring some light even in this case (see for instance Ref.  for a similar approach).

The question arises of wether systems A, B, B’ will perform any better than their well known counterparts (BSSN, Bona-Masso, KST) in numerical simulations. We have done preliminary tests, comparing system B against the Bona-Masso one in black-hole spacetimes. Our results indicate that the Bona-Masso code performs better, because of the fact that the additional quantities $`V_i`$ evolve with an ordinary differential equation (46), avoiding the truncation errors inherent to the discretization of the partial derivatives that appear instead in the evolution equation (8) for $`Z_i`$. The fact that $`Z_i`$ is a true tridimensional vector does not help very much in that context. Something similar happens when comparing system A with BSSN: the evolution equation (27) for the $`\stackrel{~}{\mathrm{\Gamma }}^i`$ is still simpler than (8). In our opinion, the BSSN, Bona-Masso and KST systems are the state of the art (each one in its class) and will be difficult to beat by any other quasiequivalent system. One would need to provide something more different, modifying even the causal structure of the system (which is given by the principal part), in order to improve in a significant way the performance (stability and accuracy) of the current numerical simulations.

Acknowledgements: This work has been supported by the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract (HPRN-CT-2000-00137), by the Spanish Ministerio de Ciencia y Tecnologia through the research grant number BFM2001-0988 and by a grant from the Conselleria d’Innovacio i Energia of the Govern de les Illes Balears.


format: text

text:
# The Chemical Composition of the 30 Doradus Nebula Derived from VLT Echelle Spectrophotometry1footnote 11footnote 1Based on observations collected at the European Southern Observatory, Chile, proposal number ESO 68.C-0149(A).

## 1 Introduction

The determination of the chemical composition of H II regions has been paramount for the study of the chemical evolution of galaxies and for the determination of the primordial helium abundance, $`Y_p`$. In recent times the determination of atomic data of higher accuracy and the detection of fainter emission lines with the use of echelle spectrophotometry have permitted to derive more accurate physical conditions for Galactic H II regions (e.g., Esteban et al., 1998, 1999a, 1999b). For these reasons it was decided to carry out echelle spectrophotometry of 30 Doradus.

Due to its proximity, its high angular dimensions, and its high surface brightness 30 Doradus, NGC 2070, is the most spectacular extragalactic H II region and thus it has been the subject of many spectrophotometric studies (e.g., Peimbert, & Torres-Peimbert, 1974; Aller et al., 1974; Dufour, 1975; Pagel et al., 1978; Boeshaar et al., 1980; Dufour, Shields, & Talbot, Dufour, et al.1982; Shaver et al., 1983; Mathis, Chu, & Peterson, 1985; Rosa & Mathis, 1987; Vermeij & van der Hulst, 2002; Tsamis et al., 2002).

The main aim of this paper is to make a new determination of the chemical abundances of 30 Doradus including the following improvements over previous determinations: the consideration of the temperature structure that affects the helium and heavy elements abundance determinations, the derivation of the O and C abundances from recombination line intensities of very high accuracy, the consideration of the collisional excitation of the triplet He I lines from the $`2^3`$S level by determining the electron density from many line intensity ratios, and the study the $`2^3`$S level optical depth effects on the intensity of the triplet lines by observing a large number of singlet and triplet lines of He I.

In sections 2 and 3 the observations and the reduction procedure are described. In section 4 temperatures and densities are derived from eight and six different methods respectively; also in this section, four independent values of the mean square temperature fluctuation, $`t^2`$, are determined by combining the electron temperatures. In section 5 ionic abundances are determined based on recombination lines that are almost independent of the temperature structure, and ionic abundances based on ratios of collisionally excited lines to recombination lines that do depend on the temperature structure of the nebula. In section 6 the total abundances are determined and compared with those of NGC 346 (the most luminous H II region in the SMC), the Orion nebula, M17, and the Sun. In section 7 $`\mathrm{\Delta }Y/\mathrm{\Delta }O`$ and $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$ are determined; these ratios are important restrictions for the study of the chemical evolution of galaxies and for the determination of the primordial helium abundance, $`Y_p`$. Also in section 7 $`Y_p`$ is determined based on the abundances of 30 Dor and a value of $`\mathrm{\Delta }Y/\mathrm{\Delta }O`$ derived from the observations of other objects and from chemical evolution models of irregular galaxies.

## 2 Observations

The observations were obtained with the Ultraviolet Visual Echelle Spectrograph, UVES (D’Odorico et al., 2000), at the VLT Kueyen Telescope in Chile. We observed simultaneously with the red and blue arms in two settings, covering the region from 3100 Å to 10360 Å (see Table 1). The wavelength regions 5783–5830 ÅÅ and 8540–8650 ÅÅ were not observed due to the separation between the two CCDs used in the red arm. There were also two small gaps that were not observed, 10084–10088 ÅÅ and 10252–10259 ÅÅ, because the two redmost orders did not fit within the CCD. In addition to the long exposure spectra we took 60 second exposures for the four observed wavelength ranges to check for possible saturation effects.

The slit was oriented east-west and the atmospheric dispersion corrector (ADC) was used to keep the same observed region within the slit regardless of the air mass value. The slit width was set to 3.0” and the slit length was set to 10” in the blue arm and to 12” in the red arm; the slit width was chosen to maximize the S/N ratio of the emission lines and to maintain the required resolution to separate most of the weak lines needed for this project. The FWHM resolution for the 30 Dor lines at a given wavelength is given by $`\mathrm{\Delta }\lambda \lambda /8800`$. The reductions were made for an area of 3” $`\times `$ 10” located 64” N and 60” E of HD 38268, these coordinates are similar to those of 30 Dor II, position observed by Peimbert, & Torres-Peimbert (1974), but the area observed by us is considerably smaller. The slit was placed in a region of high emission measure that does not include luminous stars.

The spectra were reduced using the IRAF222IRAF is distributed by NOAO, which is operated by AURA, under cooperative agreement with NSF. echelle reduction package, following the standard procedure of bias subtraction, aperture extraction, flatfielding, wavelength calibration and flux calibration. For flux calibration the standard star EG 247 was observed.

## 3 Line Intensities and Reddening Correction

Line intensities were measured integrating all the flux in the line between two given limits and over a local continuum estimated by eye. In the few cases of line-blending, the line flux of each individual line was derived from a multiple Voigt profile fit procedure. All these measurements were carried out with the splot task of the IRAF package.

An initial reddening coefficient, $`C`$(H$`\beta `$), was determined by fitting the observed $`I`$(H$`\beta `$)/$`I`$(H Balmer lines) ratios (with the exception of Balmer $`\alpha `$) to the theoretical one computed by Storey & Hummer (1995) for $`T_e`$ = 9,000 K and $`N_e`$ = 300 cm-3, (see below) and assuming the extinction law of Seaton (1979). With this initial fit the the observed $`I`$(H$`\alpha `$)/ $`I`$(H$`\beta `$) and the $`I`$(Paschen lines)/ $`I`$(H$`\beta `$) ratios become larger than predicted. A similar situation prevails in the Orion Nebula (e.g., Costero & Peimbert, 1970; Cardelli & Clayton, 1988; Greve, Castles, & McKeith, 1994).

A simultaneous fit to all the H lines was obtained by adopting the following reddening law: it was assumed that the dust distribution is well represented by Seaton’s reddening law, but that part of the dust is in front of the nebula and part is well mixed with the gas. A good fit for all the Paschen and Balmer lines is obtained by adopting a $`C`$(H$`\beta `$) in front of the nebula of 0.31 dex and a $`C`$(H$`\beta `$) from the front to the back of the H II region of 1.70 dex. The effective $`C`$(H$`\beta `$) amounts to 0.92 dex and the effective $`R`$ value to 4.2 (where $`R=A_V/E(BV)`$ ), while Seaton’s law provides an $`R`$ value of 3.2 for material in front of the source.

There are other assumptions that can produce the same or a similar extinction law. For example if all the extinction occurs in front of the nebula but is not uniform across the observing slit it is also possible to obtain effective $`R`$ values higher than 3.2 without changing the grain properties. Also it is possible to have all the extinction between us and the H II region caused by grains with a different size distribution to that implied by Seaton’s law. Consequently the adopted extinction law is empirical and it is not the purpose of this paper to estimate the relative proportion of the different causes that produce an $`R`$ value higher than 3.2.

Since we avoided the brightest stars, the stellar continuum is very weak and the equivalent width of H$`\beta `$ in emission amounts to 734 Å, consequently the correction due to underlying absorption of the H and He lines is negligible and was not taken into account.

Table 2 presents the emission line intensities of 30 Dor. The first and second columns include the observed wavelength in the framework of the solar standard of rest, $`\lambda `$, and the adopted laboratory wavelength, $`\lambda _0`$. The third and fourth columns include the ion and the multiplet number, or the Balmer, B, or Paschen, P, transitions for each line. The fifth and sixth columns include the observed flux relative to H$`\beta `$, $`F(\lambda `$), and the flux corrected for reddening relative to to H$`\beta `$, $`I(\lambda `$). The seventh column includes the fractional error (1$`\sigma `$) in the line intensities.

366 distinct emission lines were measured; of them 43 are unidentified, 31 are blends, and 292 are single lines. By combining the single and blended lines 356 line identifications were obtained; of them 269 are permitted, 85 forbidden and 2 semiforbidden.

Two of the lines presented in Table 2 correspond to a nebular component in emission that originates in our Galaxy and that is not physically associated with 30 Dor. This component has an emission measure in H$`\alpha `$ four orders of magnitude smaller than that of 30 Dor in the observed line of sight. The H$`\beta `$ line of this component was also detected. Both Balmer lines are blueshifted relative to 30 Dor by 252 km sec-1. This nebular Galactic component is not present in other wavelengths and it does not affect the line intensities presented in Table 2, with the possible exception of 4098.24 O II, that might be blended with the Galactic H$`\delta `$ line.

## 4 Physical Conditions

### 4.1 Temperatures and Densities

The temperatures and densities presented in Table 3 were derived from the line intensities presented in Table 2. Most of the determinations were carried out based on the IRAF subroutines. The contribution to the intensities of the $`\lambda \lambda `$ 7319, 7320, 7331, and 7332 \[O II\] lines due to recombination was taken into account based on the following equation:

$$I_R(7319+7320+7331+7332)/I(\mathrm{H}\beta )=9.36(T/10^4)^{0.44}\times \mathrm{O}^{++}/\mathrm{H}^+,$$
(1)
(see Liu et al., 2000). The \[Cl III\] temperature was obtained from the $`I`$(8500)/$`I`$(5518 + 5538) ratio and the computations by Keenan et al. (2000). The Balmer continuum temperature was determined from the following equation:

$$T=368\times (1+0.259y^++3.409y{}_{}{}^{+}+)(\mathrm{Bac}/\mathrm{H11})\mathrm{K},$$
(2)
(see Liu et al., 2001) where Bac/H11 is in Å-1, $`y^+`$ = He+/H+, and $`y^{++}`$ = He++/H+. Figure The Chemical Composition of the 30 Doradus Nebula Derived from VLT Echelle Spectrophotometry111Based on observations collected at the European Southern Observatory, Chile, proposal number ESO 68.C-0149(A). presents the region near the Balmer limit, where the Balmer continuum can be easily appreciated.

The \[Fe III\] density was derived from the $`I`$(4986)/$`I`$(4658) ratio and the computations by Keenan et al. (2001). In Table 4 I compare the observed $`I(\lambda _{nm})`$/$`I`$(4658) ratios for a set of \[Fe III\] lines with the predicted ratios by Keenan et al. (2001) for $`N_e`$ = 316 cm-3 and $`T_e`$ = 10,000 K. I did not list in this table the \[Fe III\] $`\lambda \lambda `$ 4008, 4925 and 4931 lines because probably 4008 has a contribution from an unidentified component, 4925 is blended with a permitted O II line, and 4931 is blended with a forbidden \[O III\] line. From Table 4 it can be seen that the observed and predicted line ratios are in excellent agreement since the differences between them are of the order of their estimated errors.

### 4.2 Temperature variations

To derive the ionic abundance ratios the average temperature, $`T_0`$, and the mean square temperature fluctuation, $`t^2`$, were used. These quantities are given by

$$T_0(N_e,N_i)=\frac{T_e(𝐫)N_e(𝐫)N_i(𝐫)𝑑V}{N_e(𝐫)N_i(𝐫)𝑑V},$$
(3)
and

$$t^2=\frac{(T_eT_0)^2N_eN_i𝑑V}{T_0^2N_eN_i𝑑V},$$
(4)
respectively, where $`N_e`$ and $`N_i`$ are the electron and the ion densities of the observed emission line and $`V`$ is the observed volume (Peimbert, 1967).

The H I and He I lines originate in the high and the low degree of ionization zones and will be represented by $`T_0`$ and $`t^2`$. I will divide the other lines in two groups those that originate mainly in the high degree of ionization zones: C++, O++, Ne++, S++, Cl++, Cl+3, Ar++, and Ar+3, and those that originate mainly in the low degree of ionization zones: N+, O+, S+, Cl+, and Fe++. These two groups will be represented by $`T_{0,H}`$, $`T_{0,L}`$, $`t_H^2`$, and $`t_L^2`$. I will also assume that $`t_H^2=t_L^2t^2`$.

To determine $`T_0`$, $`T_{0,H}`$, $`T_{0,L}`$, and $`t^2`$ we need at least three independent $`T_e`$ determinations that weigh differently the high and low temperature zones and the fraction of the emissivity originating in the high and low ionization zones, that for these observations corresponds to 85% and 15% respectively (see below). For example, it is possible to combine the temperature derived from the ratio of the \[O III\] $`\lambda \lambda `$ 4363, 5007 lines, $`T_{(4363/5007)}`$, the temperature derived from the ratio of the \[N II\] $`\lambda \lambda `$ 5755, 6584 lines, $`T_{(5755/6584)}`$, and the temperature derived from the ratio of the Balmer continuum to $`I(H\beta )`$, $`T_{(\mathrm{Bac}/H\beta )}`$, that are given by

$$T_{(4363/5007)}=T_{0,H}\left[1+\frac{1}{2}\left(\frac{90800}{T_0}3\right)t_H^2\right],$$
(5)
$$T_{(5755/6584)}=T_{0,L}\left[1+\frac{1}{2}\left(\frac{69000}{T_0}3\right)t_L^2\right],$$
(6)
and

$$T_{(\mathrm{Bac}/\mathrm{H}\beta )}=T(\mathrm{Bac})=T_0(11.70t^2).$$
(7)
Similar equations can be derived from the \[O II\], \[S II\], \[S III\], and \[Ar III\] auroral to nebular line ratios. From the abundances of O+ and O++ I estimate that 85% of the H($`\beta `$) emissivity comes from the regions of high degree of ionization and 15% from the regions of low degree of ionization. By combining the forbidden line temperatures, $`T`$(FL) with the $`T_{(\mathrm{Bac}/H\beta )}`$ temperature I obtain the $`t^2`$ value presented in Table 5.

Under the assumption that $`T`$ is constant along the line of sight it is possible to derive the abundances of a given ion from a collisionally excited line of an element $`p`$ times ionized or from a recombination line of the same element $`p1`$ times ionized. In many objects the abundances derived from the recombination lines are higher than those derived from the collisionally excited lines, possibly indicating the presence of temperature variations along the line of sight.

Taking into account that the abundances derived from the recombination lines and the collisionally excited lines are the same it is possible to derive from the ratio of a collisionally excited line to a recombination line a function of $`T_{0,H}`$ and $`t_H^2`$, or $`T_{0,L}`$ and $`t_L^2`$. By combining this relation for lines that originate mainly in regions of high degree of ionization with a temperature determined from the ratio of two collisionally excited lines, like $`T_{(4363/5007)}`$, it is also possible to derive $`T_{0,H}`$ and $`t_H^2`$.

Table 5 presents three $`t^2`$ values derived from the O+(R/C), O++(R/C), and C++(R/C) line intensity ratios. The lines used were: for O+, the three lines of multiplet 1 of O I (7773) together with the \[O II\] $`\lambda 3727`$ lines; for O++, the eight lines of multiplet 1 of O II (4659) together with the \[O III\] $`\lambda 5007`$ line; and, for C++, the C II line $`\lambda `$ 4267 of multiplet 6 together with the $`\lambda \lambda `$ 1906 and 1909 collisionally excited lines of C III. These relations were combined with the temperatures derived from the ratios of forbidden lines for the low degree of ionization zones in the first case, and with the temperatures derived from the ratios of forbidden lines for the high degree of ionization zones in the second and third cases.

## 5 Ionic Chemical Abundances

### 5.1 Helium ionic abundances

To obtain He+/H+ values we need a set of effective recombination coefficients for the He and H lines, the contribution due to collisional excitation to the helium line intensities, and an estimate of the optical depth effects for the helium lines. The recombination coefficients used were those by Storey & Hummer (1995) for H, and those by Smits (1996) and Benjamin, Skillman, & Smits (1999) for He. The collisional contribution was estimated from Sawey & Berrington (1993) and Kingdon & Ferland (1995). The optical depth effects in the triplet lines were estimated from the computations by Benjamin, Skillman, & Smits (Benjamin et al.2002).

Table 6 presents the He+/H+ values obtained from the eleven best observed helium lines for $`t^2=0.033`$. Considering the observational errors, for no value of $`\tau _{3889}`$ it is possible to find agreement in the He+/H+ values among the six triplet helium lines. The four more sensitive lines to $`\tau _{3889}`$ are $`\lambda \lambda `$ 3188, 3889, 4713, and 7065, two of them, $`\lambda \lambda `$ 4713 and 7065, increase in intensity with increasing $`\tau _{3889}`$ and two of them, $`\lambda \lambda `$ 3188 and 3889, decrease.

Two acceptable solutions for $`\tau _{3889}`$ have been obtained. By excluding $`\lambda \lambda `$ 3188 and 3889, $`\tau _{3889}`$ = 4.4 is found; alternatively, by excluding $`\lambda \lambda `$ 4713 and 7065, $`\tau _{3889}`$ = 10.5 is found. This result indicates that the computations for spherical geometry by (Benjamin et al.2002) do not apply to the observed region of 30 Dor. This region is part of a bright shell probably longer along the line of sight than in the plane of the sky, obviously a non spherical object.

A fraction of the $`np2s`$ photons is converted into $`\lambda \lambda `$ 4713, 7065, 4471, and 5876 photons; due to geometrical effects a smaller fraction than that expected, under the sphericaly simetric case, is sent in our direction. I consider that it is a good approximation to assume that the increase in the $`\lambda \lambda `$ 4713, 7065, 4471, and 5876 line intensities corresponds to that predicted by the same $`\tau _{3889}`$ value and that the decrease in the $`\lambda \lambda `$ 3188 and 3889 lines corresponds to another $`\tau _{3889}`$ value. Consequently by averaging the values of the 9 helium lines for $`\tau _{3889}`$ = 4.4 (excluding $`\lambda \lambda `$ 3188 and 3889) a value of He+/H+ = 0.08470 $`\pm 0.00068`$ is obtained.

The five singlet lines, which are not affected by the $`\tau _{3889}`$ effect, yield He+/H+ = 0.08448 $`\pm 0.00099`$ in excellent agreement with the value derived from the nine helium lines and the discussion presented above.

From the observed spectra it is found that the $`I`$(4686)/$`I`$(H$`\beta `$) value is smaller than 3.5 $`\times 10^5`$, which together with the recombination coefficients by Brocklehurst (1971) imply that $`N`$(He++)/$`N`$(H+) is smaller than 2.9 $`\times 10^6`$.

### 5.2 C and O ionic abundances from recombination lines

The C++ abundance was derived from the $`\lambda 4267`$ Å line of C II and the effective recombination coefficients computed by Davey, Storey, & Kisielius (2000) for Case A and $`T`$ = 10,000 K.

The O+ abundance was derived from the $`\lambda \lambda `$ 7771.96, and 7775.40 lines of O I and the effective recombination coefficient for the multiplet computed by Péquignot, Petitjean, & Boisson (1991). The third line of the multiplet, $`\lambda `$ 7774.18, was partially blended with a telluric line in emission, consequently it was not possible to measure its intensity, it was assumed that the three lines of the multiplet are in LS coupling and consequently that $`I(7774.18)=I(7771.96+7775.40)/2`$.

The O++ abundance was derived from the eight lines of multiplet 1 of O II (see Figure The Chemical Composition of the 30 Doradus Nebula Derived from VLT Echelle Spectrophotometry111Based on observations collected at the European Southern Observatory, Chile, proposal number ESO 68.C-0149(A).) together with the effective recombination coefficient for the multiplet computed by Storey (1994) under the assumption of Case B for $`T_e=10,000`$ K and $`N_e`$ = 300 cm-3. The result is almost independent of the case assumed, the difference in the O++/H+ value between Case A and Case B is smaller than 4%. It was found that the O II lines of multiplet 1 are in Case B based on the observed intensities of multiplets 19, 2, and 28 of O II that are strongly case sensitive; Peimbert, Storey, & Torres-Peimbert (1993) also found that the O II lines in the Orion nebula are in Case B. The line intensity ratios within multiplet 1 do not follow the LS coupling predictions. Figure The Chemical Composition of the 30 Doradus Nebula Derived from VLT Echelle Spectrophotometry111Based on observations collected at the European Southern Observatory, Chile, proposal number ESO 68.C-0149(A). provides an excelent visual reference to estimate the quality of the data, it includes 2 lines four orders of magnitude fainter than H$`\beta `$ and also shows that lines separated by 2Å are completely resolved.

### 5.3 Ionic abundances from collisionally excited lines

With the exception of C++/H+ and Fe++/H+ all the other values presented in Table 8 for $`t^2=0.00`$ were derived with the IRAF task abund, using only the low- and medium-ionization zones. The low and medium ionization zones of IRAF correspond to the low and high ionization zones of this paper.

The C++/H+ value for $`t^2=0.00`$ was derived from the collisionally excited lines of C III $`\lambda \lambda `$ 1906 and 1909 by (Dufour, et al.1982) and Garnett et al. (1995). I consider this procedure valid because the O degree of ionization derived here is in excellent agreement with theirs (O++/O equal to 85% and 83% respectively).

The Fe++/H+ value for $`t^2=0.00`$ was derived from the atomic data by Nahar & Pradhan (1996) and Zhang (1996). I did not determine the Fe+/H+ abundances because the observed \[Fe II\] lines are produced by collisions and by non negligible radiative processes difficult to estimate (Rodríguez, 1999).

To derive the abundances for $`t_H^2=t_L^2t^2=0.033`$ I used the abundances for $`t^2=0.00`$ and the formulation for $`t^2>0.00`$ presented by Peimbert & Costero (1969). To derive abundances for other $`t^2`$ values it is possible to interpolate or to extrapolate the values presented in Table 8.

## 6 Total Abundances

Table 9 presents the total abundances of 30 Doradus for $`t_L^2=t_H^2t^2=0.033t^2`$. To derive the total gaseous abundances the set of equations presented below was used, where the ionization correction factors, $`ICF`$’s, correct for the unseen ionization stages.

The total He/H value is given by:

$`{\displaystyle \frac{N(\mathrm{He})}{N(\mathrm{H})}}`$ $`=`$ $`{\displaystyle \frac{N(\mathrm{He}^0)+N(\mathrm{He}^+)+N(\mathrm{He}^{++})}{N(\mathrm{H}^0)+N(\mathrm{H}^+)}},`$ (8)
$`=`$ $`ICF(\mathrm{He}){\displaystyle \frac{N(\mathrm{He}^+)+N(\mathrm{He}^{++})}{N(\mathrm{H}^+)}}.`$
The He++/H+ ratio is completely negligible (see previous section).

In objects of low degree of ionization the presence of neutral helium inside the H II region is important and $`ICF`$(He) becomes larger than 1. To study this problem Vílchez & Pagel (1988, see also ) defined a radiation softness parameter given by

$$\zeta =\frac{N(\mathrm{O}^+)N(\mathrm{S}^{++})}{N(\mathrm{S}^+)N(\mathrm{O}^{++})};$$
(9)
for large values of $`\zeta `$ the amount of neutral helium is significant, while for low values of $`\zeta `$ it is negligible, where the critical value is around 8. From the previous equation and the values in Table 8 it is found that $`\zeta `$ = 2.82 which indicates that the amount of He0 inside the H+ region is negligible.

On the other hand, for ionization bounded objects of very high degree of ionization the amount of H0 inside the He+ Strömgren sphere becomes significant and the $`ICF`$(He) can become smaller than 1. This possibility was firstly mentioned by Shields (1974) and studied extensively by Armour et al. (1999), Viegas, Gruenwald, & Steigman (2000), Ballantyne, Ferland, & Martin (Ballantyne et al.2000) and Sauer & Jedamzik (2002).

According to (Ballantyne et al.2000) for \[O III\]$`\lambda 5007`$/\[O I\]$`\lambda 6300300`$, the $`ICF`$(He) becomes very close to unity; from Table 2 it is found that $`I(5007)/I(6300)=548`$. Consequently I conclude that the amount of H0 inside the observed H II region is negligible and in what follows I will adopt an $`ICF`$(He) = 1.000.

The gaseous abundances for O, N, and Ne were obtained from the following equations (Peimbert & Costero, 1969)

$$\frac{N(\mathrm{O})}{N(\mathrm{H})}=\frac{N(\mathrm{O}^+)+N(\mathrm{O}^{++})}{N(\mathrm{H}^+)},$$
(10)
$$\frac{N(\mathrm{N})}{N(\mathrm{H})}=\left(\frac{N(\mathrm{O}^+)+N(\mathrm{O}^{++})}{N(\mathrm{O}^+)}\right)\frac{N(\mathrm{N}^+)}{N(\mathrm{H}^+)},$$
(11)
and

$$\frac{N(\mathrm{Ne})}{N(\mathrm{H})}=\left(\frac{N(\mathrm{O}^+)+N(\mathrm{O}^{++})}{N(\mathrm{O}^{++})}\right)\frac{N(\mathrm{Ne}^{++})}{N(\mathrm{H}^+)}.$$
(12)
To obtain the total O/H gaseous abundance the O+/H+ and O++/H+ values presented in Tables 7 and 8 were weighed according to their observational errors. To obtain the total O abundances a correction of 0.08 dex was adopted to take into account the fraction of O tied up in dust grains, this fraction was estimated from the Mg/O, Si/O, and Fe/O values derived for the Orion nebula (Esteban et al., 1998).

To obtain the C gaseous abundance the following equation was adopted

$$\frac{N(\mathrm{C})}{N(\mathrm{H})}=ICF(\mathrm{C})\frac{N(\mathrm{C}^{++})}{N(\mathrm{H}^+)},$$
(13)
where the C++/H+ gaseous abundance was obtained by weighing the C++/H+ values presented in Tables 7 and 8 according to their observational errors, the $`ICF`$(C) value was obtained from Garnett et al. (1995). Following Esteban et al. (1998) 0.10 dex were added to the total C/H gaseous value to take into account the fraction of C atoms embedded in dust grains.

The gaseous abundances of S, Cl, Ar and Fe were obtained from the following equations:

$$\frac{N(\mathrm{S})}{N(\mathrm{H})}=ICF(\mathrm{S})\frac{N(\mathrm{S}^+)+N(\mathrm{S}^{++})}{N(\mathrm{H}^+)},$$
(14)
$$\frac{N(\mathrm{Cl})}{N(\mathrm{H})}=\frac{N(\mathrm{Cl}^+)+N(\mathrm{Cl}^{++})+N(\mathrm{Cl}^{3+})}{N(\mathrm{H}^+)},$$
(15)
$$\frac{N(\mathrm{Ar})}{N(\mathrm{H})}=ICF(\mathrm{Ar})\frac{N(\mathrm{Ar}^{++})+N(\mathrm{Ar}^{3+})}{N(\mathrm{H}^+)},$$
(16)
and

$$\frac{N(\mathrm{Fe})}{N(\mathrm{H})}=ICF(\mathrm{Fe})\frac{N(\mathrm{Fe}^{++})}{N(\mathrm{H}^+)}.$$
(17)
The $`ICF`$(S) value was estimated from the models by Garnett (1989). The $`ICF`$(Fe) value was estimated from the models of NGC 346 by Relaño, Peimbert, & Beckman (Relaño et al.2002) and amounts to $`8\pm 2`$. The $`ICF`$(Ar) includes the Ar+/H+ contribution and according to Liu et al. (2000) can be approximated by

$$ICF(\mathrm{Ar})=\left(1\frac{N(\mathrm{O}^+)}{N(\mathrm{O})}\right)^1.$$
(18)
## 7 Discussion and Conclusions

### 7.1  The H II regions and the solar abundances

Table 10 presents the abundances of four very well observed H II regions: Orion and M17 in the Galaxy, NGC 346 in the SMC, and 30 Doradus in the LMC. The Orion, M17, and 30 Doradus values include additions to the C and O gaseous values of 0.1 dex and 0.08 dex respectively to take into account the fraction of these atoms embedded in dust. For the same reason the NGC 346 C and O values include additions of 0.05 dex and 0.04 dex respectively (, Relaño et al.2002). The H II region abundances have been obtained adopting values of $`t^2`$ larger than 0.00. Further arguments in favor of $`t^2>0.00`$ have been presented elsewhere (Peimbert, 2002; Peimbert & Peimbert, 2002a, b).

In addition Table 10 presents also the solar photospheric values for C, N, O, Ne, and Ar, and the solar abundances derived from meteoritic data for S, Cl, and Fe. For the solar initial helium abundance the $`Y_0`$ by Christensen-Dalsgaard (1998) was adopted, and not the photospheric one because, apparently, it has been affected by settling.

The H II region values for Ne/O, S/O, and Ar/O are in excellent agreement with the solar values which implies that the production of these elements is primary and due to massive stars, and that the assumptions involved in the two types of abundance determinations are sound.

The H II regions S/O and Cl/O abundances are in better agreement with the solar meteoritic abundances than with the photospheric ones, the solar photospheric values are S/O = -1.38 $`\pm 0.11`$ dex, and Cl/H = -3.21 $`\pm 0.3`$ dex.

The 30 Doradus C/O value is intermediate between that of NGC 346 and those of Orion, M17 and the Sun. The differences are significant and imply that even if C is of primary origin part of it is due to intermediate mass stars and part is due to massive stars, and that the C yield increases with the O/H ratio (Garnett et al., 1995, 1999; Carigi, 2002).

The 30 Doradus N/O value is intermediate between that of NGC 346 and those of Orion, M17 and the Sun. The differences are significant and imply that part of the N is of primary origin and part of secondary origin (see Henry, Edmunds & Koppen, 2000, and references therein).

The accuracies of the He/H abundances of 30 Doradus, NGC 346, and M17 are higher than that of the Orion nebula because the first three objects have $`ICFs`$(He) = 1.00, while the $`ICF`$(He) for the Orion nebula is larger than 1.00 and it is not well determined. Similarly the values for 30 Doradus, NGC 346, and M17 are more accurate than the solar one because they are based on direct determinations, while the solar value is obtained from models that depend on the helium abundance in a more complex way.

### 7.2 The $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ and $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }Z`$ ratios and chemical evolution

The determination of the $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ and $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }Z`$ ratios is crucial for the determination of $`Y_p`$, and for constraining the models of galactic chemical evolution. The abundance determinations of 30 Doradus are based on emission line intensities of high quality, take into account the temperature structure of the nebula, and include a very accurate He/H value because the degree of ionization of 30 Doradus is relatively high implying an $`ICF`$(He) very close to unity.

To determine the hydrogen, helium, heavy elements, and oxygen abundance by mass, $`X`$, $`Y`$, $`Z`$, and $`O`$, presented in Table 10 I proceeded as follows: for the Sun I adopted the initial helium abundance by mass ($`Y_0`$) of Christensen-Dalsgaard (1998), the $`Z/X`$ value derived from the C/H, N/H, O/H, Ne/H, and Ar/H values presented in Table 10, and $`0.56\times O/X`$ for the ratio of the rest of the heavy elements to hydrogen, value obtained from the meteoritic abundances by Grevesse & Sauval (1998); for the H II regions I adopted the He/H value determined by the different observers, the C/H, N/H, O/H, Ne/H, and Ar/H values presented in Table 10, and $`0.56\times O/X`$ for the rest of the heavy elements.

From Table 11 it can be seen that the spread among the $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ and $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }Z`$ values derived from the H II regions is smaller for the $`t^2>0.00`$ results than for the $`t^2=0.00`$ results, which is consistent with the idea that the $`t^2>0.00`$ values are better. Moreover the theoretical computations for the chemical evolution of irregular galaxies and for the solar vicinity predict values of $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ in the 2.9 to 4.6 range with a representative value of $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ = 3.5 $`\pm 0.9`$ (see Carigi et al., 1995; Carigi, Colín, & Peimbert, 1999; Carigi, 2000; Chiappini, Matteucci, & Gratton, 1997; Peimbert et al., Peimbert, Peimbert, & Ruiz2000), again in better agreement with the results from H II regions under the assumption that $`t^2>0.00`$.

### 7.3 The primordial helium abundance

To determine the $`Y_p`$ value from 30 Doradus it is necessary to estimate the fraction of helium present in the interstellar medium produced by galactic chemical evolution. For this purpose it was assumed that:

$$Y_p=Y(30\mathrm{Dor})O(30\mathrm{Dor})\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }O}.$$
(19)
As in section 7.2 a $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ = 3.5 $`\pm 0.9`$ will be adopted, which together with the $`Y`$(30 Dor) and $`Z`$(30 Dor) values of Table 10 yield $`Y_p`$ = 0.2345 $`\pm `$ 0.0047, where most of the error comes from the adopted $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ value. This $`Y_p`$ value is in excellent agreement with the value derived from NGC 346 (Peimbert et al., Peimbert, Peimbert, & Ruiz2000). This agreement is due to the similarity between the adopted $`\mathrm{\Delta }Y`$/$`\mathrm{\Delta }O`$ value and that derived from 30 Doradus and NGC 346.

I am grateful to María Teresa Ruiz and Manuel Peimbert for carrying out the observations presented in this paper and for fruitful discussions, to César Esteban for his collaboration in the initial stages of this project, and to Anabel Arrieta for her assistance in the reduction process. It is also a pleasure to thank the Departamento de Astronomía de la Universidad de Chile for its hospitality during a visit where part of this work was made.


text:
# Expansion of an interacting Fermi gas

## Abstract

We study the expansion of a dilute ultracold sample of fermions initially trapped in a anisotropic harmonic trap. The expansion of the cloud provides valuable information about the state of the system and the role of interactions. In particular the time evolution of the deformation of the expanding cloud behaves quite differently depending on whether the system is in the normal or in the superfluid phase. For the superfluid phase, we predict an inversion of the deformation of the sample, similarly to what happens with Bose-Einstein condensates. Viceversa, in the normal phase, the inversion of the aspect ratio is never achieved, if the mean field interaction is attractive and collisions are negligible.

Since the first experiments on trapped Bose-Einstein condensates the imaging of the expanding cloud, following the sudden switching off of the confining potential, has provided crucial information on the novel features exhibited by atomic gases in conditions of quantum degeneracy. These include, in particular, the bimodal structure of an expanding Bose gas at finite temperature and the anisotropy of the asymptotic profile of the condensate . In the Thomas-Fermi limit, where the Gross-Pitaevskii equations coincide with the hydrodynamic theory of superfluids, the anisotropy of the expanded gas reflects the anisotropy of the pressure force which is stronger in the direction of tighter confinement . The predictions of the hydrodynamic equations and of the consequent scaling behaviour exhibited by Bose-Einstein condensates during the expansion have been investigated by several authors (,,), providing excellent agreement with experiments  and pointing out the difference with respect to the expansion of a non condensed gas. In the latter case, in the collisionless regime, the density profile approaches an isotropic shape, independent of the initial deformation of the gas.

In this letter we study the problem of the expansion of an ultracold sample of fermions initially trapped in an anisotropic harmonic trap. We will show that also in the case of fermions the expansion of the gas provides valuable information about the state of the system and the role of interactions. We will consider a gas of atoms interacting with attractive forces. This is a natural requirement for the realization of Cooper-pairs and hence for the achievement of the superfluid phase . Such interactions are naturally present in some fermionc species like 6Li and can otherwise be obtained by changing the scattering length profiting of the presence of a Feshbach resonance.

The description of the expansion of a cold fermionic gas in the normal and superfluid phase requires different theoretical approaches. For the normal phase we use the formalism of the Landau-Vlasov equations, while in the superfluid phase we study the expansion using the hydrodynamic theory of superfluids.

We consider the case of two different fermionic states, hereafter called 1 and 2, initially confined in a harmonic trap. We assume that the two species are present in the same amount and feel the same trapping potential, so that the densities of the two species are equal at equilibrium as well as during the expansion: $`n_1(𝐫,t)=n_2(𝐫,t)=n(𝐫,t)/2`$. The trapping potential will be chosen of harmonic type

$`V_{ho}={\displaystyle \frac{1}{2}}m(\omega _{}^2x^2+\omega _{}^2y^2+\omega _z^2z^2),`$ (1)
describing a cilyndrically simmetric trap with deformation $`\lambda =\omega _z/\omega _{}`$. The interaction between the two fermionic species is fixed by the coupling constant $`g=4\pi \mathrm{}^2a/m`$, where $`a`$ is the $`s`$-wave scattering length.

In this paper we will use the equation of state

$`\mu _\mathrm{}e(n)={\displaystyle \frac{\mathrm{}^2}{2m}}(3\pi ^2n)^{2/3}+{\displaystyle \frac{1}{2}}gn,`$ (2)
to describe the uniform phase of the gas where the first term is the kinetic energy evaluated at zero temperaure while the second one is the interaction energy evaluated in the mean field approximation. Equation (2) neglects the effects of correlations which become important for large values of the scattering length and affect in a different way the equation of state of the normal and superfluid phase. The formalism developed in this paper can be easily generalized to include a more accurate description of the equation of state. It is however worth pointing out that, even using the same equation of state for the normal and superfluid phases, the expansion of the gas behaves quite differently in the two cases being described by different kinetic equations. In the presence of the external potential (1), the equilibrium condition in the local density approximation is determined by the equation

$$\mu _\mathrm{}e(n)+V_{ho}(𝐫)=\mu $$
(3)
where $`\mu `$ is the chemical potential of the sample fixed by the normalization condition.

The relevant parameter characterizing the interaction in the fermionic system is the ratio

$$\chi =\frac{E_{int}}{E_{ho}}$$
(4)
between the interaction energy

$`E_{int}={\displaystyle \frac{g}{4}}{\displaystyle n^2(𝐫)d^3r}`$ (5)
and the oscillator energy

$`E_{ho}`$ $`=`$ $`{\displaystyle V_{ho}(𝐫)n(𝐫)d^3r}.`$ (6)
In the perturbative regime the integrals (5,6) can be evaluated using the non interacting density profile which, in Thomas-Femi approximation (3), takes the simple form

$$n(𝐫)=\frac{1}{3\pi ^2}\left(\frac{2m}{\mathrm{}^2}\right)^{3/2}\left[\mu V_{ho}(𝐫)\right]^{3/2}.$$
(7)
After integration of (5) and (6) one finds

$`{\displaystyle \frac{E_{int}}{E_{ho}}}=0.5{\displaystyle \frac{N^{1/6}a}{a_{ho}}}=0.3ak_F,`$ (8)
where $`k_F=(3\pi ^2n(0))^{1/3}`$ is the Fermi momentum evaluated at the central value of the density. In order to go beyond the perturbative regime, one determines numerically the ground state density by minimizing the energy of the system

$`E`$ $`=`$ $`{\displaystyle \frac{3(3\pi ^2)^{2/3}}{5}}{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle n^{5/3}(𝐫)d^3r}+`$ (9)
$`+`$ $`{\displaystyle V_{ho}(𝐫)n(𝐫)d^3r}+{\displaystyle \frac{g}{4}}{\displaystyle n^2(𝐫)d^3r},`$ (10)
One easily finds that the equilibrium value of the ratio $`E_{int}/E_{ho}`$ depends on the dimensionless combination $`aN^{1/6}/a_{ho}`$ also in the non perturbative regime , as reported in Fig.1. Using (10) one predicts that the compressibility of the gas becomes negative in the center of the trap if $`|a|N^{1/6}/a_{ho}>0.61`$. For large values of $`|a|`$ the resulting predictions should be however taken with care since the functional (10) ignores correlation effects beyond mean field.

Let us now discuss the expansion of a fermionic sample trapped in an elongated harmonic trap ($`\lambda <1`$). We describe first the expansion of the normal fluid and afterwards the one of the superfluid. In the ideal case, using the semiclassical description, one finds that the ratio of the square radii evolves according to the classical law

$$\frac{r_{}^2}{z^2}=\frac{1+\omega _{}^2t^2}{1+\omega _z^2t^2}\frac{\omega _z^2}{\omega _{}^2},$$
(11)
The ratio (11) approaches unity for large times, reflecting the isotropy of the momentum distribution. This result ignores the effects of collisions which are however expected to play a minor role at low temperature, due to Pauli blocking, unless the scattering length is very large.

In order to take into account the effects of the interactions, we consider the mean field description based on the Landau-Vlasov equation

$`{\displaystyle \frac{f}{t}}+𝐯{\displaystyle \frac{f}{𝐫}}{\displaystyle \frac{V_{ho}}{𝐫}}{\displaystyle \frac{f}{𝐯}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{g}{m}}{\displaystyle \frac{n}{𝐫}}{\displaystyle \frac{f}{𝐯}}=0,`$ (12)
where $`f(𝐫,𝐯,t)`$ is the distribution function, $`n=fd^3v`$ is the atomic density. Eq.(12) describes the dynamics of a normal weakly interacting gas in the collisionless regime.

An approximate solution of Eq.(12), is obtained by making a scaling ansatz for the distribution function

$`f(𝐫,𝐯,t)=f_0(\stackrel{~}{𝐫}(t),\stackrel{~}{𝐯}(t)),`$ (13)
where $`f_0`$ is the equilibrium distribution, $`\stackrel{~}{r}_i(t)=r_i/b_i`$ and $`\stackrel{~}{v}_i(t)=b_iv_i\dot{b}_ir_i`$. Under the scaling assumption the velocity field $`𝐮(𝐫,t)=𝐯fd^3v/n`$ takes the simple form $`u_i=\dot{b}_ir_i/b_i`$. This ansatz has been recently used by Guéry-Odelin  to investigate the effect of the interaction on the collective oscillation of a classical gas in the collisionless regime.

The equations for the scaling parameters $`b_i`$ can be obtained by multiplying (12) by $`\stackrel{~}{r}_i`$ and $`\stackrel{~}{v}_i`$ and integrating in phase space. Making use of the equilibrium properties of the distribution function, after some straigthforward algebra one finds

$`\ddot{b}_i+\omega _i^2b_i{\displaystyle \frac{\omega _i^2}{b_i^3}}+{\displaystyle \frac{3}{2}}\chi \omega _i^2\left({\displaystyle \frac{1}{b_i^3}}{\displaystyle \frac{1}{b_i_jb_j}}\right)=0,`$ (14)
where $`\chi `$ is the ratio (4) evaluated at equilibrium. The second term in (14) describes the restoring force of the external oscillator potential, the third one originates from the kinetic energy, while the last term, linear in $`\chi `$, accounts for the effects of mean field interaction.

An immediate application of Eq.(14) concerns the study of the oscillations of the gas. By linearizing the equations around equilibrium ($`b_i=1`$) one finds, in the presence of isotropic harmonic trapping ($`\omega _{}=\omega _z=\omega _0`$), the result:

$`\omega _M=2\omega _0\sqrt{1+3\chi /8},\omega _Q=2\omega _0\sqrt{13\chi /4},`$ (15)
for the frequencies of the monopole and quadrupole oscillations which coincide with the results already derived in  using a sum-rule approach.

The equations describing the expansion are obtained by suddenly removing the second term of Eq.(14), originating from the trapping potential. In the study of the expansion we are interested in the case of anisotropic trapping. In particular we will consider the case of cigar shaped traps. For high deformations ($`\lambda =\omega _z/\omega _{}1`$) Eq.(14) yields the asymptotic result

$`b_z^2`$ $``$ $`\omega _z^2(1{\displaystyle \frac{3}{2}}\chi )t^2,`$ (16)
$`b_{}^2`$ $``$ $`\omega _{}^2t^2,`$ (17)
showing that the aspect ratio

$`{\displaystyle \frac{R_{}(t)}{R_z(t)}}{\displaystyle \frac{1}{\sqrt{13\chi /2}}}`$ (18)
approaches a value smaller than 1 if the interaction is attractive ($`\chi <0`$). In Eq.(18), $`R_{}`$ and $`R_z`$ are the radii where the atomic density vanishes (Thomas-Fermi radii).

The results of the numeric integration of the equations of motion (14) are reported in Fig.2 and 3 as a function of time for the choices $`\chi =0`$ and $`\chi =0.4`$.

We address now the problem of the expansion of a superfluid Fermi gas. As already anticipated we will make use of the hydrodynamic equations of superfluids. Those equations have been already used to describe the collective oscillations of a superluid trapped Fermi gas  including its rotational behaviour . The hydrodynamic equations are applicable if the healing length is much smaller than the size of the sample, which implies that the energy gap should be larger than the oscillator energies $`\mathrm{}\omega _z`$, $`\mathrm{}\omega _{}`$. This non trivial condition implies that the whole system behaves like superfluid. Furthermore the hydrodynamic equations are applicable up to excitation energies of the order of the energy gap. In the problem of the expansion it is crucial that the system remains superfluid in the first instants when the hydrodynamic forces provide the relevant acceleration to the expanding atoms. One expects that this condition be satisfied if the initial temperature is small enough.

The hydrodynamic description is based on the equation of continuity

$`{\displaystyle \frac{}{t}}n+(n𝐮)=0`$ (19)
and on the Euler equation

$`m{\displaystyle \frac{}{t}}𝐮+\left(\mu _\mathrm{}e(n)+V_{ho}(𝐫)+{\displaystyle \frac{1}{2}}mu^2\right)=0`$ (20)
where $`\mu _\mathrm{}e(n)`$ is the chemical potential of a uniform gas calculated at the density $`n`$ and $`𝐮`$ is the velocity field.

If the equation of state is a power law ($`\mu _\mathrm{}en^\gamma `$) these equations admit the simple scaling solution

$`n(r_i,t)`$ $`=`$ $`{\displaystyle \frac{1}{_jb_j}}n_0\left({\displaystyle \frac{r_i}{b_i}}\right),`$ (21)
$`u_i(r_i,t)`$ $`=`$ $`{\displaystyle \frac{\dot{b}_i}{b_i}}r_i,`$ (22)
and the Thomas-Fermi radii evolve according to the law $`R_i(t)=R_i(0)b_i(t)`$. In this case, it is immediate to show that, during the expansion, the scaling parameters obey the coupled differential equations

$$\ddot{b}_i=\frac{\omega _i^2}{(b_xb_yb_z)^\gamma b_i}.$$
(23)
which reduce to

$$\ddot{b}_z=\frac{\omega _z^2}{b_{}^{2\gamma }},\mathrm{and}\ddot{b}_{}=\frac{\omega _{}^2}{b_{}^{(2\gamma +1)}}.$$
(24)
for highly elongated configurations ($`\lambda 1`$). For $`\gamma =1`$ (corresponding to a Bose-Einstein condensed gas) the equation for the radial motion is integrable analytically and one finds the result $`b_{}(t)=(1+\omega _{}^2t^2)^{1/2}`$ .

To describe the expansion of superfluid Fermi gas we will use the same equation of state (2) as for the normal phase. The case of a very dilute gas is also described by a power law with $`\gamma =2/3`$ (first term in (2)). For $`\lambda =0.1`$ and $`0.3`$ the solution is given by the full upper line in Fig.2 and 3 respectively which show that the deformation of the trap is inverted in time and the aspect ratio $`R_{}/R_z`$ reaches asymptotically a value significantly larger than 1 . Superfluidity has hence the effect of distributing the release energy in a strongly asymmetric way along the axial and radial directions. It is worth noticing that the same scaling equations (23), with $`\gamma =2/3`$, are obtained for a classic gas in the collisional regime .

In the more general case (2), a useful approximation to the solution of the hydrodynamic equations, based on the scaling ansatz (21,22), is obtained by multiplying the Euler’s equation (20) by $`r_in(𝐫)`$ and integrating over the spatial coordinates. Using the equation of state (2), one finally obtains the following set of differential equations

$`\ddot{b}_i+\omega _i^2b_i{\displaystyle \frac{\omega _i^2}{b_i}}{\displaystyle \frac{1}{(_ib_i)^{2/3}}}+`$ (25)
$`+{\displaystyle \frac{3}{2}}\chi {\displaystyle \frac{\omega _i^2}{b_i}}\left({\displaystyle \frac{1}{(_ib_i)^{2/3}}}{\displaystyle \frac{1}{_ib_i}}\right)=0,`$ (26)
with $`\chi `$ defined by Eq.(4). Equations (26) differ from the anlalogous equations (14) holding in the normal phase. By linearizing Eqs.(26) around $`b_i=1`$ one gets, in the case of a spherical trap, the result $`\omega _Q=\sqrt{2}\omega _0`$ for the quadrupole frequency , which, contrary to (15), is independent of the interaction term in $`\chi `$.

The predictions of Eqs.(26) for the expansion of the gas are reported in Figs.2,3 and show that the inclusion of the interaction term in the equation of state affects the expansion of the superfluid only in a minor way.

In conclusion we have shown that the expansion of a superfluid Fermi gas, being governed by the equations of hydrodynamics, differs in a crucial way from the one of a normal gas in the collisionless regime. From a theoretical point of view several questions remain to be investigated: among them, the effect of large scattering lengths  on the equation of state and the role of collisions which, under certain conditions, might give rise to a hydrodynamic regime, and hence to anisotropic expansion, also in the normal phase. Finally one should develop the formalism at finite temperature where both the normal and superfluid components are present. The resulting bimodal structure in the expanding cloud is expected to be affected by the transfer of atoms from the superfluid to the normal component during the first stage of the expansion.

This research is supported by the Ministero dell’Istru-
zione, dell’Università e della Ricerca (MIUR). We acknowledge R. Combescot, G. Shlyapnikov and J. E. Thomas for useful discussions.