text:
# When only two thirds of the entanglement can be distilled

## Abstract

We provide an example of distillable bipartite mixed state such that, even in the asymptotic limit, more pure-state entanglement is required to create it than can be distilled from it. Thus, we show that the irreversibility in the processes of formation and distillation of bipartite states, recently proved in \[G. Vidal, J.I. Cirac, Phys. Rev. Lett. 86, (2001) 5803-5806\], is not limited to bound-entangled states.

Distillation is one of the basic concepts in entanglement theory. As shown in the pioneering works on entanglement transformations , it is possible to use local operations and classical communication (LOCC) to convert, in the asymptotic limit ($`N\mathrm{}`$), $`N`$ copies of some bipartite mixed state $`\rho `$ into $`M`$ of copies of some reference pure state, the maximally entangled state

$$|\mathrm{\Phi }\frac{1}{\sqrt{2}}(|00+|11)$$
(1)
of a two-qubit system, which is said to contain one ebit (entangled bit). Moreover, the distillable entanglement $`E_D(\rho )`$, defined as the maximal achievable yield $`M/N`$, was shown to be often finite. This is a remarkable result with important implications in quantum information theory. It says, for instance, that a noisy channel can be used to establish perfect quantum communication between two distant parties, if these are allowed to perform LOCC. Indeed, the imperfect channel can be used to create many copies of some mixed entangled state $`\rho `$, which can then be purified into fewer copies of $`|\mathrm{\Phi }`$ and subsequently used to achieve perfect quantum communication through teleportation .

A notion dual to distillation is that of preparation of $`\rho `$ using pure-state entanglement and LOCC . Now $`M`$ copies of $`|\mathrm{\Phi }`$ are transformed into $`N`$ copies of $`\rho `$. The entanglement cost $`E_C(\rho )`$  (asymptotic version of the entanglement of formation $`E_F(\rho )`$ ) is defined as the minimal ratio $`M/N`$ asymptotically achievable by LOCC. $`E_C(\rho )`$ quantifies the amount of pure-state entanglement required to create a copy of $`\rho `$, in the above asymptotic sense.

Notice that the processes of formation and distillation can be concatenated into a cycle. Starting from $`NE_C(\rho )`$ copies of $`|\mathrm{\Phi }`$, two distant parties can use LOCC to prepare $`N`$ copies of $`\rho `$; and the $`N`$ copies of $`\rho `$ can be subsequently distilled back into $`NE_D(\rho )`$ copies of $`|\mathrm{\Phi }`$,

$$|\mathrm{\Phi }^{NE_C(\rho )}\rho ^N|\mathrm{\Phi }^{NE_D(\rho )}.$$
(2)
Already in the early contributions it was suggested that maybe sometimes this cycle can not be closed completely, in that perhaps not all the initial pure-state entanglement used in the preparation process can be recovered through distillation. That is, maybe an irreversible loss of quantum correlations takes place during the mixing of pure-state entanglement into $`\rho ^N`$ and, accordingly, the distillable entanglement $`E_D(\rho )`$ is smaller than the entanglement cost $`E_D(\rho )`$.

Very recently this phenomenon has been proved to indeed occur . In particular, it has been shown that some undistillable bipartite state $`\rho _b`$ —i.e. with $`E_D(\rho _b)=0`$ ebits— has non-vanishing entanglement cost. Notably, the irreversibility observed in the asymptotic preparation and distillation of $`\rho _b`$ remains even when LOCC are supplemented with loaned pure-state entanglement, to be returned after the manipulation, in the so-called catalytic LOCC setting.

The results in  still leave, however, an important question open. One could associate the irreversibility demonstrated there to the fact that the state $`\rho _b`$ is bound entangled, that is, to the remarkable property that no pure-state entanglement at all can be distilled from it . It could well be the case that the gap observed between $`E_C`$ and $`E_D`$ is just a characteristic feature of some bound entangled states, whereas $`E_C=E_D`$ always holds for distillable states. After all, this is the case for bipartite pure states  and some simple cases of mixed state , which exhaust all the cases where $`E_C`$ and $`E_D`$ have been computed.

In this paper we will present an example of bipartite mixed state $`\sigma `$ that can be distilled, that is $`E_D(\sigma )>0`$, and such that $`E_C(\sigma )>E_D(\sigma )`$. We extend thereby the irreversibility result of  to the case of distillable states. In particular, the extension also holds for catalytic LOCC transformations.

A widely recognized, major problem concerning the study of mixed-state entanglement is that it is very difficult to compute the asymptotic measures $`E_C`$ and $`E_D`$. Here, however, we are not interested in the actual values of $`E_C(\sigma )`$ and $`E_D(\sigma )`$. For the present purposes it is sufficient to show that $`\sigma `$ can be distilled, and to bound $`E_D(\sigma )`$ and $`E_C(\sigma )`$ tight enough from above and from below, respectively, so that the bounds already imply a gap between the two quantities. We start by collecting an amalgam of useful facts.

($`i`$) A sufficient condition for a mixed state $`\rho `$ to be distillable is that a projector $`P`$ into a $`𝒞^2𝒞^2`$ subspace (that is, a subspace which is the tensor product of two-dimensional subspaces for each of the two separated parts of the composite system) exists such that the projection $`P\rho P^{}`$ is still entangled , that is, such that the partial transposition of $`P\rho P^{}`$ has a negative eigenvalue.

($`ii`$) The logarithmic negativity $`E_N(\rho )\mathrm{log}_2(1+2𝒩(\rho ))`$ , where $`𝒩(\rho )`$ is the absolute value of the sum of negative eigenvalues of partial transposition of $`\rho `$, is an upper bound to the distillable entanglement. Werner’s bound reads ,

$$E_D(\rho )E_N(\rho ).$$
(3)
In addition, $`E_N`$ is an additive function,

$$E_N(\rho _1\rho _2)=E_N(\rho _1)+E_N(\rho _2),$$
(4)
which only vanishes for states with positive partial transposition (PPT states),

$$E_N(\rho _{PPT})=0.$$
(5)
Finally, and very important to us, $`E_N(\rho )`$ is a continuous function of $`\rho `$.

($`iii`$) The entanglement of formation $`E_F`$  of $`\rho `$ is bounded below by

$$E_F(\rho )\mathrm{log}_2\alpha ,$$
(6)
where $`\alpha `$ is the maximal overlap of a product state $`|ab`$ with the projector $`\mathrm{\Pi }`$ onto the support of $`\rho `$,

$$\alpha \underset{|ab}{\mathrm{max}}ab|\mathrm{\Pi }|ab.$$
(7)
Accordingly, the entanglement cost $`E_C(\rho )`$ is bounded below by

$$E_C(\rho )log_2\beta ,$$
(8)
if for all $`N`$ the maximal overlap of a normalized product vector $`|a_Nb_N`$ with the $`N`$-fold tensor product of $`\mathrm{\Pi }`$ is at most $`\beta ^N`$,

$$\underset{|a_Nb_N}{\mathrm{max}}a_Nb_N|\mathrm{\Pi }^N|a_Nb_N\beta ^N.$$
(9)
($`iv`$) The four-dimensional subspace $`V𝒞^3𝒞^3`$ orthogonal to the five product vectors

$`|0`$ $``$ $`(|0+|1),`$ (10)
$`(|0+|1)`$ $``$ $`|2,`$ (11)
$`|2`$ $``$ $`(|1+|2),`$ (12)
$`(|1+|2)`$ $``$ $`|0,`$ (13)
$`(|0|1+|2)`$ $``$ $`(|0|1+|2),`$ (14)
does not contain product vectors . The projector $`\mathrm{\Pi }_b`$ onto $`V`$ satisfies: ($`a`$) it has a PPT , ($`b`$) it fulfills Eq. (9) with $`\beta <0.99`$ .

We introduce now a one-parameter family of states

$$\sigma (p)(1p)\rho _b+p|\psi \psi |,$$
(15)
where $`\rho _b\mathrm{\Pi }_b/4`$ is the PPT bound entangled state introduced in  in the context of the so-called non-extendible product basis, and used in  to prove irreversibility of asymptotic manipulations, and

$$|\psi \frac{1}{\sqrt{6}}(|00|012|11)$$
(16)
is an entangled pure state that is orthogonal to all product states of Eq. (14), that is $`|\psi V`$. For $`p=0`$ we recover $`\rho _b`$, for which we know that $`E_C(\rho _b)>\mathrm{log}_20.99>D(\rho _b)=0`$. In what follows we will use facts ($`i`$)-($`iv`$) and perturbation theory to show that for $`p>0`$ we encounter states $`\sigma (p)`$ which can be distilled, and with $`E_C(\sigma _p)>E_D(\sigma _p)`$.

The family of states $`\sigma (p)`$ in Eq. (15) has been carefully chosen to fulfill two important properties. First, $`\sigma (p)`$ is supported on $`V`$, since $`V`$ is the support of $`\rho _b`$ and also $`|\psi `$ is supported in $`V`$. Using ($`iii`$) and ($`iv.b`$) this means that for any $`p[0,1]`$ we have a constant lower bound for $`E_C`$.

Property 1: The entanglement cost of $`\sigma (p)`$, $`p[0,1]`$, is bounded below by

$$E_C(\sigma (p))>log_20.99=0.015\text{ebits}.$$
(17)
Let $`\rho _b^{T_A}`$ denote the partial transposition of $`\rho _b`$, and $`P`$ the rank-four, product projector $`(|00|+|11|)(|00|+|11|)`$. Notice that, by construction, $`\rho _b^{T_A}=\rho _b0`$, $`(P\rho _bP^{})^{T_A}=P\rho _bP^{}`$, and $`P|\psi =|\psi `$. The second important feature of $`\sigma (p)`$ is that, for any $`p>0`$ the partial transposition of the projection $`P\sigma (p)P^{}`$,

$$(P\sigma (p)P^{})^{T_A}=(1p)P\rho _bP^{}+p|\psi \psi |^{T_A},$$
(18)
has a negative eigenvalue $`n`$. Therefore, because of fact ($`i`$), the corresponding state $`\sigma (p)`$ can be distilled.

Property 2: For $`p(0,1]`$, the state $`\sigma (\rho )`$ can be distilled, that is

$$E_D(\sigma (p))>0\text{ ebits.}$$
(19)
For instance, for $`p=0.015`$, $`|n|=2.7\times 10^4`$ (see also Fig. (1)).

Property 2 has been achieved by selecting a projector $`P`$ such that $`(P\rho _bP^{})^{T_A}`$ has only rank three, and thus one vanishing eigenvalue, whereas $`|\psi V`$ has been chosen so that $`(P|\psi \psi |P^{})^{T_A}=|\psi \psi |^{T_A}`$, that is, so that the negative eigenvalue of $`|\psi \psi |^{T_A}`$ entirely contributes to (18). We can use perturbation theory to check what is the effect of such choices.

Let $`M`$ and $`N`$ be finite dimensional hermitian operators, $`M=_{i=0}^lm_i|m_im_i|`$ the spectral decomposition of $`M`$, with $`m_i`$ its decreasingly ordered eigenvalues and $`m_0m_1`$, and let $`ϵ`$ be a small parameter. Then the lowest eigenvalue of $`M+ϵN`$ is, as given by perturbation theory ,

$$m_0+ϵ0|N|0+ϵ^2\underset{i=1}{\overset{l}{}}\frac{|m_0|M|m_0|^2}{m_0m_i}+𝒪(ϵ^3).$$
(20)
Making the proper identifications we realize that the negative eigenvalue $`n(p)`$ of the operator in Eq. (18) is

$$n=|k|p^2+𝒪(p^3),$$
(21)
where $`|k|>0`$ is of the order of $`1`$ and the zero and first order contributions vanish due, respectively, to the fact that the smallest eigenvalue of Eq. (18) vanishes, and to the fact that the corresponding eigenvector, $`|\tau |0(|0+|1)/\sqrt{2}`$, fulfills

$$\tau |(|\psi \psi |^{T_A})|\tau =\tau |\psi \psi |\tau =0.$$
(22)
Finally, for $`p1`$ such that contributions $`𝒪(p^3)`$ may become important, numerical calculations show that $`|n|`$ grows monotonically with $`p`$ (see Fig. (1)).

Summarizing, so far we have learn that $`\sigma (p)`$ can be distilled for any $`p>0`$, while the entanglement cost is bounded below by Eq. (17). In order to complete the result we need to prove that the distillable entanglement of $`\sigma (p)`$ is, in some regime of $`p(0,1]`$, smaller than the lower bound (17). This would already follow from the above if $`E_D(\sigma (p))`$ were a continuous function of $`p`$. For $`p=0`$ we have the bound entangled state $`\rho _b`$, that is, $`E_D(\sigma (0))=0`$ ebits, whereas at the other extreme, $`p=1`$, we have the pure entangled state $`|\psi `$, whose distillable entanglement $`E_D`$ (and entanglement cost $`E_C`$) can be easily computed and reads $`E_D(\sigma (1))=0.55`$ ebits. But, unfortunately, we can not base our argument in the continuity of $`E_D(\sigma (p))`$ as a function of $`p`$, to conclude that an intermediate $`p`$ must exist such that the distillable entanglement is non-zero and still below the bound (17). Whereas it may well be that $`E_D(\rho )`$ is a continuous function of $`\rho `$, this has not been proved. Notice that a plausible objection to continuity relies on the fact that $`E_D(\rho )`$ is actually a function of $`\rho ^N`$ in the large $`N`$ limit. Therefore, a small perturbation of $`\rho `$, which produces a large deviation in $`\rho ^N`$, may imply a discontinuous change in $`E_D(\rho )`$.

Nevertheless, following fact ($`ii`$), the logarithmic negativity $`E_N(\sigma (p))`$ is a continuous upper bound for $`E_D(\sigma (p))`$ (see Fig. (1)). A direct calculation of $`E_N(\sigma (p))`$ finally proves the irreversibility of the preparation-distillation cycle for distillable bipartite states. In particular, for $`p=0.0015`$ we have

$$E_D(\sigma (0.0015))<E_N(\sigma (0.0015))=0.012\text{ ebits}.$$
(23)
Thus, $`\sigma (0.0015)`$ is an example of distillable state with a finite gap $`E_CE_D>0.003`$ ebits.

We can now further use the properties of the logarithmic negativity $`E_N`$ to show that such a gap remains even when pure-state entanglement is loaned to assist in the transformations, as it was done with $`\rho _b`$ in . This is achieved by considering a distillation process starting from $`N`$ copies of $`\sigma (p)`$ together with $`L`$ copies of $`|\mathrm{\Phi }`$,

$$\sigma ^N|\mathrm{\Phi }\mathrm{\Phi }|^L|\mathrm{\Phi }^{L+NE_D^c(\sigma )},$$
(24)
Where $`E_D^c(\sigma )`$ denotes the distillable entanglement of $`\sigma `$ in the catalytic setting. For any $`N`$ and $`L`$, we have

$`E_N(\sigma ^N|\mathrm{\Phi }\mathrm{\Phi }|^L)`$ $`=`$ $`E_N(\sigma ^N)+E_N(\mathrm{\Phi }^L)`$ (25)
$`=`$ $`NE_N(\sigma )+L,`$ (26)
where we have used additivity of $`E_N`$ and the fact that $`E_N(\mathrm{\Phi })=1`$. This means that even in the large $`N`$ limit, and once the $`L`$ loaned states $`|\mathrm{\Phi }`$ have been discounted from the distillation outcome, at most $`NE_N(\sigma )`$ ebits of entanglement has been distilled, so that even in the catalytic scenario the bound $`E_D^c(\sigma )<E_N(\sigma )`$ holds. Therefore the irreversibility result of the paper also applies to this case.

We have shown that the irreversibility in the asymptotic manipulation of bipartite mixed states is not a phenomena restricted to bound entangled states, by providing an specific example of distillable state with a finite gap between its entanglement cost $`E_C`$ and its distillable entanglement $`E_D`$. Notice that these results legitimate the use of different measures of entanglement, such as $`E_C`$ and $`E_D`$, to quantify, in the asymptotic limit, the resources of entangled mixed-states. The search for an intrinsic irreversibility in the asymptotic manipulation of bipartite systems has motivated, through several contributions —see for instance —, the development of many techniques for the study of entanglement and has certainly implied an important gain in insight. Paradoxically, a remaining open question is now whether a non-trivial example  of bipartite mixed states exist for which the processes of preparation and distillation can be performed in a fully reversible fashion.

The authors acknowledge motivating discussions with C.H. Bennett, D.P. Divincenzo, G. Giedke, D.W. Leung, J. Preskill, J.A. Smolin, B.M. Terhal and R.F. Werner. G.V. thanks C.H. Bennett for the invitation to visit his group at IBM, May 2001, and all the members of the group for hospitality during the visit. This work was supported by the Austrian Science Foundation under the SFB “control and measurement of coherent quantum systems” (Project 11), the European Community under the TMR network ERB–FMRX–CT96–0087, project EQUIP (contract IST-1999-11053), and contract HPMF-CT-1999-00200, the European Science Foundation, and the Institute for Quantum Information GmbH.


format: text

text:
# Enhacement in the dymanic response of a viscoelastic fluid flowing through a longitudinally vibrating tube

## I Introduction

A quantitative theory of propagation of elastic waves in a fluid-saturated porous solid was formulated in the classic papers by Biot . One of the major findings of Biot’s work was that there was a breakdown in Poisseuille flow above a certain characteristic frequency

specific to the fluid-saturated porous material. Biot theoretically studied this phenomenon by considering the flow of a viscous fluid in a tube with longitudinally oscillating walls under an oscillatory pressure gradient. Apart from the fundamental interest, the investigation of the dynamics of fluid in porous media under oscillatory pressure gradient and oscillating pore walls is of prime importance for the recently emerged technology of acoustic stimulation of oil reservoirs . For example, it is known that the natural pressure in an oil reservoir generally yields no more than approximately 10 percent oil recovery. The residual oil is difficult to produce due to its naturally low mobility, and the enhanced oil recovery operations are used to increase production. It has been experimentally proven that there is a substantial increase in the net fluid flow through porous space if the latter is treated with elastic waves. However, there is a fundamental lack of understanding of the physical mechanisms of fluid mobilization in saturated rock through the effect of elastic waves; the theory of such mobilization virtually does not exist. Biot’s theory can be used to describe the interaction of fluid-saturated

solid with the sound for a classic Newtonian fluid; however, oil and other hydrocarbons exhibit significant non-Newtonian behavior . The aim of this paper is therefore to incorporate non-Newtonian effects into the classical study of Biot .

Recently, del Rio, Lopez de Haro, Whitaker  presented a study of enhancement in the dynamic response of a viscoelastic (Maxwell) fluid flowing in a stationary (non-oscillating) tube under the effect of an oscillatory pressure gradient. We combine this theory with the effect of the acoustic oscillations of the walls of a tube introduced by Biot , providing a complete description of the interaction of Maxwell fluid, filling the pores, with acoustic waves.

Finally, in order to emphasize that the consept of dynamic permeability is an adequate way of desciprtion of the phenomenon, we note that this concept has been widely used before .

In the next section we formulate our model, whereas the last section will conclude with the discussion of the results.

## II The model

In this section we present our model of Maxwell fluid flowing in a cylindrical tube whose walls are oscillating longitudinally and the fluid is subject to an oscillatory pressure gradient. We give analitical solutions of the problem in the frequency domain.

The governing equations of the problem consist of the continuity equation for the incompressible fluid,

$$\stackrel{}{v}=0,$$
$`(1)`$
and the linearized momentum equation,

$$\rho \frac{\stackrel{}{v}}{t}=p\stackrel{~}{\tau }.$$
$`(2)`$
Here, $`\stackrel{}{v}`$, $`p`$, $`\rho `$ denote velocity, pressure and mass density of the fluid, whereas $`\stackrel{~}{\tau }`$ represents the viscous stress tensor. We describe the viscoelatic effects of the fluid using Maxwell’s model, which assumes that

$$t_m\frac{\stackrel{~}{\tau }}{t}=\eta \stackrel{}{v}\stackrel{~}{\tau },$$
$`(3)`$
where $`\eta `$ is the viscosity coefficient and $`t_m`$ is the relaxation time.

Now, let $`u`$ be a velocity of the wall of the tube which oscillates in time as $`e^{i\omega t}`$, where $`\omega `$ is the angular frequency. The flow of fluid in a cylindrical tube with longitudinally oscillating walls can be descibed by a single component of the velocity, namely, its $`z`$-component $`v_z`$ ($`z`$ axis is along the ceterline of the tube). We use the cylindrical coordinate system $`(r,\varphi ,z)`$ in the treatment of the problem. We introduce the relative velocty $`U_1`$ as $`U_1=v_zu`$. Thus, assuming that all physical quantities vary in time as

$`e^{i\omega t}`$, we arrive at the following master equation for $`U_1`$

$$^2U_1+\frac{\omega ^2t_m+i\omega }{\nu }U_1=\frac{X}{\nu }(1i\omega t_m).$$
$`(4)`$
Here, we have introduced the following notations:

$$X=(p+\rho \frac{u}{t}),$$

which is a sum of the applied pressure gradient and force exerted on the fluid from the oscillating wall of the tube and $`\nu `$, which is $`\nu =\eta /\rho `$. Note that a no-slip boundary condition at the wall is assumed.

The solution of Eq.(4) can be found to be

$$U_1=\frac{X}{i\omega }+CJ_0(\beta r),$$

where $`J_0`$ is the Bessel function and $`\beta =\sqrt{(\omega ^2t_m+i\omega )/\nu }`$.

Applying a no-slip boundary condition $`U_1(a)=0`$ at the wall of the tube, where $`a`$ is its radius, we finally get

$$U_1(r)=\frac{X}{i\omega }\left[1\frac{J_0(\beta r)}{J_0(\beta a)}\right]=\frac{Xa^2(1i\omega t_m)}{\nu }\frac{1}{(\beta a)^2}\left[1\frac{J_0(\beta r)}{J_0(\beta a)}\right].$$
$`(5)`$
Defining the cross-section averaged velocity as

$$\overline{U}_1=\frac{2}{a^2}_0^aU_1(r)r𝑑r,$$

we obtain

$$\overline{U}_1=\frac{Xa^2(1i\omega t_m)}{\nu }\frac{1}{(\beta a)^2}\left[1\frac{2J_1(\beta a)}{(\beta a)J_0(\beta a)}\right]K(\omega )X.$$
$`(6)`$
Here $`K(\omega )`$ is the dynamic permeability  that describes the frequency dependent response of the tube to the applied total force on the fluid. Simple comparison reveals that Eq.(6) resembles closely Eq.(6) from Ref. , with the only difference being that we have $`X`$ in place of their $`p/z`$. Note that in the case of stationary tube walls ($`\rho (u/t)0`$) Eq.(6) exactly coincides with Eq.(6) from Ref.. Simple calculation (applying L’Hospital rule for the $`0/0`$ uncertanty) shows that $`lim_{\omega 0}K(\omega )=a^2/(8\nu )`$. Thus, following Ref., we will introduce the dimensionless dynamic permeability as $`K^{}(\omega )=K(\omega )/K(0)`$, which will be used later (see Fig.7). Note, that we were able to easily reproduce Fig. 1 from Ref. , confirming the existence of sharp resonances of $`K^{}(\omega )`$ in the elastic regime (see below for the definition of this regime).

Following the work of Biot  we calculate the stress at the wall $`\tau `$,

$$\tau =\frac{\eta }{1i\omega t_m}\left(\frac{U_1(r)}{r}\right)_{r=a}=\frac{\eta \beta X}{i\omega (1i\omega t_m)}\frac{J_1(\beta a)}{J_0(\beta a)}.$$
$`(7)`$
Note, that when $`t_m0`$ this expression obviously coincides with the corresponding Newtonian form.

The total friction force is $`2\pi a\tau `$. Following Biot we calculate the ratio of total friction force to the average velocity, i.e.

$$\frac{2\pi a\tau }{\overline{U}_1}=\frac{2\pi \eta (\beta a)[J_1(\beta a)/J_0(\beta a)]}{(1i\omega t_m)}\left[1\frac{2J_1(\beta a)}{(\beta a)J_0(\beta a)}\right]^1.$$
$`(8)`$
Simple analysis reveals that

$$\underset{\omega 0}{lim}\frac{2\pi a\tau }{\overline{U}_1}=8\pi \eta ,$$

which corresponds to the limiting case of Poiseuille flow. Following Biot , we also introduce a function $`F(\kappa )`$ with $`\kappa `$ being $`\kappa =a\sqrt{\omega /\nu }`$ in the following manner

$$\frac{2\pi a\tau }{\overline{U}_1}=8\pi \eta F(\kappa ),$$

thus,

$$F(\kappa )=\frac{1}{4}\frac{\kappa \sqrt{i+\kappa ^2/\alpha }\left[J_1(\kappa \sqrt{i+\kappa ^2/\alpha })/J_0(\kappa \sqrt{i+\kappa ^2/\alpha })\right]}{(1i\kappa ^2/\alpha )}\left[1\frac{2J_1(\kappa \sqrt{i+\kappa ^2/\alpha })}{\kappa \sqrt{i+\kappa ^2/\alpha }J_0(\kappa \sqrt{i+\kappa ^2/\alpha })}\right]^1.$$
$`(9)`$
Note, that $`F(\kappa )`$ measures the deviation from Poiseuille flow friction as a function of frequency parameter $`\kappa `$, as introduced by Biot .

In Eq.(9), $`\alpha `$ denotes the Deborah number , which is defined as the ratio of characteristic time of viscous effects $`t_v=a^2/\nu `$ to the relaxation time $`t_m`$, i.e. $`\alpha =t_v/t_m=a^2/(\nu t_m)`$.

As noted in Ref. , the value of the parameter $`\alpha `$ determines in which regime the system resides. Beyond a certain critical value ($`\alpha _c=11.64`$), the system is dissipative, and viscous effects dominate. On the other hand, for small $`\alpha `$ ($`\alpha <\alpha _c`$) the system exhibits viscoelastic behavoir which we call the elastic regime.

Note, that the Newtonian flow regime can be easily recovered from Eq.(9) by putting $`\alpha \mathrm{}`$. We plot this limiting case in Fig. 1. As it can be seen from the figure, both $`Re[F(\kappa )]`$ and $`Im[F(\kappa )]`$ exactly coincide with the Newtonian limiting case studied in Biot’s work (see Fig. 4 in Ref.). This graph demonstrates a breakdown in Poiseuille flow as frequency increases (recall that $`\kappa \sqrt{\omega }`$). In all our calculations we have used polynomial expansions of $`J_0`$ and $`J_1`$ with absolute error not exceeding $`10^6`$ %. Thus, our calculation results are accurate to this order.

A finite-but-large $`\alpha `$ regime is shown in the next two figures. Fig. 2 corresponds to the case when $`\alpha =10^4`$. We see in Fig. 2 that the real and imaginary parts of $`F(\kappa )`$ start to deviate from the Newtonian fluid behavior at large frequencies. In Fig. 3 solutions correspond to the case when $`\alpha =100.0`$; thus we see how viscoelastic effects become pronounced already at low frequencies.

Fig. 4 presents the behavior of $`Re[F(\kappa )]`$ and $`Im[F(\kappa )]`$ when $`\alpha =10.0`$. As it can be seen from the graph, sharp resonances appear on the curves. This feature can be explained by the fact that in this case $`\alpha `$ is less than $`\alpha _c`$, which means that the system switched from the dissipative (viscous) regime to the elastic one.

Fig. 5 shows the solutions when $`\alpha =1.0`$. In this case we see more irregular behavior of $`Re[F(\kappa )]`$ and $`Im[F(\kappa )]`$ with a number of irregular spikes.

The extreme non-Newtonian (elastic) regime is studied in Fig. 6 where we plot the solutions for the case when $`\alpha =10^3`$. In this case notable change is that there are fewer but more pronounced spikes. $`Re[F(\kappa )]`$ is close to zero for most of the frequencies, and only at certain frequencies we see sharp resonances. In this regime the system acts a window for these frequencies.

Another noteworthy observation is that, on one hand, as long as $`\alpha >\alpha _c=11.64`$ (Figs. 2 and 3), $`Re[F(\kappa )]`$ is always greater than its initial value, i.e. $`Re[F(\kappa )]>1`$, and for large frequencies it reaches a ceratin asymptotic value. On the other hand, when $`\alpha <\alpha _c`$ (Figs. 4-6), we observe an overall decrease in $`Re[F(\kappa )]`$ with the increase in frequency, i.e., $`Re[F(\kappa )]<1`$ for all $`\kappa `$’s.

In Fig. 7 we have also studied the behavior of the dimensionless dynamic permeability $`K^{}(\omega _{})`$ as a function of $`\omega _{}`$, the dimensionless frequency defined as $`\omega _{}=t_m\omega `$, for the case when $`\alpha =0.1`$. Since $`\alpha <\alpha _c`$, we observe sharp resonances in the dynamic permability at certain frequences, which can be explained by the non-Newtonian behavior of the fluid.

## III discussion

In this paper we have studied non-Newtonian effects on the dynamics of fluids in porous media by calculating a flow of Maxwell fluid in a tube, which oscillates longitudinally and is subject to an oscillatory pressure gradient. The present study investigates novelties brought about into the classic Biot’s theory  of propagation of elastic waves in a fluid-saturated porous solid by inclusion of non-Newtonian effects. We have used time Fourier transform and converted the problem to the frequency domain. We have calculated the dynamic permeability, thus modifying the work of del Rio, Lopez de Haro, Whitaker  by inclusion of the effect of longitunally oscillating tube walls. We investigated how the function $`F(\kappa )`$, which measures the deviation from Poiseuille flow friction as a function of frequency parameter $`\kappa `$, is modified by non-Newtonian effects. Our work thus provides a combined theory of non-Newtonian flow in a longitudinally oscillating tube, which constitutes the basis for a realistic model of the effects of elastic waves in a fluid-saturated porous space. The present analisys clearly demonstrates the existence of a transition from a dissipative to elastic regime (as $`\alpha `$ decreases), in which sharp enhancements of flow (resonances) occur.

The importance of the current work is two-fold:

(A) We studied modifications brought about by non-Newtonian effects into Biot’s theory. The investigation of the function $`F(\kappa )`$ is important for a number of applications, since $`F(\kappa )`$ uniquely determines the response of a realistic fluid-saturated porous medium to the elastic waves. Thus, determination of $`F(\kappa )`$ for non-Newtonian (Maxwell) fluid is necessary to guide e.g., oil-field exploration applications.

(B) As we have seen from Fig. 7, non-Newtonian effects cause substantial enhancements in dynamic permeability. We were not able to determine in the literature what value of $`\alpha `$ a natural crude oil would have. However, as it can be seen from Fig.7 for $`\alpha =0.1`$, we can get up to 60 times increase in permeability at certain resonant frequences. Lower $`\alpha `$’s yield even more drastic emhancements. At any rate, we obtained an analytical expression for $`F(\kappa )`$ (Eq.(9)), which can provide the behavior of this function for any given $`\alpha `$.

In the introduction section, the practical impact of the possibility of acoustic stimulation of oil reservoirs has been outlined. This result clearly demonstrates that in the crude oil, that can be modelled as a Maxwell fluid, there are certain resonant frequencies at which oil production can be increased significantly if the well is irradiated with elastic waves at these frequencies.

## ACKNOWLEDGMENTS

This work was supported by the Iowa State University Center for Advanced Technology Development and ETREMA Products, Inc.

Figure captions

Fig. 1 Behavior of $`Re[F(\kappa )]`$ (solid line) and $`Im[F(\kappa )]`$ (dashed line), as function of $`\kappa `$ according to Eq.(9). Here, $`\alpha =\mathrm{}`$.

Fig. 2 Same as in Fig. 1 but for $`\alpha =10^4`$.

Fig. 3 Same as in Fig. 1 but for $`\alpha =100.0`$.

Fig. 4 Same as in Fig. 1 but for $`\alpha =10.0`$.

Fig. 5 Same as in Fig. 1 but for $`\alpha =1.0`$.

Fig. 6 Same as in Fig. 1 but for $`\alpha =10^3`$.

Fig. 7 Behavior of $`Re[K^{}(\omega _{})]=Re[K(\omega _{})/K(0)]`$ (solid line) and $`Im[K^{}(\omega _{})]=Im[K(\omega _{})/K(0)]`$ (dashed line), as a function of $`\omega _{}`$ according to Eq.(6). Here $`\alpha =0.1`$.


text:
# Application of Chaos Induced Near-Resonance Dynamics to Locate the Global Optimum of Functions

## Abstract

The problem of locating the global optimum of functions is studied in a dynamic setting. The dynamics of simple multistable systems under the influence of chaotic forcing is investigated. When the magnitude of the forcing signal decays slowly, it is shown that the system attains an equilibrium state, which corresponds to the global optimum of the corresponding multimodal potential function. The role of bifurcations in facilitating the approach to the global optimum is discussed.

The problem of locating the global optimum (specifically the minimum) $`𝐱=𝐱^{}`$, of a continuous, convex function $`f(𝐱)`$ of $`n1`$ variables $`𝐱`$ is considered. Solving the equation(s) representing the condition that the derivatives vanish at the extrema of the function, i.e.

$$\frac{f}{x_i}=0\mathrm{for}i=1,2,\mathrm{},n$$
(1)
is one among several traditional methods available for solving this problem. However. the solutions obtained using almost all of the traditional methods depend on the choice of the initial guess and often correspond to a local optimum.

Alternately, the given objective function can be recast as a system of ordinary differential equations (ODE’s). This is possible if we consider the set of equations (1) to be a set representing the equilibrium of a set of ODE’s in the variables $`𝐱`$, i.e.

$$\frac{dx_i}{dt}=\frac{f}{x_i}$$
(2)
with $`x_i(t=0)=x_i(0),i=1,2,\mathrm{},n`$. Depending on the choice of the initial conditions, the equilibrium point attained by the system will correspond to one of the minima of (1) and one of these will be the global minimum of the objective function (1). In what follows, the equilibrium point corresponding to the global minimum will be referred to as the desired equilibrium point. Even for relatively simple problems, the initial conditions which guarantee an approach to the desired equilibrium point are often difficult to estimate.

Here we study the dynamics of (2) when it is subjected to chaotic forcing. Specifically, we focus on the case where the magnitude of forcing is reduced in a gradual manner. Our objective is to investigate the effect of the magnitude of forcing and the rate of its decay on the transient dynamics and the eventual equilibrium point attained by (2). The system we study is

$$\frac{dx_i}{dt}=\frac{f}{x_i}+C_ie^{k_2t}\mathrm{for}i=1,2,\mathrm{},n$$
(3)
where $`C_i`$ is the chaotic forcing term and $`k_2`$ controls the rate of decay of the forcing. The Lorenz equation with parameters $`\sigma =10;r=60;b=8/3`$ has been used to generate the chaotic signal.

In this article we investigate the evolution of simple multistable systems under the influence of chaotic forcing. The main result is that with appropriate tuning of the chaotic forcing signal, the system attains an equilibrium state, which corresponds to the global optimum of the corresponding multimodal potential function. This formulation can be considered to be a deterministic analog of the simulated annealing optimization technique , with $`C_i`$ and $`k_2`$ being analogous to the temperature and the rate of cooling, respectively.

The optimization problem (1) is thus transformed to a form which has been studied in the context of stochastic resonance (SR) . SR refers to the phenomenon wherein the response of a randomly perturbed nonlinear system to a very weak periodic signal is significantly increased with appropriate tuning of the noise intensity. The mechanism of SR was initially introduced as a possible explanation of the long time climate changes. In the past few years, this phenomenon has been discovered and applied to several situations in physics (in a ring laser and a variety of other bistable systems) and biology (neurological processes). A typical example is in the evolution of the following system

$$\frac{dx}{dt}=\frac{V(x,t)}{x}+\sigma \eta $$
(4)
where $`\eta `$ is white noise, $`\sigma `$ is its magnitude and $`V`$ is a time-periodic double well potential

$$V(x,t)=\frac{x^4}{4}+\frac{x^2}{2}+Ax\mathrm{cos}(\omega t)$$
(5)
Two time scales characterize the evolution of this system. The first corresponds to oscillations in one of the two potential wells (intra-well dynamics), while the second corresponds to oscillations between the wells (inter-well dynamics). The latter characterizes the mean time of barrier crossings and accounts for the global dynamics of noise induced transitions between the potential wells. It has to be noted that the amplitude of the periodic modulation signal is sufficiently small in the sense that the system cannot switch between the potential wells in the absence of noise. For the example system given above, it has been shown that there exist a range of values of $`A`$, $`\sigma `$ and $`\omega `$ where the jumps between the two oscillating wells are strongly synchronized, as a consequence of a resonance between the periodic forcing and random perturbation.

In analogy with SR, chaotic resonance (CR) refers to the transitions induced by chaotic forcing of nonlinear systems . The presence of deterministic chaos plays the role of random perturbation. Depending on the magnitude of the chaotic forcing signal, a multistable system can show intra-well dynamics or inter-well dynamics. This is similar to that discussed earlier in the context of SR. Thus under appropriate conditions, the dynamics of the system can be strongly synchronized with the time variation of the periodically varying control parameter. In several systems, CR has been shown to occur near the onset of chaos and near crisis bifurcation points. Briefly, a crisis bifurcation point refers to the bifurcations following which a system undergoes a transition from intra-well to inter-well dynamics .

Inducing a crisis bifurcation is therefore the key step in our objective of locating the global optimum of functions. The proposed scheme involves inducing crisis bifurcations in the system by driving the system with an appropriate amount of a chaotic input. We have verified that signals from a variety of chaotic systems can be used to induce crisis bifurcations in several multistable systems of interest to us. Transitions from intra-well to inter-well dynamics can be easily obtained by increasing the magnitude of the forcing. Once the system starts to show inter-well dynamic behavior, the magnitude of the forcing signal is reduced in a gradual manner. An important result of this work is that for a vast choice of initial conditions. the system (2) stabilizes at the state which corresponds to the global minimum of the potential function $`f(𝐱)`$. The first of the four example minimization problems presented below has been selected with the purpose of illustrating the nature of the bifurcations which make it possible for the system to reach the global minimum. The remaining three examples are well known test functions and have been selected to show the validity of the technique.

The first example is:

$$\mathrm{min}f(x)=\frac{x^4}{4}16x^2+5x$$
(6)
The ODE system which is solved for determining the global minimum is:

$$\dot{x}=f^{}(x)+k_1ce^{k_2t};f^{}(x)=(x^332x+5)$$
(7)
Here $`c`$, $`k_1`$ and $`k_2`$ represent the chaotic forcing signal, its magnitude and the decay constanT, respectively. The negative sign is needed to ensure that the two minima of (6) are the stable equilibrium points of (7).

When $`k_1=0`$, the unforced system (7) has an unstable equilibrium point near $`x=0.40`$ and this corresponds to the maximum value of (6). This equilibrium acts as a potential barrier and prevents an approach to the negative equilibrium point for almost all choices of initial conditions that are positive. The crossover from the local to the global minimum can occur only if the forcing input is sufficient to overcome the potential barrier associated with the presence of this maximum. For small magnitude of forcing and depending on the initial condition, the system (7) oscillates in the vicinity of one of the two stable equilibrium points. As the magnitude is increased past a critical value, the oscillatory state associated with the positive equilibrium point disappears. This is due to the occurrence of a crisis bifurcation This value of the magnitude of forcing is the least value which makes it possible to reach the desired equilibrium state. Following this bifurcation and irrespective of the initial conditions, the system (7) intermittently evolves across its former basin of attraction, which no longer exists. Now the only stable state is that associated with the negative equilibrium point. As the magnitude of forcing is gradually reduced, the system spends more time in the vicinity of this equilibrium point and hence has an increased probability of eventually stabilizing there.

Fig 1(a) and 1(b) show the dynamics of the system before it reaches the local equilibrium (Fig. 1(a)) and global equilibrium (Fig. 1(b)) points. As can be seen, the system oscillates near each minima before stabilizing at the global minimum. For sufficiently low decay rates, the systems always reaches the negative equilibrium solution, which is the global minimum of (6). Finally, it must be noted that the magnitudes of forcing and the decay constant are the crucial parameters which together decide whether or not the system actually reaches the desired state.

To summarize, the ability of this system (a particle in a double-well potential) to jump over the energy barrier is acquired following the occurrence of a crisis bifurcation. Increased sensitivity of the system following the occurrence of a crisis bifurcation (near-resonance dynamics) allows stabilization the system at a state corresponding to the global minimum of the function.

The second example is a slightly more complex function:

$$\mathrm{min}f(x)=(x_1^2+x_211)^2+(x_1+x_2^27)^2+(x3)^2+(x2)^2$$
(8)
This function has four minima with the global minima at $`(x_1=3,x_2=2)`$. We solve the following equation

$`{\displaystyle \frac{dx_1}{dt}}`$ $`=`$ $`(4x_1^3+2x_2^2+4x_1x_240x_120)+k_1ce^{k_2t}`$ (9)
$`{\displaystyle \frac{dx_2}{dt}}`$ $`=`$ $`(4x_2^3+2x_1^2+4x_1x_224x_226)+k_1ce^{k_2t}`$ (10)
The terms within the braces in (9) and (10) are the derivatives with respect to the dependent variables $`x_1`$ and $`x_2`$ respectively. For low decay rates $`(k_2<0.80)`$ and when $`(k_1>1.50)`$, the solutions of these equations evolve (randomly) between the different minima before finally stabilizing at the state corresponding to the global minimum (Fig. 2). This again shows the ease with which the global minimum of multistable systems can be located.

Fig. 3 shows results of a parametric study carried out in $`k_1`$-$`k_2`$ space. 1000 sets of randomly chosen initial conditions were used and the number of initial conditions leading to the global minimum was computed. This figure shows that for sufficiently high values of $`k_1`$ and sufficiently low values of $`k_2`$, the system always stabilizes at the desired equilibrium point.

The third example considered is that of the multimodal Rastrigin function with one and three variables

$$\mathrm{min}f(x)=10n+\underset{i=1}{\overset{n}{}}(x_i^210\mathrm{cos}(2\pi x_i)),n=1\mathrm{or}3$$
(11)
with a global minimum at $`x_i=0,i=1,\mathrm{},n`$. Fig. 4 shows the evolution of the forced system and the approach to the desired equilibrium point. Values of $`k_1>20`$ and $`k_2<0.005`$ were found to be sufficient to stabilize the system at the desired equilibrium point.

A common feature of the functions considered previously is that the global minima lies in the deepest potential well. In such systems, addition of a chaotic forcing signal systematically destroys the local minimum in the order of increasing depth of potential wells. This makes it possible to locate the global minimum of the functions. A similar argument shows that it will not be possible to reach the global minimum using the scheme (3) for systems in which the global minimum and the deepest potential well do not coincide. An example of such a situation is the function (Fig. 5)

$$\mathrm{min}f(x)=10x(x0.2)(x1)(x2)(x3)(x4)+3x^3+x^2$$
(12)
In such cases, on/off tuning of the chaotic forcing signal using a suitably chosen acceptance or rejection criterion based on the magnitude of the function (12) provides a way of locating the global minimum. Fig. 6 shows a typical example of the temporal evolution of the dynamical system obtained from (12).

In conclusion, it has been shown that it is possible to locate the global optimum of convex functions by solving the problem in a dynamic setting. Chaotic forcing of appropriate magnitudes lead to the occurrence of crisis bifurcation(s) in the multimodal system. This enables the system to sequentially overcome the (local) energy barriers. A suitably low rate of decrease of the forcing signal stabilizes the system at the state in which its energy is the lowest.

LIST OF FIGURES

FIG. 1. Evolution of the system (7) with time showing the approach to the (a) local minimum for $`k_1=2.750`$, $`k_2=0.01`$; (b) global minimum for $`k_1=5`$, $`k_2=0.01`$. The chaotic drive is the $`y`$ variable of the Lorenz equation. The dashed lines show the location of each of the two stable equilibrium points.

FIG. 2. Evolution of $`x_1`$ and $`x_2`$ of Equation (9) and (10) with time showing the approach to the global minimum, $`k_1=5`$ and $`k_2=0.05`$. The chaotic drive is the $`y`$ variable of the Lorenz equation. The dashed lines show the location of the desired equilibrium point.

FIG. 3. Parametric dependence of the number fraction of initial conditions $`F`$ stabilizing at the local minima of the Himmelblau function (8).

FIG. 4. Temporal evolution of the dynamic system obtained from (11) showing the approach to the global minimum for $`k_1=25`$, $`k_2=0.01`$. The chaotic drive is the $`y`$ variable of the Lorenz equation. The dashed line shows the location of the desired equilibrium point.

FIG. 5. Graphical representation of the function (12). This is an example of the case where the global minimum does not coincide with the deepest potential well.

FIG. 6. Temporal evolution of the dynamic system obtained from (12) showing the approach to the global minimum for $`k_1=10`$, $`k_2=0.008`$. The simulation was carried out using an on/off control strategy. The chaotic drive is the $`y`$ variable of the Lorenz equation. The dashed line shows the location of the three stable equilibrium points.


text:
# Velocity-induced collapses of stable neutron stars

## 1 Introduction

In the early eighties Shapiro and Teukolsky (1980) addressed the question of formation of black holes in astrophysical collapses. Among other, they asked the (theoretical) question of the minimum mass of a black hole formed by the adiabatic collapse of a stellar core and “In particular, can the effective mass-energy potential barrier associated with equilibrium configurations be penetrated by low-mass cores with substantial inward, radial kinetic energy?”. In astrophysical scenarios, black holes can also form from accretion-induced collapses of neutron star, if the neutron star is part of a binary system or during a supernova event, when part of the envelop fails to reach escape velocity and falls back onto the new born neutron star. Therefore, one can also ask the question: how much of inward kinetic energy has to be put to a neutron star to make it collapse to a black hole? Can a neutron star always collapse to a black hole, provided that it gets enough kinetic energy?

If one looks at static neutron star models, these are stable against perturbations if their masses are lower than some maximal mass and the density is also lower than the “critical” density corresponding to this maximal mass (they are located on the so-called stable branch). They can collapse to form a black hole when the central density is higher than the critical one (see e.g. Gourgoulhon (1991)). The question is then whether stable neutron stars (i.e. with lower density and mass than the critical ones) can be “pushed”, with a certain amount of inward kinetic energy, to form a black hole, and what is the minimal mass of the formed black hole and/or progenitor? The problem of the minimal mass of a black hole has been solved, first from a more mathematical point of view, by Choptuik (1993), who discovered the “critical collapse” phenomena in General Relativity (see also Sec. 4.1 for a discussion). Still, the problem of a stellar core/neutron star with inward kinetic energy have never been completely studied; to our knowledge, only a partial study was done by Gourgoulhon (1992). The aim of this paper is to numerically follow the collapse of stable neutron stars, with inward velocity profile, and determine the initial conditions necessary to obtain a black hole as result of the collapse; the purpose is also to allow for a realistic equation of state (EOS), to study dependence on the EOS, as well as on the initial velocity profile, giving kinetic energy. Therefore, emphasis shall be put on neutron star properties in General Relativity and the question of the physical mechanism giving this kinetic energy will not be addressed.

The complete model in General Relativity is described in Sec. 2, including the system of partial derivative equations (2.1), the two EOS used (2.2), and the procedure for obtaining initial numerical models (2.3). Time evolution is studied in Sec. 3, where two different possible evolutions for neutron stars are given (3.1 and 3.2). Numerical results are displayed in Sec. 4: together with a link with the “critical collapse” paradigm (4.1), properties of initial data (4.2) and dependence on the parameters (4.3) are discussed. Finally, Sec. 5 summarizes the results and gives some concluding remarks.

## 2 Evolution of spherically symmetric neutron stars

### 2.1 Field and Matter equations

The equations for the evolution of the matter and gravitational fields are derived from the Einstein equations using the simplifying assumption of spherical symmetry. Space-time is foliated into spacelike hypersurfaces $`\mathrm{\Sigma }_t`$ and equation are written in terms of the classical 3+1 formalism of General Relativity (see e.g. Arnowitt et al. (1962)). The real parameter $`t`$ is called the coordinate time and it can be shown that, making the choice of polar time slicing and radial gauge, in spherical symmetry the metric can be expressed in the form (generalization of the Schwarzschild metric, see also Gourgoulhon (1991)):

$`g_{\mu \nu }\mathrm{d}x^\mu \mathrm{d}x^\nu `$ $`=`$
$`N^2(r,t)\mathrm{d}t^2+A^2(r,t)\mathrm{d}r^2`$ $`+`$ $`r^2(\mathrm{d}\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\varphi ^2),`$ (1)
where $`N(r,t)`$ is called the lapse function. The two functions $`N`$ and $`A`$ will often be replaced by $`\nu (r,t)`$ and $`m(r,t)`$, defined as111We use geometrized units in which the speed of light $`c`$ and Newton’s gravitational constant $`G`$ are equal to unity:

$$\nu =\mathrm{ln}(N)\mathrm{and}A=\left(1\frac{2m}{r}\right)^{1/2}.$$
(2)
The 4-velocity of the fluid is denoted $`v^\mu `$, the fluid radial coordinate velocity being thus:

$$\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{v^r}{v^0};$$

and the “physical” fluid radial velocity $`U`$, as measured locally by the hypersurfaces observer is defined by:

$$U=\frac{A}{N}\frac{\mathrm{d}r}{\mathrm{d}t}.$$
(3)
Following Romero et al. (1996), the stress-energy tensor was taken to be that of a perfect fluid; so hydrodynamical equations in General Relativity ($`_\mu T_\nu ^\mu `$ = 0) have been written in the form of a system of conservation laws

$$\frac{u}{t}+\frac{1}{r^2}\frac{}{r}\left[r^2\frac{N}{A}f(u)\right]=s(u),$$
(4)
where $`u=\{D,\mu ,\tau \}`$ is the vector of evolved quantities, $`f(u)`$ the vector of fluxes and $`s(u)`$ the vector of sources (for details see Romero et al. (1996)). The evolved quantities are defined from hydrodynamical variables, baryon and total energy densities in the fluid frame ($`n_B`$ and $`e`$) by:

$`D`$ $`=`$ $`A\mathrm{\Gamma }n_B,`$
$`\mu `$ $`=`$ $`(E+p)U,`$ (5)
$`\tau `$ $`=`$ $`ED,`$
where $`p`$ is the fluid pressure given by the EOS (see 2.2), $`\mathrm{\Gamma }=(1U^2)^{1/2}`$ is the Lorentz factor of the fluid, and $`E=\mathrm{\Gamma }^2(e+p)p`$.

The gravitational field equations in radial gauge, polar slicing and spherical symmetry reduce to two equations (with no time evolution of the gravitational field):

$`{\displaystyle \frac{m}{r}}`$ $`=`$ $`4\pi r^2E,`$ (6)
$`{\displaystyle \frac{\nu }{r}}`$ $`=`$ $`A^2\left({\displaystyle \frac{m}{r^2}}+4\pi r(p+(E+p)U^2)\right).`$ (7)
The system therefore consists of three evolution equations for the hydrodynamical variables (Eqs. 4) and two constraint equations for the gravitational field (Eqs. (6)–(7)). Finally, an equation of state $`p(n_B,e)`$ closes the system.

### 2.2 Equation of state

In this work two different EOS for nuclear matter have been used, to describe microscopic properties of neutron stars. Apart from the initial condition models, which were calculated with an equation of state $`p(n_B)`$ for cold, catalyzed matter, the used equations of state were of the form $`p(n_B,e)`$. The first EOS used is the well-known ideal gas model:

$$p=(\gamma 1)(em_Bn_B)$$
(8)
where $`m_B=1.66\times 10^{27}`$ kg is the baryon mass and $`\gamma `$ is the index. We chose $`\gamma =2`$ which may mimic relatively well the properties of neutron star matter but always keeps a sound speed lower than $`c`$ (causal EOS). This EOS has been used not only because it is very convenient from a numerical point of view, but also because it allows for the calibration and test of the code. For obtaining initial conditions, the following expression has been used

$$p(n_B)=Kn_0m_B\left(\frac{n_B}{n_0}\right)^\gamma ;$$
(9)
with $`K=0.1`$ (as in the study of rotating neutron star models by Salgado et al. (1994) ) and

$$n_0=0.1\mathrm{fm}^3.$$
(10)
To get a more realistic EOS for nuclear matter interactions, we used the EOS described in Pons et al. (2000). It is a relativistic field theory model supplemented by nonlinear scalar self interactions. Nucleons $`(n,p)`$ interact via the exchange of $`\sigma `$, $`\omega `$ and $`\rho `$mesons. The contribution from leptons is given by its non-interacting form, since their interactions give negligible contributions. Moreover, the star being at chemical equilibrium with respect to the weak processes, the neutrino chemical potentials are zero. However, as a difference from Pons et al. (2000), kaon interactions are not considered in this work. The EOS was tabulated and, during the time integration, the interpolation has been done using bi-cubic splines. This may not be as precise as the method presented by Swesty (1996) (using bi-quintic interpolation), but thermodynamical consistency is still preserved. At low densities, this EOS has been smoothly joined with a polytrope. As far as initial conditions are concerned, the same model has been used, but the temperature has been set to zero.

### 2.3 Initial Conditions

To get initial numerical conditions, we first obtained stable spherical neutron star models and then, to add a velocity profile. A stable neutron star is defined by the fact that its central density is lower than the critical one ($`n_B^{\mathrm{crit}}`$), defined by ($`M_g`$ being the total gravitational mass of the star and $`n_B^c=n_B(r=0)`$):

$$\frac{\mathrm{d}M_g}{\mathrm{d}n_B^c}|_{n_B^c=n_B^{\mathrm{crit}}}=0.$$
(11)
The stable model is easily computed, integrating the well known Tolman-Oppenheimer-Volkoff (TOV) system, which is the static (all $`/t`$ terms and $`U`$ are set to zero) limit of momentum evolution equation in system (4), plus the equations for gravitational fields Eqs. (6)–(7). For the polytropic equation of state ($`\gamma =2,K=0.1`$), one has $`n_B^{\mathrm{crit}}=3.18`$ $`n_0`$ and the maximal stable mass is $`M_g=3.16\mathrm{M}_{}`$. For the tabulated EOS of Pons et al. (2000) (see previous section), one has $`n_B^{\mathrm{crit}}=10.5`$ $`n_0`$ and the maximal stable mass is $`M_g=2.08\mathrm{M}_{}`$. The rather high value for maximal mass for the $`\gamma =2`$ polytropic EOS is linked to an also high value chosen for $`K`$. Thus mass scales between both EOS used in this work were different. However for polytropes of the form (9) with $`\gamma =1+1/n`$, the constant $`K^{n/2}`$ has units of length in geometrized units. Following Cook et al.  (1994), one can use this constant to set the fundamental length scale of the system. Defining dimensionless quantities as in Cook et al.  (1994) ($`\overline{t}=K^{n/2}t`$, $`\overline{p}=K^np`$, $`\overline{n}_B=K^nn_B`$, …), one can see that in our case all masses (even in dynamical evolution) scale like $`K^{0.5}`$.

Once the density and gravitational fields are computed, the inward velocity profile is added in the r.h.s. of gravitational field equations (6)–(7). The metric coefficients are then relaxed in order to take into account the contribution from the kinetic energy in the total gravitational mass. Therefore, initial conditions are consistent with the gravitational field equations. The considered profiles were of the form

$$U(r)=\frac{A(r)}{N(r)}V\left(\frac{r}{R_{\mathrm{surface}}}\right),$$

where $`V(x)`$ has one of these forms ($`x=r/R_{\mathrm{surface}}`$):

$`V(x)`$ $`=`$ $`{\displaystyle \frac{V_{\mathrm{amp}}}{2}}(x^33x)`$ (12)
$`V(x)`$ $`=`$ $`{\displaystyle \frac{27V_{\mathrm{amp}}}{10\sqrt{5}}}\left(x^3{\displaystyle \frac{5x}{3}}\right)`$ (13)
($`V_{\mathrm{amp}}`$ is the parameter defining the amplitude of the profile).

Both are such that the minimal value in  is $`V_{\mathrm{amp}}`$, but (12) verifies $`V^{}(1)=0`$ and (13) has a null divergence at the surface. These profiles correspond to what is usually observed in collapses of neutron stars (see e.g. Gourgoulhon (1991)) or in the formation of neutron stars (see e.g. Romero et al. (1996)). The initial minimal value of the velocity $`U(r)`$ in the interior of the star ($`r[0,R_{\mathrm{surface}}]`$) is noted $`U_{\mathrm{min}}`$. It shall also be called the initial velocity profile amplitude.

## 3 Dynamical scenarios

Numerical time integration has been done using High Resolution Shock-Capturing schemes (HRSC, see Banyuls et al. (1997)) for the hydrodynamical system (4). The metric constraint equations (6)–(7) were integrated using standard finite-differences methods. The need of using numerical methods able to handle shocks comes from the fact that strong discontinuities can form (see Sec. (3.2)), as it has first been observed by Gourgoulhon (1992). Unfortunately, this last study has been limited by the use of spectral methods, unable to handle shocks.

### 3.1 Direct collapse

For a stable star close to the maximal mass, but with an inward velocity profile, there may be two final issues. If the velocity is relatively small in amplitude, the star enters an (theoretically) infinite series of oscillations, for no viscous nor radiative damping is present. With higher velocity, a collapse to a black hole occurs “normally”, almost like for an unstable neutron star configuration (as in e.g. Gourgoulhon (1991)).

Figure 1 shows energy ($`E`$) and velocity profiles at different moments of the collapse of such a neutron star. It corresponds to an initial configuration of 2.74 $`\mathrm{M}_{}`$, where the initial velocity profile is of the type (13), with $`V\mathrm{amp}=0.4c`$, which gives an initial minimal radial velocity of $`U_{\mathrm{min}}=0.6c`$. This collapse has been followed up to $`t=0.48`$ ms, when then central value of the lapse became $`5.13\times 10^5`$, with a radial velocity at the outer edge of -0.99 $`c`$. Due to the choice of polar slicing which avoids the appearance of singularities, the horizon of a black hole cannot be numerically described. However, the final state is a black hole, but from the “frozen star” viewpoint, with a collapse of the lapse. These criteria have been used in Sec. 4 to determine whether the result of the collapse was a black hole or not, namely a central value of the lapse $`N(r=0)<10^4`$ with an ingoing radial velocity. If both of these conditions are fulfilled, it seems unlikely that the collapse lead to anything else than a black hole.

Still, a difference from Gourgoulhon (1991) in the dynamical evolution can be seen in Fig. 1: a strong gradient appears in the velocity profile around $`r=9`$ km, even quite before the “frozen star” regime (where the metric potentials $`A`$ and $`N`$ also exhibit strong gradients). This had already been observed by Gourgoulhon (1992), where it limited the study. It is shown in Sec. 3.2 that, with higher initial velocities, a shock can form and part of the infalling matter can be ejected. In the present “direct collapse” regime all the matter ends in the black hole so that the mass of the resulting object is that of the initial neutron star. In Sec. 4, this has been used to distinguish between both possible types of collapses.

### 3.2 Shock and bounce

For some initial velocities, the “strong gradient” of previous section turns into a shock, and part of the infalling material is ejected.

This can be seen on Fig. 2 where energy density and radial velocity profiles are displayed, for an initial neutron star of 1.91 $`\mathrm{M}_{}`$ (the mass is lower than that of Fig. 1 due to a lower initial central baryon density). The shock appears around $`t=t_3`$, $`r=10`$ km and then moves out of the numerical grid ($`t_4`$) due to the accumulation of matter still falling at $`r>10`$ km. The velocity on the left side of the discontinuity grows and reaches positive values. Later during the collapse, one can see that part of the falling matter starts moving outward and reaches escape velocity from the central object $`t_5`$. This of course reduces the possibility of forming a black hole, since matter is spatially less concentrated. Still, the dynamics of the central region (where the velocity is still inward) can proceed to a black hole, which shall then accrete matter which has not reached escape velocity. This has been observed in the collapse displayed on Figs. 2 where, at the end of the computation ($`t=0.63`$ ms), the central lapse was $`N(r=0)=1.18\times 10^9`$ and the velocity on the edge of the “frozen star” equal to -0.999 $`c`$. The central region would collapse “directly”, as described before in Sec. 3.1. The final mass of the black hole was 1.61 $`\mathrm{M}_{}`$ (and the areal radius in RGPS coordinates equal to 4.77 km), which shows that almost 16% of the initial matter of the neutron star has been ejected.

## 4 Numerical results

The dynamical fate of neutron stars with initial velocity are now studied, keeping in mind the different scenarios described in previous section and varying the mass of the progenitor, as well as the amplitude of the initial velocity profile. This exploration of the parameter space shows the possibility of formation of very low mass black holes, which can be seen as a feature of the “critical collapse” phenomena.

### 4.1 Link with critical collapses

The critical collapse phenomenon (for an interesting review see Gundlach (1999)) was discovered in the early 90s by Choptuik (1993), who was numerically studying the gravitational collapse of spherically symmetric massless scalar field. Depending on some parameter of the initial conditions (which is generically noted $`p`$), the final result of the collapse would be a black hole ($`p`$ large) or the dispersion of the field ($`p`$ small). He discovered that, when fine-tuning this parameter, he could get black holes of arbitrarily small masses. Moreover, the relation giving the mass of the resulting black hole, close to $`p_{}`$ the minimal value to form a black hole, appeared to be universal in the form:

$$M_{BH}C(pp_{})^\alpha ,$$
(14)
where $`C`$ is a constant and $`\alpha `$ is called the critical exponent. The space-time obtained with $`p=p_{}`$ shows the very interesting geometrical property of self-similarity.

There have been many works on critical collapses since this pioneering one and, in particular, the study concerning perfect fluids by Neilsen and Choptuik (2000) is of particular interest for our study. It is shown that perfect fluid collapses, in ultra-relativistic regime, also exhibit “critical” behavior. Our work can easily be connected to this one, if one considers, at fixed central density of the initial neutron star, the parameter $`p`$ to be e.g. the parameter $`V_{\mathrm{amp}}`$ or the amplitude of the initial velocity profile $`U_{\mathrm{min}}`$. One can then see that, by fine-tuning the parameter $`V_{\mathrm{amp}}`$ ($`U_{\mathrm{min}}`$ being a monotonic and continuous function of $`V_{\mathrm{amp}}`$) of (12) or (13), one could get black holes of arbitrarily small masses, as a result of the velocity-induced collapse. We have not investigated the domain of very small black hole masses, for our code was not designed for it. Still, starting with a $`1.16\mathrm{M}_{}`$ neutron star and an initial velocity profile given by (13) (and $`V_{\mathrm{amp}}=0.79732`$), the result of the evolution was a $`3.7\times 10^2\mathrm{M}_{}`$ black hole. This corresponds to the lowest black hole mass range (about $`10^2\mathrm{M}_{}`$) obtained with our code. In order to test the code against the relation (14), the central density of the initial configuration has been kept constant, and only the parameter $`V_{\mathrm{amp}}`$ of (13) has been varied (and thus $`U_{\mathrm{min}}`$).

The results displayed on Fig. 3 correspond to the masses of black holes resulting from collapses, if the central density of the initial configurations is held fixed to $`0.1n_0`$. These masses are drawn (crosses) as a function of the minimal value of the velocity $`U`$ inside the neutron star. The relation (14) is verified, at least for “small” masses and the best fit (shown on the figure in solid line) corresponds to $`\alpha =0.52`$, $`C=10.4\mathrm{M}_{}`$ and $`p_{}=0.8705c`$ ($`U_{\mathrm{min}}=p`$). These values has shown to depend on the particular central density used for the initial configurations. The fit is poor for large masses and, in particular, the relation (14) breaks down completely near the maximal mass for stable static neutron stars (when central density reaches $`n_B^{\mathrm{crit}}`$). The mass scaling exponent $`\alpha =0.52`$ differs from the result obtained by Neilsen and Choptuik (2000) ($`\alpha =1`$) who used an ultrarelativistic EOS with $`\gamma =2`$. Although the result from Neilsen and Choptuik (2000) is rather universal, the study here starts from very different initial conditions (neutron star/exponential distribution of matter in their case). In particular, here, the maximal mass for stable neutron stars gives a mass limit to the problem. Another difference is that the parameter used to get the relation (14) is the amplitude of a velocity profile, which has not been studied by Neilsen and Choptuik (2000). Finally, it has to be stressed out that the family of initial conditions used for this study depends on two parameters (velocity amplitude and central density), the work of Neilsen and Choptuik (2000) assumed only a one-parameter family of initial data. The role played by the central density as a parameter is not that of a “critical” one: if one considers the family of initial data given by static neutron stars (without any velocity profile) of increasing central densities, then the value $`n_B^{\mathrm{crit}}`$ defined by (11) would seem to be the same as $`p_{}`$ in controlling the formation of a black hole, but there is no critical behavior at this point. Still it is interesting to note that the mass scaling relation is valid for velocity induced neutron star collapses, at least in the limit of low masses (and therefore relatively low central densities).

### 4.2 Masses of progenitors - Initial velocities

Once the EOS and type of velocity profile are chosen, a neutron star model serving as initial condition for velocity-induced collapse is completely specified by two parameters: the central baryonic density $`n_B^c`$ and the parameter $`V_{\mathrm{amp}}`$ or, equivalently, the neutron star gravitational mass and the amplitude of the initial velocity profile $`U_{\mathrm{min}}`$. For each point ($`n_B^c`$, $`U_{\mathrm{min}}`$) in this parameter space the question of knowing whether the corresponding neutron star would collapse to a black hole have been addressed. First, it has been found that there exists a maximal value for $`U_{\mathrm{min}}`$, for which an initial model of neutron star could be computed. For higher velocities, it appeared to be numerically impossible to compute the metric coefficients taking into account this large amount of kinetic energy: the relaxation would fail to converge, leading to velocity $`U>c`$. This may come from the way these initial conditions are set: the addition of a velocity profile to a static model does not give a well defined result. Matter is not at equilibrium, whereas the metric potentials are static. Therefore an improvement of this study would be to consider more realistic initial conditions (as for example in Hawley and Choptuik (2000)), where the velocity profile is not set ad hoc but comes from a dynamical interaction.

This maximal value of $`U_{\mathrm{min}}`$ depends on the central density of the neutron star, as displayed on Fig. 4; one can see that for central densities going to zero (and therefore also masses), this maximal value goes to one. This figure also shows the various fates of a neutron star in the ($`n_B^c`$,$`U_{\mathrm{min}}`$) plane, for the ($`\gamma =2`$) polytropic EOS of Sec. 2.2 and an initial velocity profile of type (13).

Curves are described from top-right to bottom-left part of the figure. A neutron star set to a point of this parameter space can become a black hole or not. Still, some region is not really interesting to study. The second curve to be displayed selects neutron stars with a mass lower than the maximal one, as defined in Sec. 2.3. Studying collapses for neutron stars with higher masses is not relevant since one knows that in that case the neutron star shall end as a black hole; the additional velocity profile of the neutron star only gives kinetic energy so that the total mass of the star become larger than the maximal one for stars at rest. It was then observed that this collapse generally proceeds “directly”, i.e. no shock is present. One more curve gives the limit between the direct collapse domain and the part of the parameter space for which the star undergoes a shock and bounce, as described in Sec. 3.2. Finally, if the parameters of the initial neutron star model are below the last curve, the black hole does not form, the matter being ejected to infinity.

The same regions are displayed in the $`(M_g,U_{\mathrm{min}})`$ plane, on Fig. 5. In particular, one can see that there exists a minimal mass for the initial neutron star, to form a black hole. This mass is, in the case of a polytrope $`\gamma =2`$ and an initial velocity profile (13):

$$M_{\mathrm{min}}=1.164\mathrm{M}_{}.$$
(15)
This figure also shows the minimal neutron star mass for which direct collapse (i.e. no ejection of matter) can occur, it found to be 2.276 $`\mathrm{M}_{}`$.

### 4.3 Dependence on EOS and velocity profiles

Although the velocity profile (13) is supposed to occur in most of physical scenarios of gravitational collapse of compact stars, one has to check the dependence of results of previous sections on a change of this velocity profile. Therefore the study of Sec. 4.2 has been undertaken with the initial velocity profile for neutron stars being given by (12). For initial configurations, one can see that this change in the velocity profile induces, for given values of $`n_B^c`$ and $`U_{\mathrm{min}}`$, a change in the distribution of kinetic energy in the star and thus a change in the total gravitational mass $`M_g`$, with respect to the velocity distribution given by (13). Still, qualitative results concerning collapses were retrieved: direct collapse or bounce as dynamical scenarios and formation of low-mass black hole by fine-tuning the parameter $`U_{\mathrm{min}}`$. The scaling relation (14) has been studied for a fixed central density (chosen to be $`0.1n_0`$ as in Sec. 4.1) and the best fit of results gives the relation:

$$M_{BH}6.98\times (U_{\mathrm{min}}0.8991)^{0.52}\mathrm{M}_{},$$

for which the exponent ($`\alpha =0.52`$) is the same as that obtained using velocity profile (13). The two other coefficients (noted $`C`$ and $`p_{}`$ in Sec. 4.1) are different for the reasons stated above, concerning the difference in kinetic energy distribution. The borders of different regions of Fig. 5 also slightly change, nevertheless we get as minimal mass for the initial neutron star to form a black hole for a $`\gamma =2`$ polytrope and an initial velocity profile (12):

$$M_{\mathrm{min}}=1.155\mathrm{M}_{}.$$
(16)
This is less than 1% different from the minimal mass for velocity profile (13), that is within global numerical error. One then may suppose that this quantity is little dependent on the particular type of velocity profile chosen.

The second global “parameter” one would like to change is the EOS for nuclear matter, since a rather important uncertainty exists on properties of neutron star matter. In particular, the $`\gamma =2`$ polytrope used in this study had the advantage of rapid computation as well as numerical stability, combined with a rather good description of neutron matter properties. But if one wants to go further, it is then necessary to use a more realistic EOS, as for example the one described in Sec. 2.2. As stated in that section, maximal mass for static neutron stars is found to be $`M_g^{\mathrm{crit}}=2.08\mathrm{M}_{}`$, for a corresponding central density of $`n_B^c=10.5n_0`$. Therefore, results from computations using this second EOS would show lower masses and higher central densities than those using the analytical EOS. The parameter space ($`n_B^c`$, $`U_{\mathrm{min}}`$) has been studied, using the velocity profile given by (13). As for the $`\gamma =2`$ polytrope, neutron stars with realistic EOS would collapse to form a black hole, provided that the amplitude of the velocity profile be large enough. This collapse could occur directly, as in the case of unstable neutron stars or with the bounce of a part of the matter, which would be ejected, allowing only for the central region to form a black hole. The mass scaling relation is recovered, for relatively small initial central densities; for example, with $`n_B^c=1.5n_0`$, the following relation is found:

$$M_{BH}5.79\times (U_{\mathrm{min}}0.7519)^{0.71}\mathrm{M}_{},$$

which has been tested for black holes in the mass range $`5\times 10^30.7\mathrm{M}_{}`$. The global topology of regions described in figures 4 and 5 is kept and the minimal mass for a neutron star (with matter described by this realistic EOS) to form a black hole is found to be

$$M_{\mathrm{min}}=0.36\mathrm{M}_{}.$$
(17)
During the collapses leading to very low-mass black holes (a few $`10^2\mathrm{M}_{}`$), the central regions of the star would reach very high densities and it has been checked that sound velocity had never become higher than $`c`$.

## 5 Summary and conclusions

With an important amount of inward kinetic energy, stable neutron stars may collapse to a black hole, overcoming the potential barrier which separates both types of objects. However, there seem to exist an absolute lower mass limit, depending on the particular equation of state used to describe neutron star matter properties, below which neutron stars cannot pass this barrier. This lower-limit mass has been found to be $`1.16\mathrm{M}_{}`$ for a $`\gamma =2`$ polytropic EOS, and $`0.36\mathrm{M}_{}`$ for a more realistic model of nuclear matter. It seems also that these values are independent from the particular initial velocity profile added to static neutron star models. But, as it has been noted in Sec. 2.3, for the polytropic EOS, like all masses this lower-limit mass scales like $`K^{0.5}`$.

For both studied EOS, in the case when the star collapses to a black hole it may either proceed directly, as if it were an unstable neutron star, all the matter ending in the black hole; or there may appear a shock and bounce, ejecting a part of the matter to infinity, so that only a fraction of the initial neutron star forms a black hole. With such a mechanism, the resulting black hole mass can be arbitrarily small if one is fine-tuning the amount of kinetic energy added to the initial neutron star. This result is in accordance with works by Choptuik (1993) on critical collapses, more precisely the mass-scaling relation (14) also applies here, at least for black hole masses not too close to the neutron star critical mass. The mass scaling exponent $`\alpha `$ found for the polytropic EOS is rather different from that found, with quite different initial conditions, by Neilsen and Choptuik (2000) for ultrarelativistic $`\gamma =2`$ EOS.

Nevertheless, the velocity necessary to achieve the collapse to a black hole of a typical $`1.4\mathrm{M}_{}`$ neutron star is enormous (see Fig. 5). Even in the case of a smaller $`K`$, leading to a smaller maximal mass (see Sec. 2.3), this velocity always remains larger than $`0.1c`$. It is therefore difficult to imagine any physical process which could inject such an amount of kinetic energy to a neutron star. Tidal (anisotropic) effects when the neutron star is in a binary system with another compact object may result in the destruction of the star (compressed to a “pancake”, as ordinary stars passing near a black hole in Marck et al. (1996)) or in its oscillations. It has been argued (e.g. by Mathews and Wilson (2000)) that in a binary neutron star system, tidal effects could induce such a compression (i.e. increase of central density) that the stars could collapse to black holes before their merging. This is very difficult to achieve, since in our calculations it happened that during the collapse the central density could reach a value substantially higher than the critical density $`n_B^{\mathrm{crit}}`$ (in some cases several times $`n_B^{\mathrm{crit}}`$), whereas the star would not end in a black hole, but would rather be dispersed or enter an infinite series of oscillations. Even in the case of the neutron star instability transition in tensor-scalar theory (see Novak (1998)), the velocity is never larger than a few percents of $`c`$, therefore quite far from values displayed on Fig. 5. Finally, let us mention the recent interesting work by Hawley and Choptuik (2000) who studied the critical collapses of boson stars, where the interaction with a massless real scalar field can result in a significant transfer of energy from the field to the star, allowing for the collapse to a black hole of stable configurations. In the case of neutron stars the interaction with a gravitational wave seems unlikely to transfer enough energy to lead to the collapse to a black hole.

###### Acknowledgements.

I am indebted to José A. Pons for providing with the realistic equation of state. I also would like to thank Eric Gourgoulhon for his help and very fruitful discussions, José M$`^{\underset{¯}{\text{a}}}`$ Ibáñez who suggested this study and the anonymous referee for his comments.