text:
# Distilling Gaussian states with Gaussian operations is impossible

## Abstract

We will show that no distillation protocol for Gaussian quantum states exists that relies on (i) arbitrary local unitary operations that preserve the Gaussian character of the state and (ii) homodyne detection together with classical communication and postprocessing by means of local Gaussian unitary operations on two symmetric identically prepared copies. This analysis shows that unlike the finite-dimensional case, where entanglement can be distilled in an iterative protocol using two copies at a time, there is no such procedure in the case of continuous variables for Gaussian initial states and the above Gaussian operations. The ramifications for the distribution of Gaussian states over large distances will be outlined. We will also comment on the generality of the approach and sketch the most general form of a Gaussian local operation with classical communication in a bi-partite setting.

In most practical implementations of information processing devices sophisticated methods are necessary in order to preserve the coherence of the involved quantum states. Even the mere preparation of an entangled state of spatially distributed quantum systems requires such techniques: once prepared locally and then distributed, an entangled state will to some extent deteriorate from a highly entangled state to a less correlated state through the process of decoherence. This process can quite obviously not be entirely avoided. However, one may prepare and distribute several identical entangled states, and then apply appropriate partly measuring local quantum operations and classical communication to obtain states that are similar to the highly entangled original state. This is only possible at the expense that one has fewer identically prepared systems or copies at hand, but this is a small price to pay. Appropriately indeed, this process has been given the name distillation , as fewer more highly entangled states are ’distilled’ from a supply of many less entangled states. It was one of the major successes of the field of quantum information science to realize that for two-level systems such a distillation procedure may be performed on only two copies at a time, and it requires only two steps: (i) a local unitary operation, and (ii) a local measurement, together with the classical communication about the measurement outcome. Based on the measurement outcome further local unitary operations are then implemented.

Such distillation protocols may also be of crucial importance in the infinite-dimensional setting. Quantum information science over continuous variables has seen an enormous progress recently, both in theory and experiment, mostly involving Gaussian states of field modes in a quantum optical setting . Quite naturally, one should expect that a similar distillation procedure also works for Gaussian states in the infinite-dimensional case, also under the preservation of the Gaussian character of the state. If one transmits two pure two-mode squeezed Gaussian states through lossy optical systems such as fibers, the corresponding modes being from now on labeled $`A1`$, $`A2`$, $`B1`$, and $`B2`$, one obtains two identical copies of less entangled symmetric states . A feasible distillation protocol preserving the Gaussian character may consist of the subsequent steps (see Fig. 1):

(i) Application of any local Gaussian unitary operation. That is, one may implement any unitary operations $`U_A`$ and $`U_B`$ on both $`A1`$ and $`A2`$ on one hand and $`B1`$ and $`B2`$ on the other hand corresponding to symplectic transformations  $`S_A,S_BSp(4,)`$ . This set includes all two-mode and one-mode squeezings, mixing at beam splitters and phase shifts. To specify these operations 20 real parameters are necessary. Note that we do not require both parties to realize the same transformation.

(ii) A homodyne measurement on the modes $`A2`$ and $`B2`$. The parties communicate classically about the outcome of the measurement, and may postprocess the states of modes $`A1`$ and $`B1`$ with unitary Gaussian operations.

The main result of this Letter is that very much as a surprise, none of these protocols amounts to a distillation protocol. No matter how ingeniusely the local unitary operation is chosen, the degree of entanglement can not be increased. The optimal procedure is simply to do nothing at all, which means that at least no entanglement is lost . The degree of entanglement will be measured in terms of the log-negativity, which is defined as $`E_N(\rho )=\mathrm{log}_2\rho ^{T_A}`$ for a state $`\rho `$, where $`.`$ denotes the trace norm, and $`\rho ^{T_A}`$ is the partial transpose of $`\rho `$. The negativity has been shown to be an entanglement measure in the sense that it is non-increasing on average under local operations with classical communication , and is to date the only known feasible measure of entanglement for Gaussian states. For pure (and for symmetric mixed) Gaussian states it is related to the degree of squeezing in a monotone way (see, e.g., ). This means that as a corollary of the main result, it follows that with Gaussian operations as specified above one cannot transform two identically prepared two-mode squeezed vacua into a single two-mode squeezed vacuum state with a higher degree of squeezing.

We will start by fixing the notation. Gaussian states  of an $`n`$-mode system are completely characterized by their first and second moments. The first moments are the expectation values of the canonical coordinates. The second moments can be collected in the real symmetric covariance matrix $`\mathrm{\Gamma }C(2n)M(2n,)`$, where $`M(2n,)`$ denotes the set of real $`2n\times 2n`$-matrices, and $`C(2n)`$ the subset of matrices obeying the Heisenberg uncertainty principle . The linear transformations from one set of canonical coordinates to another which preserve the canonical commutation relations form the group of real linear symplectic transformations $`Sp(2n,)`$ . A symplectic transformation $`S`$ changes the covariance matrix according to $`\mathrm{\Gamma }S\mathrm{\Gamma }S^T`$, while states undergo a unitary operation $`\rho U\rho U^{}`$. The $`n=4`$ modes $`A1`$, $`A2`$, $`B1`$, $`B2`$ will be equipped with the canonical operators $`(X_{A1},P_{A1},\mathrm{},X_{B2},P_{B2})`$. To make the notation more transparent, both tensor products and direct sums will carry a label indicating the underlying split, meaning either $`A,B`$ or $`1,2`$. We state the main result of this Letter in form of a theorem:

Theorem. – Let $`\rho \rho `$ be two identically prepared symmetric Gaussian states of two-mode systems consisting of the parts $`A1`$, $`A2`$, $`B1`$, and $`B2`$, respectively, each of which having the covariance matrix

$$\mathrm{\Gamma }^{(0)}=\left(\begin{array}{cccc}a& 0& c& 0\\ 0& a& 0& c\\ c& 0& a& 0\\ 0& c& 0& a\end{array}\right),a1,\mathrm{\hspace{0.17em}\hspace{0.17em}0}c(a^21)^{1/2},$$
(1)
let $`S_A,S_BSp(4,)`$ be any symplectic transformations with associated unitaries $`U_A`$ and $`U_B`$, and let

$$\rho ^{}=(U_A_{A,B}U_B)(\rho _{1,2}\rho )(U_A_{A,B}U_B)^{}.$$
(2)
Then any state $`\rho ^{\prime \prime }`$ that is obtained from $`\rho ^{}`$ via a selective homodyne measurement on systems $`A2`$ and $`B2`$ satisfies $`E_N(\rho ^{\prime \prime })E_N(\rho )`$, that is, the degree of entanglement can only decrease.

The proof of this statement will turn out to be technically involved, and while the statement itself is concerned with practical quantum optics, the techniques used in the proof will be mostly taken from matrix analysis . In order to give the general argument more structure, the proof is split into several lemmata. The entire proof will be formulated in terms of covariance matrices, rather than in terms of the states.

The log-negativity of a state $`\sigma `$ of a two-mode system can be easily expressed in terms of the entries of the associated covariance matrix $`\gamma C(4)`$. The latter can be partitioned in block form according to

$$\gamma =\left(\begin{array}{cc}\gamma _A& \gamma _C\\ \gamma _C^T& \gamma _B\end{array}\right),\gamma _A,\gamma _B,\gamma _CM(2,R).$$
(3)
The log-negativity $`E_N(\sigma )`$ is then given by

$`E_N(\sigma )=\{\begin{array}{cc}(\mathrm{log}f)(\gamma )/2,& \text{if }f(\gamma )<1,\\ 0& \text{otherwise.}\end{array}`$ (6)
where the function $`f:C(4)^+`$ is defined as

$`f(\gamma )`$ $`=`$ $`((\text{det}[\gamma _A]+\text{det}[\gamma _B])/2\text{det}[\gamma _C])`$ (7)
$``$ $`\left[((\text{det}[\gamma _A]+\text{det}[\gamma _B])/2\text{det}[\gamma _C])^2\text{det}[\gamma ]\right]^{1/2}.`$ (8)
The covariance matrix associated with the Gaussian state $`\rho ^{}`$ in the Theorem will be denoted as $`\mathrm{\Gamma }^{}C(8)`$. For any $`S_A,S_BSp(4,)`$ this covariance matrix of the modes $`A1`$, $`A2`$, $`B1`$, and $`B2`$ becomes $`\mathrm{\Gamma }^{}:=(S_A_{A,B}S_B)(\mathrm{\Gamma }^{(0)}_{1,2}\mathrm{\Gamma }^{(0)})(S_A_{A,B}S_B)^T`$ The first step is to relate the covariance matrix $`\mathrm{\Gamma }^{\prime \prime }`$ associated with the state after the measurement to a Schur complement . This Schur complement structure is a general feature of Gaussian operations and will be further discussed at the end of the letter.

Lemma 1. – Let $`\mathrm{\Gamma }^{}C(8)`$ be a covariance matrix of systems $`A1`$, $`A2`$, $`B1`$, and $`B2`$ associated with a state $`\rho ^{}`$, which can be written in block form as

$$\mathrm{\Gamma }^{}=\left(\begin{array}{cc}C_1& C_3\\ C_3^T& C_2\end{array}\right),$$
(9)
where $`C_1,C_2,C_3M(4,)`$. The covariance matrix of the state that is obtained by a projection in $`A2`$ and $`B2`$ on the pure Gaussian state with covariance matrix $`D_d:=\text{diag}(1/d,d,1/d,d)C(4)`$, $`d>0`$, is then given by

$$M_d=C_1C_3(C_2+D_d^2)^1C_3^T.$$
(10)
Proof. This statement can be most conveniently be shown in terms of the characteristic function $`\chi `$ . By employing the Weyl (displacement) operator, the state $`\rho ^{}`$ associated with the covariance matrix $`\mathrm{\Gamma }^{}`$ can be written in terms of the characteristic function according to $`\rho ^{}=(1/\pi ^4)d^8\xi W(\xi )\chi (\xi )`$ (see, e.g., Ref. ). The projection corresponds on the level of the characteristic function therefore to an incomplete Gaussian integration. The characteristic function associated with the modes $`A1`$ and $`B1`$ can then be written as

$`\chi (\xi _1,\mathrm{},\xi _4)`$ $`=`$ $`{\displaystyle \frac{d\xi _5\mathrm{}d\xi _8}{\pi ^2}e^{\xi ^T\mathrm{\Gamma }^{}\xi /2}e^{\frac{1}{2d^2}(\xi _5^2+\xi _7^2)\frac{d^2}{2}(\xi _6^2+\xi _8^2)}}`$ (11)
$`=`$ $`|C_2+D_d^2|^{1/2}e^{(\xi _1,\mathrm{},\xi _4)\mathrm{\Gamma }^{\prime \prime }(\xi _1,\mathrm{},\xi _4)^T/2},`$ (12)
with $`M_d`$ defined as in Eq. (10).

Hence, the resulting covariance matrix is given by the Schur complement $`C_1C_3(C_2+D_d^2)^1C_3^T`$ of the matrix

$$\mathrm{\Gamma }_d^{}:=\left(\begin{array}{cc}C_1& C_3\\ C_3^T& C_2+D_d^2\end{array}\right)$$
(13)
with respect to the leading principal submatrix $`C_1`$. The additional matrix $`D_d^2`$ originates from the projection in the modes $`A2`$ and $`B2`$. Note that although this Lemma has been formulated in terms of the projection on a certain class of pure Gaussian states, it applies to the projection on any pure Gaussian state in the modes $`A2`$ and $`B2`$: the projection on any other pure Gaussian state can be realized by an appropriate choice of the symplectic transformations $`S_A`$ and $`S_B`$. Ideal homodyne detections can now be formulated as projections on ‘infinitely squeezed’ pure Gaussian states . The central feature is that the initial first moments do not affect the form of the covariance matrix after the measurement. Lemma 2 gives the form of the resulting covariance matrix in case of a homodyne detection in modes $`A2`$ and $`B2`$. In the limit $`d0`$ the matrix $`D_d`$ gives rise to a projection operator, and the inverse becomes a Moore Penrose inverse (MP) :

Lemma 2. – In the notation of Lemma 1, the covariance matrix of modes $`A1`$ and $`B1`$ after a selective homodyne measurement in modes $`A2`$ and $`B2`$ is given by

$$\mathrm{\Gamma }^{\prime \prime }:=\underset{d0}{lim}\mathrm{\Gamma }_d=C_1C_3(\pi C_2\pi )^{\text{MP}}C_3^T,$$
(14)
where $`\pi =\text{diag}(1,0,1,0)`$.

Equipped with these preparatory considerations, we will now turn to the core of the proof. In order to be able to evaluate the logarithmic negativity according to Eq. (6), one needs to know the values of the invariants under local symplectic transformations, i.e., the determinants of four submatrices. To find an expression for all these determinants is however a quite difficult task. Instead, we will later make use of an upper bound of the logarithmic negativity that only involves determinants of principal submatrices  of $`\mathrm{\Gamma }^{\prime \prime }`$.

Lemma 3. – Let $`\mathrm{\Gamma }^{\prime \prime }C(4)`$ be defined as in Lemma 2. Then, independent of $`S_A,S_BSp(4,)`$,

$$\text{det}[\mathrm{\Gamma }^{\prime \prime }]=\text{det}[\mathrm{\Gamma }^{(0)}]=(a^2c^2)^2.$$
(15)
Proof. According to Lemma 2, $`\mathrm{\Gamma }^{\prime \prime }`$ is given by $`\mathrm{\Gamma }^{\prime \prime }=lim_{d0}M_d`$. The Schur complement of the matrix $`\mathrm{\Gamma }_d^{}`$ as defined in Eq. (13) is related to $`\mathrm{\Gamma }_d^{}`$ and one of its principal submatrices via the congruence

$$\left(\begin{array}{cc}𝟙_\mathrm{𝟜}& X\\ 0& 𝟙_\mathrm{𝟜}\end{array}\right)\mathrm{\Gamma }_d^{}\left(\begin{array}{cc}𝟙_\mathrm{𝟜}& 0\\ X^T& 𝟙_\mathrm{𝟜}\end{array}\right)=\left(\begin{array}{cc}\mathrm{\Gamma }_d& 0\\ 0& C_2+D_d^2\end{array}\right),$$
(16)
where $`X:=C_3(C_2+D_d^2)^1`$. Hence, according to the determinant multiplication theorem we obtain $`\text{det}[\mathrm{\Gamma }_d^{}]=\text{det}[\mathrm{\Gamma }_d]\text{det}[C_2+D_d^2]`$, which yields in the limit $`d0`$

$$\text{det}[\mathrm{\Gamma }^{\prime \prime }]=\text{det}[P\mathrm{\Gamma }_d^{}P+(𝟙_\mathrm{𝟠})]/\text{det}[\mathbb{\Gamma }_𝕕^{}+(𝟙_\mathrm{𝟠})],$$
(17)
where the projections $`P`$ and $`Q`$ are defined as $`P:=\text{diag}(1,1,1,1,0,1,0,1)`$ and $`Q:=\text{diag}(0,0,0,0,0,1,0,1)`$. With these tools, it is feasible to directly prove the statement of Lemma 3 by parameterizing $`S_A,S_BSp(4,)`$. Every $`SSp(4,)`$ can be written as a product $`S=VDW`$, where $`V,WSp(4,)SO(4)`$, and $`D:=\text{diag}(d_1,1/d_1,d_2,1/d_2)`$ with $`d_1,d_2`$ .

Lemma 4. – Let $`\mathrm{\Gamma }^{\prime \prime }C(4)`$ be defined as in Lemma 2, and let $`\mathrm{\Gamma }_A^{\prime \prime }`$ and $`\mathrm{\Gamma }_B^{\prime \prime }`$ be the principal submatrices belonging to mode $`A1`$ and $`B1`$. Then, for all $`S_A,S_BSp(4,)`$,

$`\text{det}[\mathrm{\Gamma }_A^{\prime \prime }]\text{det}[\mathrm{\Gamma }_A^{(0)}]=a^2,\text{det}[\mathrm{\Gamma }_B^{\prime \prime }]\text{det}[\mathrm{\Gamma }_B^{(0)}]=a^2.`$ (18)
Proof. $`\mathrm{\Gamma }^{\prime \prime }`$ is defined as the covariance matrix corresponding to modes $`A1`$ and $`B1`$ after the projective measurements in both $`A2`$ and $`B2`$. Let us assume that one first performs the projective measurement in $`A2`$, leading to a the covariance matrix $`N_AC(2)`$ of the reduced state of $`A1`$. The covariance matrix $`\mathrm{\Gamma }_A^{\prime \prime }`$ after the projection in $`B2`$ is then obtained as a Schur complement. In particular, $`\mathrm{\Gamma }_A^{\prime \prime }`$ can be written as $`\mathrm{\Gamma }_A^{\prime \prime }=N_AR`$, where $`RM(2,)`$ is a real symmetric positive matrix. Hence, as $`\mathrm{\Gamma }_A^{\prime \prime }`$ and $`N_A`$ are also positive, $`\text{det}[\mathrm{\Gamma }_A^{\prime \prime }]\text{det}[N_A]`$ . In other words, one obtains an upper bound for $`\text{det}[\mathrm{\Gamma }_A^{\prime \prime }]`$ when considering only a projective measurement in $`A1`$. The statement of Lemma 4 follows from Lemma 3 in the special case that $`c=0`$: one can after a few steps conclude that then $`\text{det}[N_A]=a^2`$, independent of $`S_A,S_BSp(4,)`$. The same reasoning applies to $`\mathrm{\Gamma }_B^{\prime \prime }`$.

The most important step is now an appropriate upper bound of the log-negativity of the resulting state. The actual bound might appear somewhat arbitrary, but it will turn out that it is exactly the tool that we need in the last step of the proof.

Lemma 5. – Let $`\gamma C(4)`$, partitioned as in Eq. (3). Then

$`f(\gamma )g(\gamma )`$ $`:=`$ $`[((\text{det}[\gamma _A]+\text{det}[\gamma _B])/2)^{1/2}`$ (19)
$``$ $`((\text{det}[\gamma _A]+\text{det}[\gamma _B])/2\text{det}[\gamma ]^{1/2})^{1/2}]^2.`$ (20)
Proof. $`g(\gamma )`$ can be expressed in terms of $`f`$ as $`g(\gamma )=f(\gamma ^{})`$,

$$\gamma ^{}=\left(\begin{array}{cc}\gamma _A^{}& \gamma _C^{}\\ \gamma _{C}^{}{}_{}{}^{T}& \gamma _B^{}\end{array}\right),\gamma _A^{}=\gamma _B^{}=a^{}𝟙_\mathrm{𝟚},\gamma _{}^{}=\left(\begin{array}{cc}𝕔^{}& \mathrm{𝟘}\\ \mathrm{𝟘}& 𝕔^{}\end{array}\right),$$
(21)
where $`a^{}:=((\text{det}[\gamma _A]^2+\text{det}[\gamma _B]^2)/2)^{1/2}`$ and $`c^{}:=(a_{}^{}{}_{}{}^{2}\text{det}[\gamma ]^{1/2})^{1/2}`$. Hence, one has to prove that $`f(\gamma ^{})f(\gamma )`$. Firstly, note that $`\text{det}[\gamma ]=\text{det}[\gamma ^{}]`$. Secondly, $`(\text{det}[\gamma _A]+\text{det}[\gamma _B])/2=a_{}^{}{}_{}{}^{2}`$. Therefore, it remains to be shown that $`c_{}^{}{}_{}{}^{2}|\text{det}[\gamma _C]|`$. This inequality is equivalent with $`\left[(\text{det}[\gamma _A]+\text{det}[\gamma _B])/2|\text{det}[\gamma _C]|\right]^2\text{det}[\gamma ]0`$, which is a valid inequality, as $`\gamma C(4)`$.

Proof of the Theorem. Let $`\mathrm{\Gamma }^{\prime \prime }C(4)`$ be the matrix defined as in Lemma 2. The log-negativity of the corresponding state of modes $`A1`$ and $`B1`$ is given by $`(\mathrm{log}f)(\mathrm{\Gamma }^{\prime \prime })`$, if the final state is entangled at all, as we will assume from now on. Lemma 5 yields the bound $`f(\mathrm{\Gamma }^{\prime \prime })g(\mathrm{\Gamma }^{\prime \prime })`$. In $`g(\mathrm{\Gamma }^{\prime \prime })`$, however, only the determinants of the principal submatrices are needed, bounds of which are available by virtue of Lemma 3 and 4. The function $`h:[y,\mathrm{})^+`$ with $`h(x)=(x^{1/2}(xy)^{1/2})^2`$, $`y>0`$, is a strictly monotone decreasing function of $`x`$. Therefore, using Lemma 3 and 4 one can conclude that $`g(\mathrm{\Gamma }^{\prime \prime })g(\mathrm{\Gamma }^{(0)})`$. Moreover, $`g(\mathrm{\Gamma }^{(0)})=f(\mathrm{\Gamma }^{(0)})`$, due to the special form of $`\mathrm{\Gamma }^{(0)}`$, as can be easily verified. Hence, $`f(\mathrm{\Gamma }^{\prime \prime })g(\mathrm{\Gamma }^{(0)})=f(\mathrm{\Gamma }^{(0)})`$, which leads to $`(\mathrm{log}f)(\mathrm{\Gamma }^{\prime \prime })(\mathrm{log}f)(\mathrm{\Gamma }^{(0)})`$. This is finally the desired result: it means that the degree of entanglement can only decrease.

We will finally comment on the generality of the approach. A general Gaussian operation is a quantum operation that maps all Gaussian states on Gaussian states . Any general Gaussian local operation with classical communication (LOCCG) – trace-preserving or non-trace-preserving – can be decomposed into the subsequent steps: (i) Appending locally additional modes that have been prepared in a Gaussian state . (ii) Application of any local unitary Gaussian operation on both the original and the additional system. These comprise operations corresponding to symplectic transformations and displacements in phase space. (iii) Projections on pure Gaussian states or ideal homodyne detections, which give rise to Schur complements on the level of covariance matrices as described above, together with the classical communication about the outcome (real numbers in case of homodyne detection, bits in case of dichotomic measurements including the projection on a pure Gaussian state in one outcome), (iv) mixing, such that the resulting state is Gaussian, and (v) a partial trace, which corresponds to considering a certain principal submatrix of the covariance matrix only . The proof is therefore restrictive in the sense that only two copies at a time are considered, other projections on Gaussian states are excluded, and no additional modes are allowed for. The statement of the present paper proves that iterative protocols in strict analogy to the corresponding methods in finite dimensional settings certainly do not work. Indeed, the findings strongly suggest that Gaussian states cannot be distilled at all with Gaussian operations. Then (less feasible) non-linear physical effects  would have to be made use of in order to distill from a supply of Gaussian two-mode states . Such techniques would then also be necessary for the realistic implementation of quantum repeaters  for continuous-variable systems when it comes to the distribution of highly entangled Gaussian states over large distances.

We would like to thank K. Audenaert, J. Fiurášek, P. van Loock, C. Silberhorn, G. Giedke, J.I. Cirac, S.D. Bartlett, B.C. Sanders, D.-G. Welsch, and N. Cerf for discussions. This work has been supported by the European Union (EQUIP, QUEST) and the A.-v.-Humboldt-Foundation.


format: text

text:
# A Computational Approach to the 𝒟-module of Meromorphic Functions

## 1 Introduction

Let us denote by $`𝒪=𝒪_{𝐂^n}`$ the sheaf of holomorphic functions on $`X=𝐂^n`$. Consider a point $`p𝐂^n`$. $`Der(𝒪_p)`$ is the $`𝒪_p`$-module of $`𝐂`$-derivations of $`𝒪_p`$. The elements in $`Der(𝒪_p)`$ are called vector fields.

Let $`DX`$ be a divisor (i.e. a hypersurface) and $`pD`$. A vector field $`\delta Der(𝒪_p)`$ is said to be logarithmic with respect to $`D`$ if $`\delta (f)=af`$ for some $`a𝒪_p`$, where $`f`$ is a local (reduced) equation of the germ $`(D,p)(𝐂^n,p)`$. The $`𝒪_p`$-module of logarithmic vector fields (or logarithmic derivations) is denoted by $`Der(\mathrm{log}D)_p`$. This yields an $`𝒪`$-module sheaf denoted by $`Der(\mathrm{log}D)`$ (see ).

Let us denote by $`𝒟=𝒟_X`$ the sheaf (of rings) of linear differential operators with holomorphic coefficients on $`X=𝐂^n`$. A local section $`P`$ of $`𝒟`$ (i.e. a linear differential operator) is a finite sum $`P=_\alpha a_\alpha ^\alpha `$ where $`\alpha =(\alpha _1,\mathrm{},\alpha _n)𝐍^n`$, $`a_\alpha `$ is a local section of $`𝒪`$ and $`=(_1,\mathrm{},_n)`$ with $`_i=\frac{}{x_i}`$ in some local chart.

For any divisor $`D𝐂^n`$ we denote by $`𝒪[D]`$ the sheaf of meromorphic functions with poles along $`D`$. It follows from the results of Bernstein-Björk (, ) on the existence of the $`b`$-function for each local equation $`f`$ of $`D`$, that $`𝒪[D]`$ is a left coherent $`𝒟`$-module. Kashiwara proved that the dimension of its characteristic variety is $`n`$ and then that $`𝒪[D]`$ is holonomic, .

We will consider some $`𝒟`$-modules associated to any divisor $`D`$:

* The (left) ideal $`I^{\mathrm{log}D}𝒟`$ generated by the logarithmic vector fields $`Der(\mathrm{log}D)`$.
* The (left) ideal $`\stackrel{~}{I}^{\mathrm{log}D}𝒟`$ generated by the set $`\{\delta +a|\delta I^{\mathrm{log}D}\text{ and }\delta (f)=af\}`$. More generally, the ideals $`\stackrel{~}{I}^{(k)\mathrm{log}D}`$ generated by the set

$$\{\delta +ka|\delta I^{\mathrm{log}D}\text{ and }\delta (f)=af\}$$
* The modules $`M^{\mathrm{log}D}=𝒟/I^{\mathrm{log}D}`$, $`\stackrel{~}{M}^{\mathrm{log}D}=𝒟/\stackrel{~}{I}^{\mathrm{log}D}`$ and more generally $`\stackrel{~}{M}^{(k)\mathrm{log}D}=𝒟/\stackrel{~}{I}^{(k)\mathrm{log}D}`$.

The inclusion $`\stackrel{~}{I}^{(k)\mathrm{log}D}Ann_𝒟(1/f^k)`$ yields to a natural morphism $`\varphi _D^k:\stackrel{~}{M}^{(k)\mathrm{log}D}𝒪[D]`$ defined by $`\varphi _D^k(\overline{P})=P(1/f^k)`$ where $`\overline{P}`$ denotes the class of the operator $`P𝒟`$ modulo $`\stackrel{~}{I}^{(k)\mathrm{log}D}`$. The image of $`\varphi _D^k`$ is $`𝒟\frac{1}{f^k}`$, i.e. the $`𝒟`$-submodule of $`𝒪[D]`$ generated by $`1/f^k`$.

Considering the general ideals $`\stackrel{~}{I}^{(k)\mathrm{log}D}`$ is a suggestion of Prof. Tajima. The point is the well known chain of inclusions

$$𝒟f^1𝒟f^2\mathrm{}𝒟f^k=𝒟f^{k1}=\mathrm{}=𝒪[D],$$

where $`k`$ is least integer root of the $`b`$-function.

We are interested in the germ $`(D,p)(𝐂^n,p)`$ for a fixed point $`pD`$. So we will work in the ring of germs of linear differential operators $`𝒟_p`$. We suppose that $`p=0𝐂^n`$ and from now on, we will denote $`𝒟=𝒟_0`$. In this context we will use the Weyl algebra $`A_n(𝐂)`$ as a subring of $`𝒟`$.

Under a computational point of view the divisor $`D`$ will be defined by a polynomial $`f𝐂[x_1,\mathrm{},x_n]`$. The $`b`$-function of $`f`$ is computable by  and we give a direct method to present $`𝒪[D]`$. If the calculation is intractable –as sometimes happens in the examples– we also present an indirect method to deduce that $`𝒪[D]`$ and the modules $`\stackrel{~}{M}^{(k)\mathrm{log}D}`$ do not coincide. The method is strongly based in the following result of :

###### Theorem 1.1.

The restriction to $`D`$ of the sheaf $`Ext_𝒟^i(𝒪[D],𝒪)`$ is zero for $`i0`$.

More precisely, we prove that, under certain algorithmic conditions, some cohomology groups are not zero. The interest of this second method has been tested in ,  and . We point out how the algorithms presented in  –that only calculate cohomology groups in the algebraic case– could also be useful in some analytic situations.

It is important to underline that our methods manage the analytic case. As the inclusion $`A_n(𝐂)𝒟`$ is flat, the computation of syzygies and free resolutions in the Weyl algebra yields to the analogous computations in $`𝒟`$.

## 2 Comparison algorithms

We propose in this section two methods to compare the logarithmic modules presented above. It is important to remark that the computation of the analytic $`Der(\mathrm{log}D)`$ can be made for a divisor $`D`$ if its local equation is a polynomial $`f𝐂[x_1,\mathrm{},x_n]`$ . Simply compute, using Gröbner basis, a system of generators of $`(f_1,\mathrm{},f_n,f)`$ where $`f_i=\frac{f}{x_i}`$ because the inclusion of the Weyl algebra in $`𝒟`$ is flat.

### 2.1 Direct comparison

The first method is complete but needs the calculation of the $`b`$-function.

Experimental evidences show that if the divisor is not locally Euler homogeneous (i.e. there is no $`\delta Der(\mathrm{log}D)`$ such that $`\delta (f)=f`$) the $`b`$-function is hard to compute. More precisely, the problem seems to be the calculation of $`Ann_{𝒟[s]}(f^s)`$ and the use of certain elimination orders during the calculation of Gröbner basis (here $`𝒟[s]`$ stands for the polynomial ring, with the indeterminate $`s`$ commuting with $`𝒟`$ ). Nevertheless, the method is applicable to all the Euler homogeneous divisors we have considered (and to some not Euler homogeneous).

###### Algorithm 2.1.

INPUT: A local equation $`f=0`$ of a divisor $`D`$;

1. Compute the $`b`$-function of $`f`$. Let $`\alpha _0`$ be the least integer root.
2. Compute the ideal $`Ann_𝒟(1/f^{\alpha _0})`$.
3. Compute a set of generators $`\{𝐬_1,\mathrm{},𝐬_r\}`$ of $`Syz(f_1,\mathrm{},f_n,f)`$. The ideal $`\stackrel{~}{I}^{(\alpha _0)\mathrm{log}D}`$ is generated by the elements

$$𝐬_j\left(\begin{array}{c}_1\\ \mathrm{}\\ _n\\ \alpha _0\end{array}\right)𝒟.$$
4. Compare $`Ann_𝒟(1/f^{\alpha _0})`$ and $`\stackrel{~}{I}^{(\alpha _0)\mathrm{log}D}`$.

OUTPUT: $`𝒪[D]\stackrel{~}{M}^{(\alpha _0)\mathrm{log}D}Ann_𝒟(1/f^{\alpha _0})=\stackrel{~}{I}^{(\alpha _0)\mathrm{log}D}`$.

The correctness of the algorithm is obvious as

$$𝒪[D]𝒟\frac{1}{f^{\alpha _0}}𝒟/Ann_𝒟(1/f^{\alpha _0}).$$

### 2.2 Indirect deduction: a sufficient condition

This second method is an alternative way when you can not obtain the $`b`$-function. In the worst case, it only needs the computation of a free resolution of $`\stackrel{~}{M}^{(\alpha )\mathrm{log}D}`$ for any integer $`\alpha 1`$. More precisely, the algorithm looks for a technical condition in some step of the free resolution. In many examples, it is enough to compute only the first syzygies111Taking into account that computing a complete free resolution can be a problem of great complexity, this option is very interesting..

###### Definition 2.2.

If

$$0𝒟^{r_s}\stackrel{\phi _s}{}\mathrm{}𝒟^{r_2}\stackrel{\phi _2}{}𝒟^{r_1}\stackrel{\phi _1}{}𝒟\stackrel{\pi }{}M0$$

is a free resolution of a module $`M`$, we will say that the Successive Matrices Condition (SMC) holds at level $`i`$ if the two succesive morphisms $`\phi _i,\phi _{i+1}`$ have matrices verifying:

1. The matrix of $`\phi _{i+1}`$ has no part in $`𝒪`$ in some column $`j`$, i.e. all the elements in the $`j`$-th column are in the (left) ideal generated by $`_1,\mathrm{},_n`$.
2. The matrix of $`\phi _i`$ has no constants in the in the $`j`$-th row. That is, for each $`P`$ in the $`j`$-th row $`P(1)`$ is a function $`h`$ with $`h(0)=0`$.

###### Algorithm 2.3.

INPUT: A local equation $`f=0`$ of a divisor $`D`$;

1. Compute a set of generators $`\{𝐬_1,\mathrm{},𝐬_r\}`$ of $`Syz(f_1,\mathrm{},f_n,f)`$. The ideal $`\stackrel{~}{I}^{(\alpha )\mathrm{log}D}`$ is generated by the elements

$$𝐬_j\left(\begin{array}{c}_1\\ \mathrm{}\\ _n\\ \alpha \end{array}\right)𝒟.$$
2. Compute a free resolution of $`M=\stackrel{~}{M}^{(\alpha )\mathrm{log}D}`$

$$0𝒟^{r_s}\stackrel{\phi _s}{}\mathrm{}𝒟^{r_2}\stackrel{\phi _2}{}𝒟^{r_1}\stackrel{\phi _1}{}𝒟\stackrel{\pi }{}M0.$$

OUTPUT: IF SMC holds THEN $`𝒪[D]\stackrel{~}{M}^{(\alpha )\mathrm{log}D}.`$

We need a lemma to justify the algorithm. It explains the role of the SMC. The idea is obtaining an element in $`Ker\phi _{i+1}`$ that is not in $`Im\phi _i`$.

###### Lemma 2.4.

Let $`D`$ be a divisor and

$$0𝒟^{r_s}\stackrel{\phi _s}{}\mathrm{}𝒟^{r_2}\stackrel{\phi _2}{}𝒟^{r_1}\stackrel{\phi _1}{}𝒟\stackrel{\pi }{}\stackrel{~}{M}^{(\alpha )\mathrm{log}D}0()$$

a free resolution of $`\stackrel{~}{M}^{(\alpha )\mathrm{log}D}`$ that verifies SMC at level $`i`$. Then

$$Ext_𝒟^i(\stackrel{~}{M}^{(\alpha )\mathrm{log}D},𝒪)0.$$

###### Proof.

To obtain the $`Ext`$ groups, we have to apply the functor $`Hom_𝒟(,𝒪)`$ to $`()`$. Using that

$$Hom_𝒟(𝒟^r,𝒪)𝒪^r$$

we obtain the complex

$$0𝒪\stackrel{\phi _1^t}{}𝒪^{r_1}\stackrel{\phi _2^t}{}𝒪^{r_2}\mathrm{}\stackrel{\phi _s^t}{}𝒪\stackrel{\phi _s^t}{}𝒪^{r_s}0,$$

where $`\phi _i^t`$ denotes the morphism with matrix the transposed of $`\phi _i`$. The derivatives now act naturally.

Then

$$Ext_𝒟^i(\stackrel{~}{M}^{(\alpha )\mathrm{log}D},𝒪)=Ker\phi _{i+1}^t/Im\phi _i^t.$$

If the matrix of $`\phi _{i+1}`$ has no part in $`𝒪`$ in the $`j`$-th row, then $`𝐞=(0,\mathrm{},1,\mathrm{},0)`$ –where $`1`$ is in the $`j`$-th position– is in $`Ker\phi _{i+1}^t`$, as the derivatives applied to 1 are zero.

This element can not be in $`Im\phi _i^t`$ if the matrix has no constants. Applying the operators of the matrix it is not possible to obtain elements of degree 0. ∎

We have the key to state the main result of this section: the correctness of 2.3.

###### Proposition 2.5.

Let $`D`$ be a divisor with a free resolution of $`\stackrel{~}{M}^{(\alpha )\mathrm{log}D}`$ that verifies SCM at some level. Then $`𝒪[D]\stackrel{~}{M}^{(\alpha )\mathrm{log}D}.`$

###### Proof.

Evident from 2.4 and 1.1. ∎

A special case of SCM appears when you have a resolution of type

$$0𝒟^{r_n}\stackrel{\phi _n}{}\mathrm{}𝒟^{r_2}\stackrel{\phi _2}{}𝒟^{r_1}\stackrel{\phi _1}{}𝒟\stackrel{\pi }{}M0$$

of length $`n`$. Then

$$Ext_𝒟^n(M,𝒪)=\frac{𝒪^{r_n}}{Im\phi _n^t},$$

and SCM means, at level $`n`$, that you can find in the matrix of $`\phi _n`$ a row with no constants.

###### Remark 2.6.

Of course, natural generalizations of the SCM condition has to do with finding explicit elements in some $`Ker\phi _{i+1}^t`$ with special properties. It is not that easy in general! Nevertheless the results of  can be applied in this situation as follows:

$$Ext_{A_n(𝐂)}^0(R[D],R)0Ext_𝒟^0(𝒪[D],𝒪)0,$$

where $`R`$ is the ring of polynomials and $`R[D]`$ is its localization with respect to the equation $`f`$ of $`D`$.

## 3 Application to the Spencer case.

In this section we explain how to apply the sufficient condition to a special case in which a tailored free resolution is provided.

###### Definition 3.1.

() The divisor $`D`$ is said to be free at the point $`pD`$ if the $`𝒪_p`$-module $`Der(\mathrm{log}D)_p`$ is free. The divisor $`D`$ is called free if it is free at each point $`pD`$.

Smooth divisors and normal crossing divisors are free. By  any reduced germ of plane curve $`D𝐂^2`$ is a free divisor.

By Saito’s criterium , $`D(f=0)𝐂^n`$ is free at a point $`p`$ if and only if there exist $`n`$ vector fields $`\delta _i=_{j=1}^na_{ij}_j`$, $`i=1,\mathrm{},n`$, such that $`det(a_{ij})=uf`$ where $`u`$ is a unit in $`𝒪_p`$. Here $`_j`$ is the partial derivative $`\frac{}{x_j}`$ and $`a_{ij}`$ is a holomorphic function in $`𝒪_p`$.

###### Definition 3.2.

We say that a free divisor $`D`$ is of Spencer type if the complex

$$𝒟_𝒪^{}Der(\mathrm{log}D)M^{\mathrm{log}D}0$$

(introduced in ) is a (locally) free resolution of $`M^{\mathrm{log}D}`$ and if this last $`𝒟`$-module is holonomic.

There are analogous resolutions for the family of modules $`\stackrel{~}{M}^{(k)\mathrm{log}D}.`$

For this family of divisors, the solution complex $`Sol(M^{\mathrm{log}D})`$ (that is, the complex $`𝐑om_𝒟(M^{\mathrm{log}D},𝒪)`$) is naturally quasi-isomorphic to $`\mathrm{\Omega }^{}(\mathrm{log}D)`$ (as we pointed in  as a deduction of ). On the other hand, a duality theorem proved in  has important consequences comparing $`\stackrel{~}{M}^{\mathrm{log}D}`$ and $`𝒪[D]`$, namely

###### Theorem 3.3.

() In dimension 2, the morphism $`\varphi _D^1`$ is an isomorphism if and only if $`D`$ is a quasi-homogeneous plane curve.

###### Theorem 3.4.

 Suppose the divisor $`D𝐂^n`$ is free and locally quasi-homogeneous. Then the morphism $`\varphi _D^1`$ is an isomorphism (so, $`\stackrel{~}{M}^{\mathrm{log}D}`$ and $`𝒪[D]`$ are isomorphic as $`𝒟`$-modules).

The methods presented in section 2 give us computational tools to check the comparison between $`\stackrel{~}{M}^{\mathrm{log}D}`$ and $`𝒪[D]`$ .

###### Remark 3.5.

Once you have the duality of , you also have a strategy to study the Logarithmic Comparison Theorem (LCT), that is, the complex $`\mathrm{\Omega }^{}(D)`$ of meromorphic differential forms and the complex $`\mathrm{\Omega }^{}(\mathrm{log}D)`$ are quasi isomorphic. You have to travel round the following chain of isomorphisms:

$$\mathrm{\Omega }^{}(D)DR(𝒪[D])DR(𝒟/Ann_𝒟(1/f))DR(\stackrel{~}{M}^{\mathrm{log}D})$$

$$DR((M^{\mathrm{log}D})^{})Sol(M^{\mathrm{log}D})\mathrm{\Omega }^{}(\mathrm{log}D),$$

where for each coherent $`𝒟`$–module $`M`$ we denote by $`DR(M)`$ the de Rham complex of $`M`$ (see ).

###### Remark 3.6.

There are two interesting experimental suggestions:

* We don’t know examples of free divisors with integer roots of their $`b`$-function less than -1.
* We only know free divisors of Spencer type.

Finally, we have the following result:

###### Proposition 3.7.

Let $`D(f=0)`$ be a Spencer divisor. Let $`\delta _1,\mathrm{},\delta _n`$ be a basis of $`Der(\mathrm{log}D)`$ with $`\delta _i=_{j=1}^na_ij_j`$ for $`1jn`$. If $`_{j=1}^n_j(a_{ij})=0`$, then $`\stackrel{~}{M}^{\mathrm{log}D}`$ and $`𝒪[D]`$ are not isomorphic.

###### Proof.

The last matrix of the Spencer free resolution of $`\stackrel{~}{M}^{\mathrm{log}D}`$ is of a very special type . Its elements (due to the duality formulas of ) are of the form

$$\delta _i+\underset{j=1}{\overset{n}{}}_j(a_{ij}).$$

So, applying lemma 2.4 $`Ext_𝒟^n(\stackrel{~}{M}^{\mathrm{log}D},𝒪)0,`$ so $`\stackrel{~}{M}^{\mathrm{log}D}`$ and $`𝒪[D]`$ are not isomorphic. Thus there is no possible quasi-isomorphism between $`DR(𝒟/Ann_𝒟(1/f))`$ and $`DR(\stackrel{~}{M}^{\mathrm{log}D})`$. ∎

## 4 Examples

In the following examples, the computation of syzygies among polynomials have been made with CoCoA (see ).

The computations of syzygies in the Weyl Algebra, global $`b`$-functions and ideals of type $`Ann_𝒟(1/f^\alpha )`$ have been made with kan/sm1, , that is, using the algorithms of .

### 4.1 Example 1: $`D(x(x^2y^3)(x^2zy^3)=0)`$

We will treat here the divisor $`D𝐂^3`$ whose local equation at $`(0,0,0)`$ is given by $`f=0`$ with

$$f=x(x^2y^3)(x^2zy^3).$$

This divisor is (globally) free and $`\delta _1,\delta _2,\delta _3`$ form a (global) basis of $`Der(logD)`$, where

$$\begin{array}{ccc}\delta _1& =& \frac{3}{2}x_x+y_y\hfill \\ \delta _2& =& (y^3zx^2)_z\hfill \\ \delta _3& =& (\frac{1}{2}xy^2)_x\frac{1}{3}x^2_y+(y^2z^2y^2z)_z,\hfill \end{array}$$

whose coefficients verify that

$$\left|\begin{array}{ccc}\frac{3}{2}x& y& 0\hfill \\ 0& 0& y^3zx^2\hfill \\ \frac{1}{2}xy^2& \frac{1}{3}x^2& y^2z^2y^2z\hfill \end{array}\right|=\frac{1}{2}f.$$

To begin with, we have to follow two steps:

* Step 1: Verify that $`M^{\mathrm{log}D}`$ is holonomic222This computation could be made with . The interest of this question is evident: if $`M^{\mathrm{log}D}`$ is not holonomic, the computation of its dual could not be managed as we do.
* Step 2: Compute a free resolution of $`M^{\mathrm{log}D}`$ with Gröbner basis computation of syzygies. Check if $`D`$ has a free resolution of Spencer type. If this happens then duality holds by .

The example verifies these properties:

1. The module $`Syz(\delta _1,\delta _2,\delta _3)`$ is generated by the syzygies obtained from the commutators $`[\delta _i,\delta _j]`$. We have $`Syz(\delta _1,\delta _2\delta _3)=𝐬_{12},𝐬_{13},𝐬_{23}`$ where

$$\begin{array}{ccc}𝐬_{12}& =& (\delta _2,\delta _13,0)\hfill \\ 𝐬_{13}& =& (\delta _3,0,\delta _12)\hfill \\ 𝐬_{12}& =& (0,\delta _3y^2z,\delta _2).\hfill \end{array}$$
2. On the other hand, the module $`Syz(𝐬_{12},𝐬_{13},𝐬_{23})`$ is generated by the element $`𝐫`$:

$$𝐫=(y^2z^2_z+y^2z_z+\frac{1}{2}xy^2_xy^2z+\frac{1}{3}x^2_y,y^3z_zx^2_z,y_y\frac{3}{2}x_x+5).$$

This is the element required to have the Spencer type resolution so, as we have said, duality holds.

We calculate the $`b`$-function of $`f`$. Its least integer root is -1, so

$$𝒪[D]𝒟1/fAnn_𝒟(1/f).$$

To finish, we check that $`\stackrel{~}{I}^{\mathrm{log}D}=Ann_𝒟(1/f)`$.

### 4.2 Example 2: $`D(x^4+y^4+z^4+x^2y^2z^2=0)𝐂^3`$

This is divisor is not free. A set of generators of $`\stackrel{~}{M}^{\mathrm{log}D}`$ is

$$\begin{array}{ccc}\delta _1& =& (1/8x^2y^3z^2+y)+(1/32x^3y^3z^21/16xyz^41/4xy)_x+\hfill \\ & & +(1/32x^2y^4z^2+1/8x^2z^21/4y^2)_y+(1/4yz)_z,\hfill \\ & & \\ \delta _2& =& (1/16x^4y^4z+z)+(1/64x^5y^4z+1/32x^3y^2z^31/4xz)_x+\hfill \\ & & +(1/64x^4y^5z1/16x^4yz1/4yz)_y+(1/8x^2y^21/4z^2)_z,\hfill \\ & & \\ \delta _3& =& (1/16x^4y^3z^31/8x^2y^5z)+\hfill \\ & & +(1/64x^5y^3z^3+1/32x^3y^5z+1/32x^3yz^5+1/16xy^3z^3)_x+\hfill \\ & & +(1/64x^4y^4z^3+1/32x^2y^6z1/16x^4z^31/8x^2y^2z1/4z^3)_y+\hfill \\ & & +(1/8x^2yz^2+1/4y^3)_z,\hfill \\ & & \\ \delta _4& =& (1/16xy^4z^4+x)+(1/64x^2y^4z^4+1/32y^2z^6+1/8y^2z^21/4x^2)_x+\hfill \\ & & +(1/64xy^5z^41/16xyz^41/4xy)_y+(1/4xz)_z,\hfill \\ & & \\ \delta _5& =& (1/8x^2y^2z+1/4z^3)_x+(1/8xy^2z^21/4x^3)_z,\hfill \\ & & \\ \delta _6& =& (1/8xy^5z^2+1/4x^3y^3)+\hfill \\ & & +(1/32x^2y^5z^21/16x^4y^31/16y^3z^41/8x^2yz^21/4y^3)_x+\hfill \\ & & +(1/32xy^6z^21/16x^3y^4+1/8xy^2z^2+1/4x^3)_y.\hfill \end{array}$$

The free resolution is huge, but anyway computable with kan/sm1. It is of type

$$0𝒟\stackrel{\phi _3}{}𝒟^r\stackrel{\phi _2}{}𝒟^s\stackrel{\phi _1}{}𝒟\stackrel{\pi }{}M0$$

To use proposition 2.4 you can check that there the elements in the last matrix $`\phi _3`$ has no constants. So

$$Ext_𝒟^3(\stackrel{~}{M}^{\mathrm{log}D},𝒪)=𝒪/im(\phi _3)^t0.$$

You have $`𝒪[D]\stackrel{~}{M}^{\mathrm{log}D}`$.

### 4.3 Example 3: $`D(x^4+y^5+xy^4+zx^6=0)𝐂^3`$

In this example we detect that we have, in fact, a (non-Euler homogeneous) product that is a free divisor333By the way, this situation is impossible in dimension 2.. The basis of $`Der(\mathrm{log}D)`$ is

$$\begin{array}{ccc}\delta _1& =& (x^6z^2+5/4x^5yz^22x^4zx^25/4xy)_x+\hfill \\ & & +(5/4x^5yz^2+3/2x^4y^2z^25/2x^3yz1/2x^2y^2z3/4xyy^2)_y+\hfill \\ & & +(2x^3z^22xz)_z,\hfill \\ & & \\ \delta _2& =& (4/5x^4yz8/5x^3y^2z3/4x^2y^3z+25/4x^4z+125/16x^3yzxy^21/4y^3+\hfill \\ & & +125/16xy)_x+(x^3y^2z39/20x^2y^3z9/10xy^4z+125/16x^3yz+75/8x^2y^2z\hfill \\ & & 3/4y^3+1/4x^25/16xy+25/4y^2)_y+(8/5xyz6/5y^2z+25/2xz)_z,\hfill \\ & & \\ \delta _3& =& (3/10x^3y^2z3/8x^2y^3z+25/8x^4z+\hfill \\ & & +125/32x^3yz+8/25x^31/8y^3+125/32xy)_x+(3/8x^2y^3z9/20xy^4z+\hfill \\ & & +125/32x^3yz+75/16x^2y^2z+2/5x^2y1/50xy^2+1/40y^3+1/8x^25/32xy+\hfill \\ & & +25/8y^2)_y+(3/5y^2z+25/4xz+16/25)_z.\hfill \end{array}$$

¿From the monomial $`16/25_z`$ we deduce by the Flow Theorem that $`D`$ is a product. It is of Spencer type and the third matrix can be used to show that $`Ext_𝒟^3(\stackrel{~}{M}^{\mathrm{log}D},𝒪)0`$ so $`𝒪[D]\stackrel{~}{M}^{\mathrm{log}D}`$.

### 4.4 Example 4: $`D((x+y)(xz+y)(x^4+y^5+xy^4)=0)𝐂^3`$

In this last example the divisor $`D𝐂^3`$ has as a local equation at $`(0,0,0)`$ the form

$$f=(x+y)(xz+y)(x^4+y^5+xy^4)=0.$$

The divisor is globally free with $`\delta _1,\delta _2,\delta _3`$ as a global basis of $`Der(logD)`$.

$$\begin{array}{ccc}\delta _1& =& (x^25/4xy)_x+(3/4xyy^2)_y+(1/4xz^2+1/4xz)_z\hfill \\ & & \\ \delta _2& =& (xz+y)_z\hfill \\ & & \\ \delta _3& =& (2x^2y^2+5/2xy^3+1/2y^4+5/2x^37x^2y35/4xy^211x^255/4xy)_x+\hfill \\ & & +(3/2xy^3+3/2y^41/2x^3+2x^2y21/4xy^27y^333/4xy11y^2)_y+\hfill \\ & & +(1/2xy^2z^2+1/2y^3z^21/2xy^2z1/2y^3z1/4xyz^2+7/4xyz+3/2y^2z+\hfill \\ & & +3/4xz^2+1/2x^2+1/2xy+5/2xz+7/2yz1/4y)_z.\hfill \end{array}$$

In this case, $`D`$ is of Spencer type and the third matrix has no constants so, again we have that $`Ext_𝒟^3(\stackrel{~}{M}^{logD},𝒪)0`$.

## 5 Acknowledgements

We thank Prof. Tajima and Prof. David Mond for very helpful ideas and comments.


text:
# Instabilities of infinite matter with effective Skyrme-type interactions

## 1 Introduction

Nucleon-nucleon effective interactions with an explicit density dependence have been largely employed since the seventies for studies of nuclear properties. Once the form of the interaction is chosen, either a zero-range force of Skyrme-type  or a finite-range force of Gogny-type , the parameters are determined by fitting some selected properties of doubly-magic nuclei and symmetric nuclear matter at saturation. Within a Hartree-Fock (HF) scheme these interactions are able to describe quantitatively the properties of stable as well as unstable exotic nuclei. The interactions are thus well controlled around the saturation density $`\rho _0`$ of symmetric nuclear matter, for moderate isospin asymmetries and zero temperature. This type of effective interactions has also been employed to study nuclear matter in conditions of astrophysical interest, for instance neutron matter with a finite proton fraction at densities up to several times $`\rho _0`$ and finite temperatures. This system at such conditions is relevant for studying proto-neutron stars. The effective interactions are thus extrapolated to conditions of density and isospin asymmetry which are not experimentally accessible, and one should then ask for the limits of validity of such an extrapolation.

Of course, the description of nuclear matter consisting of nucleons is no longer valid beyond some value of the density, as other degrees of freedom appear. For instance, strange matter is expected beyond $`34\rho _0`$ . For the present discussion, we shall arbitrarily accept that for densities up to $`4\rho _0`$, only nucleons are needed to describe nuclear matter. Even in this range of densities the effective interactions are not well determined. For instance, different effective interactions can give similar equations of state for symmetric nuclear matter and very different results for pure neutron matter. Indeed, the extreme asymmetry of isospin is not part of the usual input to determine the interaction parameters. In some cases, the determination of the parameters includes a fit to the equation of state of neutron matter calculated with microscopic methods and realistic interactions. The Skyrme interaction RATP  takes into account the variational microscopic results of Ref. . More recently, the set of SLy interactions  have included the variational results of Ref.  among the conditions required to fix the parameters. It is worth noting that these microscopic variational calculations are not exempt from uncertainties since realistic interactions are not completely known, especially their three-body part. Indeed, the same theoretical model  employing different realistic interactions leads to different equations of state, both in symmetric nuclear matter and in pure neutron matter, and sizeable differences appear for densities beyond the saturation value $`\rho _0`$. Although this type of approach gives a precious guide for the determination of the effective interaction parameters, it cannot replace the empirical data which are the natural input for the phenomenological interactions.

The extrapolation of effective interactions can result in an unphysical behavior of nuclear matter . For instance, most Skyrme parametrizations predict that the isospin asymmetry energy $`ϵ_I`$ becomes negative when the density is increased. Consequently, the symmetric system would be unstable beyond the saturation density, preferring a largely asymmetric system made by an excess of either protons or neutrons. Another type of instability refers to the magnetic properties of neutron matter. The possibility of a ferromagnetic transition at high densities has been studied long ago, employing different theoretical ingredients, and the results were contradictory. Currently used Skyrme interactions predict that neutron matter becomes spin polarized at densities between $`1.13.5\rho _0`$ . This transition would have important consequences for the evolution of a proto-neutron star since the mean free path of neutrinos would be zero . However, Gogny-type effective interactions either exclude such a ferromagnetic transition or predict it at very high densities . Relativistic mean-field calculations  predict that such a transition could appear at densities beyond $`4\rho _0`$. Finally, recent Monte Carlo simulations  as well as Brueckner-Hartree-Fock calculations  using modern two- and three-body realistic interactions do exclude such an instability.

The purpose of this work is to analyze whether the presence of instabilities is inherent or not to the Skyrme-type interactions. In Section 2 we explain the restrictions on the Skyrme parameters imposed by the requirement that no instabilities should appear at densities in the range $`(14)\rho _0`$. Furthermore, the sound velocity $`v_s/c`$ in matter must remain smaller than unity, and this adds a new constraint. Section 3 contains the discussion of the results. Conclusions are drawn in Section 4. Useful expressions are given in Appendix A and B.

## 2 Constraints from the Landau parameters

The condition that the spin-unpolarized neutron matter should be the most stable phase at any density was used in Ref.  to constrain the Skyrme parameters. As the potential energy of the fully spin-polarized phase only depends on the combination $`t_2(1+x_2)`$, ferromagnetic collapse is necessarily avoided if this combination is positive. Taking into account current Skyrme parametrizations, the authors of Ref.  concluded that $`5/4<x_2<1`$. However, this condition is not sufficient: the family of SLy parametrizations  imposes this constraint but at densities of about 0.37 fm-3, a ferromagnetism instability is predicted in neutron matter .

Here, we make a systematic use of the stability criteria related to the adimensional Landau parameters $`F_l`$, $`G_l`$, $`F_l^{}`$ and $`G_l^{}`$ in symmetric nuclear matter, and $`F_l^{(n)}`$ and $`G_l^{(n)}`$ in neutron matter. These criteria simply establish that any parameter of multipolarity $`l`$ should be greater than $`(2l+1)`$. Current effective interactions satisfy these conditions at saturation, and we shall investigate if they can be maintained for densities up to $`4\rho _0`$.

The general form of the Skyrme interaction contains ten parameters, eight of them denoted as $`t_i`$, $`x_i`$ (with $`i=0,1,2,3`$), the other two being the power $`\sigma `$ of the density dependence and the strength $`W_0`$ of the zero-range spin-orbit term. These parameters are usually determined by fitting some selected properties of finite nuclei and symmetric nuclear matter. In this study we follow a different procedure, using some accepted nuclear matter values to fix as many interaction parameters as possible and using the stability criteria to put bounds on the remaining parameters.

Our starting point is given by the empirical values of the following quantities in symmetric nuclear matter at the saturation point: density $`\rho _0`$, energy per particle $`ϵ_0`$, effective mass $`m_0^{}`$, compression modulus $`K_0`$, surface energy $`ϵ_S`$, and isospin asymmetry energy $`ϵ_I`$. Imposing that a Skyrme-type interaction reproduces the values of $`\rho _0`$, $`ϵ_0`$, $`m_0^{}`$, and $`K_0`$ completely determines the parameters $`t_0`$, $`t_3`$, $`\sigma `$ and the combination $`T_0[3t_1+(5+4x_2)t_2]/8`$. The explicit expressions are given in Appendix A. With a Skyrme interaction, the surface energy $`ϵ_S`$ depends on the previously fixed parameters, on the combination $`T_S[9t_1(5+4x_2)t_2]/8`$ and on the strength $`W_0`$ of the spin-orbit term . The value of $`W_0`$ is usually determined by fitting the spin-orbit splitting of some selected levels in finite nuclei, and for the present discussion we shall assume the often used value $`W_0=120`$ MeV$``$fm5. Furthermore, the symmetry energy $`ϵ_I`$ can be expressed in terms of the above parameters and of the four combinations $`t_ix_i`$ ($`i=03`$). In summary, from the ten Skyrme parameters, one is kept fixed ($`W_0`$), six parameters or combinations are determined by the six empirical inputs, and three combinations are free. We found convenient to choose the free combinations $`xt_1x_1`$, $`yt_2x_2`$, $`zt_3x_3`$.

From here on, our task is to analyze the possible domains of the $`(x,y,z)`$-space for different densities. For a given value of the density $`\rho `$ we would like to know what is the $`(x,y,z)`$ domain where no instability can occur. Stability implies that any adimensional Landau parameter of multipolarity $`l`$ must be larger that $`(2l+1)`$. Skyrme forces only contain monopolar and dipolar contributions to the particle-hole interaction so that all Landau parameters are zero for $`l>1`$. Thus, we have twelve inequalities, eight coming from symmetric nuclear matter, and four from neutron matter. Explicit expressions of the twelve Landau parameters are given in Appendix B.

Two of the twelve inequalities, however, play a special role. The condition $`F_1>3`$ is trivially satisfied for any density as long as $`0<m^{}(\rho _0)<m`$, which is the case. The Landau parameter $`F_0`$ only depends on $`t_1`$, $`t_3`$, $`\sigma `$ and the combination $`T_0`$, which are determined by the initial inputs. Exploring the inequality $`F_0>1`$ as a function of the density, it is found that it is violated at a density smaller than $`\rho _0`$. It corresponds to a density where the compression modulus becomes negative, which is the spinodal point, i.e, the occurrence of a liquid-gas transition. Note that within a Skyrme interaction framework, the spinodal density is completely determined by the saturation properties $`ϵ_0`$, $`\rho _0`$, $`K_0`$ and $`m_0^{}`$. After excluding the inequalities related to $`F_1`$ and $`F_0`$, we are left with ten inequalities to be satisfied by the yet free combinations $`x,y,z`$.

### 2.1 Symmetric nuclear matter

The parameters $`G_0^{}`$ and $`G_1^{}`$ give two constraints for $`y`$:

$`y`$ $`<`$ $`{\displaystyle \frac{10C_1(\rho )}{3\alpha _1\rho ^{2/3}}}+{\displaystyle \frac{2}{3}}(T_02T_S),`$ (1)
$`y`$ $`>`$ $`{\displaystyle \frac{10C_0(\rho )}{\rho }}+{\displaystyle \frac{2}{3}}(T_02T_S),`$ (2)
where the $`C_i`$ and $`\alpha _i`$ are defined in Appendix A. The parameters $`F_1^{}`$, $`G_1`$ give two constraints on combinations of $`x`$ and $`y`$:

$`x{\displaystyle \frac{3}{5}}y`$ $`>`$ $`{\displaystyle \frac{4C_0(\rho )}{\rho }}+{\displaystyle \frac{4}{15}}(T_02T_S),`$ (3)
$`x+{\displaystyle \frac{3}{5}}y`$ $`<`$ $`{\displaystyle \frac{4C_0(\rho )}{\rho }}{\displaystyle \frac{4}{15}}(T_02T_S).`$ (4)
The remaining Landau parameters $`G_0`$ and $`F_0^{}`$ constrain combinations of all three $`x`$, $`y`$ and $`z`$ unknowns:

$`{\displaystyle \frac{1}{9\alpha _1}}(\rho ^\sigma \rho _0^\sigma )z+(\rho ^{2/3}\rho _0^{2/3})x+{\displaystyle \frac{3}{5}}(\rho ^{2/3}+\rho _0^{2/3})y>`$ (5)
$`{\displaystyle \frac{4}{3\alpha _1}}\left(C_1(\rho )+C_1(\rho _0)+{\displaystyle \frac{2ϵ_I}{\rho _0}}\right){\displaystyle \frac{4}{15}}(\rho ^{2/3}+\rho _0^{2/3})(T_02T_S),`$
$`{\displaystyle \frac{1}{9\alpha _1}}(\rho ^\sigma \rho _0^\sigma )z(\rho ^{2/3}\rho _0^{2/3})x+{\displaystyle \frac{3}{5}}(\rho ^{2/3}\rho _0^{2/3})y>`$ (6)
$`{\displaystyle \frac{4}{3\alpha _1}}\left(C_1(\rho )C_1(\rho _0){\displaystyle \frac{2ϵ_I}{\rho _0}}\right){\displaystyle \frac{4}{15}}(\rho ^{2/3}\rho _0^{2/3})(T_02T_S).`$
However, it is possible to take the linear combination $`F_0^{}+G_0`$ to get another constraint independent of $`z`$. Note that Eq. (6) reduces to the trivial condition $`ϵ_I>0`$ when $`\rho =\rho _0`$. It is due to the fact that the symmetry energy $`ϵ_I`$ has been used to write $`x_0`$ in terms of the remaining parameters.

### 2.2 Neutron matter

The parameters $`F_1^{(n)}`$, and $`G_1^{(n)}`$ give two constraints on combinations of $`x`$ and $`y`$:

$`x{\displaystyle \frac{3}{5}}y`$ $`<`$ $`{\displaystyle \frac{4C_0(\rho )}{\rho }}+{\displaystyle \frac{4}{15}}(T_02T_S),`$ (7)
$`x{\displaystyle \frac{1}{5}}y`$ $`<`$ $`{\displaystyle \frac{2C_0(\rho )}{\rho }}{\displaystyle \frac{2}{15}}(T_02T_S).`$ (8)
The remaining Landau parameters $`F_0^{(n)}`$, and $`G_0^{(n)}`$ constrain combinations of all three $`x,y,z`$ parameters:

$`{\displaystyle \frac{(\sigma +1)(\sigma +2)\rho ^\sigma 2\rho _0^\sigma }{6(4\alpha _2\rho ^{2/3}3\alpha _1\rho _0^{2/3})}}z+x{\displaystyle \frac{3}{5}}y<`$ (9)
$`{\displaystyle \frac{4}{4\alpha _2\rho ^{2/3}3\alpha _1\rho _0^{2/3}}}\left(C_1(\rho _0)C_2(\rho )+{\displaystyle \frac{2ϵ_I}{\rho _0}}\right)+{\displaystyle \frac{4}{15}}(T_02T_S),`$
$`{\displaystyle \frac{1}{9}}(\rho ^\sigma \rho _0^\sigma )z+{\displaystyle \frac{1}{3}}(2\alpha _2\rho ^{2/3}3\alpha _1\rho _0^{2/3})x+{\displaystyle \frac{1}{5}}(2\alpha _2\rho ^{2/3}+3\alpha _1\rho _0^{2/3})y>`$ (10)
$`{\displaystyle \frac{4}{3}}\left(C_1(\rho _0)+C_3(\rho )+{\displaystyle \frac{2ϵ_I}{\rho _0}}\right){\displaystyle \frac{4}{45}}(7\alpha _2\rho ^{2/3}+3\alpha _1\rho _0^{2/3})(T_02T_S).`$
Note that the linear combination $`\alpha _1F_0^{}+G_0^{(n)}`$ gives another constraint independent of $`z`$.

### 2.3 Sound velocity constraint

In addition to the constraints from stability requirements it is important to check that for each density the sound velocity $`v_s`$ remains smaller than the speed of light, i.e., superluminosity does not occur. The sound velocity is directly related to the compression modulus $`K(\rho )`$ which can itself be expressed in terms of the Landau parameters $`F_0`$ and $`F_1`$:

$`mv_s^2`$ $`=`$ $`{\displaystyle \frac{1}{9}}K`$
$`=`$ $`{\displaystyle \frac{\mathrm{}^2k_F^2}{3m}}{\displaystyle \frac{1+F_0}{1+\frac{1}{3}F_1}}.`$
In symmetric nuclear matter it turns out that $`v_s`$ depends only on the parameters $`t_0`$, $`t_3`$, $`\sigma `$ and the combination $`T_0`$, all of which being already determined by the six initial empirical inputs. Therefore, there will be no additional bounds on the allowed volume in $`(x,y,z)`$-space. However, Eq. (2.3) shows that $`v_s`$ increases with increasing density and for a given choice of initial empirical inputs there is always a value $`\rho _{sound}`$ beyond which $`v_s/c`$ is greater than 1. The value of $`\rho _{sound}`$ depends essentially on the adopted value for $`K_0`$. For instance, $`\rho _{sound}`$ is around $`3.5\rho _0`$ if $`K_0`$ = 350 MeV and it may become $`6\rho _0`$ if $`K_0`$ = 250 MeV. Thus, there is no Skyrme interaction which can be reasonably be used beyond $`\rho _{sound}`$ because it would predict unphysical values of the sound velocity.

The situation is somewhat different in pure neutron matter. Now, the Landau parameters in Eq. (2.3) depend also on the $`(x,y,z)`$ parameter combinations. The requirement that $`v_s/c`$ remains smaller than unity adds one more constraint to the determination of the allowed volume in the $`(x,y,z)`$ space. We shall see in the next section that the sound velocity in neutron matter does not bring any effective restriction on the allowed volume in the parameter space, at least around the values of the six empirical inputs used in the present study.

We must mention that a study of the sound velocity in asymmetric nuclear matter was made in Ref. . The aim was to find the conditions to be satisfied by the Skyrme parameters so that superluminosity would never appear at any density whatsoever. The conditions which were found are very restrictive. In the present work, our point of view is quite different since we don’t think that the Skyrme effective approach should be valid at densities beyond $`4\rho _0`$.

## 3 Results and discussion

The inputs used in the present study are: $`ϵ_0=16.0`$ MeV, $`\rho _0=0.16`$ fm-3, $`K_0=230`$ MeV, $`m_0^{}/m=0.70`$, $`ϵ_S=18.0`$ MeV, and $`ϵ_I=32.0`$ MeV.

For a given value of the density $`\rho `$, the various constraints define in the ($`x,y,z`$) parameter space an allowed volume that we call $`\mathrm{\Omega }`$. The surface of $`\mathrm{\Omega }(\rho )`$ is made of planes because all constraints are linear in $`x,y,z`$. Any point outside $`\mathrm{\Omega }`$ is forbidden because some of the inequalities would not be fulfilled. To present the results and facilitate the discussion we will explore the $`z`$-axis. The intersection of a volume $`\mathrm{\Omega }(\rho )`$ with a $`z`$ = constant horizontal plane gives a polygon which can be represented in the ($`x,y`$)-plane. By varying $`\rho `$ in the range $`\rho _04\rho _0`$ one can follow the evolution of the polygons. The vanishing of the area of the polygon at some critical density $`\rho _{cr}`$ indicates that there is no Skyrme interaction (having the chosen value of $`z`$) which can fulfill the chosen set of constraints at densities beyond $`\rho _{cr}`$.

We have performed two types of calculations: a) using only the six constraints from the Landau parameters of symmetric nuclear matter and disregarding those of neutron matter; b) using all eleven constraints from symmetric nuclear matter and neutron matter plus the sound velocity constraint in neutron matter. In Fig. 1 are displayed some typical results. We have chosen a positive and a negative value of $`z`$ since there is no a priori limitation on $`z`$. The value $`z=210^4`$ MeV fm3(1+σ) correspond to Figs. 1a) and 1b), and $`z=210^4`$ to Figs. 1a’) and 1b’). Results involving only Landau parameters of symmetric matter are displayed in Figs. 1a) and 1a’), whereas Figs. 1b) and 1b’) correspond to results using all eleven constraints. The contours are drawn every 0.5$`\rho _0`$ starting from $`\rho _0`$.

Let us first examine the results without neutron matter constraints (Figs. 1a and 1a’). At $`\rho =\rho _0`$ the contour is a polygon whose upper and lower horizontal sides are determined by the constraints on $`G_0^{}`$ and $`G_0`$, respectively, whereas the left and right sides correspond to $`F_1^{}`$ and $`G_1`$. When $`\rho `$ increases the lower side becomes tilted and the constraint on $`F_0^{}`$ gradually appears as a new side on the right hand side of the polygon. At $`\rho 2\rho _0`$ and above the $`F_0^{}`$ constraint is dominating over the $`G_1`$ constraint. The surface of the polygon shrinks as $`\rho `$ increases. The value of $`\rho _{cr}`$ is above $`4.5\rho _0`$. This value is reached for $`x`$ and $`y`$ both negative if $`z`$ is positive, whereas for $`z`$ negative $`\rho _{cr}`$ corresponds to $`x0`$ and $`y0`$.

The situation changes somewhat when neutron matter constraints are added (Figs. 1b and 1b’). At $`\rho =\rho _0`$ the lower side is given by the $`G_0^{(n)}`$ constraint while the right side of the polygon is mainly determined by the $`F_0^{(n)}`$ constraint which limits severely the allowed area. It must be noted that the sound velocity does not bring any limitation on the ($`x,y,z`$) parameters in the domain explored here. One can see that the critical density $`\rho _{cr}`$ has now a lower value as compared to the case of Figs. 1a-1a’.

There is no natural limitation to the domain of the $`z`$ parameter and therefore, we must explore the dependence of the results on $`z`$. In Fig. 2 we present calculations where all eleven constraints are included, exploring the $`z`$-axis from negative to positive values. The interpretation of the different sides delimitating the contours is the same as in the case of Fig. 1. One can see that $`\rho _{cr}`$ is highest when $`z`$ is between 0 and $`310^4`$. For negative values of $`z`$ the allowed area at $`\rho =2.5\rho _0`$ is already fairly small. Below $`z310^4`$ the polygons of larger $`\rho `$ are contained inside those of smaller $`\rho `$, i.e., interactions stable at density $`\rho `$ are also stable at all densities between $`\rho _0`$ and $`\rho `$. For $`z`$ larger than $`310^4`$ the allowed areas at larger $`\rho `$ tend to move outside the area of $`\rho =\rho _0`$, which means that that such interactions might seem acceptable at large densities but they have the defect of having instabilities at normal density.

So far we have discussed the results obtained with the values of empirical inputs as adopted in Section 2. It is interesting to see how the conclusions may depend on this choice. In Fig. 3 we study, for a fixed value of $`\rho =2\rho _0`$ and a chosen $`z=10^4`$ MeV.fm3(1+σ), the evolution of the allowed polygons when one changes one of the empirical inputs: the effective mass $`m_0^{}/m`$ (Fig. 3a), the incompressibility $`K_0`$ (Fig. 3b), the surface energy $`ϵ_S`$ (Fig. 3c). It can be seen that the results are sensitive to the value of $`m_0^{}/m`$. A smaller value of $`m_0^{}/m`$ makes the allowed domain larger and therefore, the critical density $`\rho _{cr}`$ is shifted to higher values. On the contrary, the allowed domain is reduced when higher values of $`m_0^{}/m`$ are considered, and instabilities appear for densities below $`2\rho _0`$. Hence, to avoid low density instabilities, one should prefer low $`m_0^{}/m`$. The sensitivity to $`ϵ_S`$ is also non negligible, whereas the dependence on $`K_0`$ is moderate.

In Fig. 4 we present the calculated values of the critical density $`\rho _{cr}`$ (in units of $`\rho _0`$). Fig. 4a shows $`\rho _{cr}`$ as a function of the parameter $`z`$ for various choices of the empirical input $`m_0^{}/m`$ around the value 0.70 . As pointed out above, one can see clearly that it is more favorable to choose $`z`$ positive and less than $`310^4`$ in order to have $`\rho _{cr}`$ not too small, and that larger values of $`m_0^{}/m`$ tend to yield lower $`\rho _{cr}`$. The three curves of Fig. 4b represent $`\rho _{cr}`$ as a function of the relative variations (in percentage) around the standard choice of $`m_0^{}/m`$, $`K_0`$ and $`ϵ_S`$. The value of $`z`$ is fixed at $`10^4`$ MeV.fm3(1+σ). These curves show that $`\rho _{cr}`$ depends little on $`K_0`$ whereas moderate increases of $`m_0^{}/m`$ (from 0.70 to 0.75), or of $`ϵ_S`$ (from 18 MeV to 20 MeV) can lower $`\rho _{cr}`$ below 3$`\rho _0`$.

Now, the question arises whether these bounds for Skyrme parameters are useful or not for calculations in pure nucleon matter and finite nuclei. In principle, we expect a qualitative positive answer for finite nuclei because the used inputs guarantee that the first terms in a leptodermous expansion of the mass formula are well described. To be more quantitative, we have constructed several parametrizations using the following criteria. We first fix an arbitrary value of $`z`$, going from $`10^4`$ to $`410^4`$, and then we take the values of $`x`$ and $`y`$ as the coordinates of the point $`\mathrm{\Omega }(\rho _{cr})`$. In this way instabilities are pushed to the highest values of the density. The sets of parameters are dubbed Skzn and they are given in Table 1. Of course, the criteria of highest value of $`\rho _{cr}`$ is purely arbitrary but our intention here is just to explore the resulting parametrizations and not to give the best parameter set. In Fig. 5 are displayed the results for pure neutron matter. One can see that choosing negative values of $`z`$ results in an increasing neutron effective mass for increasing densities. On the other side, values of $`z`$ greater than $`210^4`$ produce too low binding energies. We have then restricted our calculations in finite nuclei to parametrizations Skz0, Skz1 and Skz2. In Table 2 are displayed the ground state binding energies and charge radii of some doubly-magic nuclei together with the results obtained with SIII and SLy4 interactions. The results are reasonably satisfactory. We have noticed that convergence problems occur with Skz0 when looking for the HF solution of 208Pb. To obtain better results one should relax the condition of highest value for the critical density, and of course include some finite nuclei results in the fine tuning of the parameters.

## 4 Conclusion

In this work we have explored the parameter domains of Skyrme-type forces where the stability conditions related to Landau parameters inequalities can be satisfied up to densities of the order of $`4\rho _0`$. Stability of both nuclear matter and pure neutron matter are considered. We have taken advantage of the possibility to characterize the domains by analytical expressions. Starting from a general Skyrme force containing 10 parameters, we have shown that 7 parameters or combinations thereof can be approximately fixed by physical quantities which can be considered experimentally known to some extent. Thus, the problem reduces to the study of allowed domains in a 3-dimensional space spanned by 3 parameter combinations.

The results show that, for any Skyrme-type interaction there is a critical density $`\rho _{cr}`$ above which one cannot insure all stability conditions. This critical density does not exceed $`3.54\rho _0`$ for a reasonable choice of empirical inputs. It is nevertheless an interesting result to know that it is possible to find Skyrme-type interactions which can give stable nuclear matter and neutron matter up to such densities. The parameter domains are well identified and it would be worthwhile to look inside those domains for Skyrme interactions which can also describe accurately finite nuclei. Our exploration with parametrizations Skzn has shown that reasonable ground state binding energies and charge radii can be obtained. A more systematic study over a wide number of nuclei would be necessary to get the Skyrme-type interaction which is stable over the largest range of densities. This new interaction would be very useful for neutron star calculations, like for instance the neutrino mean free path or the URCA process.

The present analysis is restricted to Skyrme interactions, because the simplicity of the resulting expressions allows for an algebraic study. The presence of instabilities in nuclear matter and neutron matter beyond some critical density could be attributed to the zero-range form of the interaction, and its $`k^2`$-dependence of the effective masses. Finite range effective interactions  are an alternative choice for describing the equation of state of nuclear and neutron matter. We know that the actual parametrizations (D1, D1S, D1P) do not predict instabilities  for densities below $`5\rho _0`$. We have not considered a systematic study of Gogny effective interactions because they contain more parameters, and they do not allow for a simple algebraic analysis.

### Acknowledgments

We have greatly benefited of fruitful discussions with our late colleague Dominique Vautherin on this subject and we acknowledge his invaluable help. This work has been partly supported by DGI (Spain), grant BFM2001-0262. J.N. and N.V.G. thank GANIL for its hospitality during completion of the work.

## Appendix A Nuclear matter relations

Employing a Skyrme interaction, the energy density of semi-infinite symmetric nuclear matter is written as

$$=\frac{\mathrm{}^2}{2m^{}}\tau +\frac{3}{8}t_0\rho ^2+\frac{1}{16}t_3\rho ^{\sigma +2}+\frac{1}{8}T_S|\rho |^2\frac{m^{}}{\mathrm{}^2}V_{\mathrm{so}}\rho |\rho |^2,$$
(12)
where $`V_{\mathrm{so}}=9W_0^2/16`$.

Denoting by $`\rho _0`$, $`ϵ_0`$, $`K_0`$ and $`m_0^{}`$ respectively the density, energy per particle, compressibility, and effective mass at saturation of the bulk symmetric nuclear matter, it turns out that the Skyrme parameters $`t_0`$, $`t_3`$, $`\sigma `$, and the combination $`T_0`$ can be written as:

$`T_0=\left({\displaystyle \frac{\mathrm{}^2}{m_0^{}}}{\displaystyle \frac{\mathrm{}^2}{m}}\right){\displaystyle \frac{1}{\rho _0}},`$ (13)
$`\sigma ={\displaystyle \frac{{\displaystyle \frac{1}{9}}K_0+ϵ_0+\left({\displaystyle \frac{\mathrm{}^2}{10m}}{\displaystyle \frac{2\mathrm{}^2}{15m_0^{}}}\right)k_F^2(0)}{ϵ_0+\left({\displaystyle \frac{3\mathrm{}^2}{10m}}{\displaystyle \frac{\mathrm{}^2}{5m_0^{}}}\right)k_F^2(0)}},`$ (14)
$`t_3={\displaystyle \frac{16}{\rho _0^{1+\sigma }}}{\displaystyle \frac{1}{\sigma }}\left[ϵ_0+\left({\displaystyle \frac{3\mathrm{}^2}{10m}}{\displaystyle \frac{\mathrm{}^2}{5m_0^{}}}\right)k_F^2(0)\right],`$ (15)
$`t_0={\displaystyle \frac{8}{3\rho _0}}\left[ϵ_0{\displaystyle \frac{3\mathrm{}^2}{10m_0^{}}}k_F^2(0){\displaystyle \frac{1}{16}}t_3\rho _0^{1+\sigma }\right],`$ (16)
where $`\rho _0=2/(3\pi )^2k_F^3(0)`$. Once the parameters $`t_0,t_3,\sigma `$ and the combination $`T_0`$ have been fixed, the surface energy only depends on the combination $`T_S`$ and the spin-orbit strength $`W_0`$. The surface energy can be written as an integral over the density

$`ϵ_S`$ $`=`$ $`8\pi r_0^2{\displaystyle _0^{\rho _0}}𝑑\rho \left[{\displaystyle \frac{\mathrm{}^2}{36m}}{\displaystyle \frac{5}{36}}T_0\rho +{\displaystyle \frac{1}{8}}T_S\rho {\displaystyle \frac{m^{}}{\mathrm{}^2}}V_{\mathrm{so}}\rho ^2\right]^{1/2}`$
$`\left[{\displaystyle \frac{3\mathrm{}^2}{10m^{}}}k_F^2+{\displaystyle \frac{3}{8}}t_0\rho +{\displaystyle \frac{1}{16}}t_3\rho ^{\sigma +1}\right]^{1/2},`$
where $`r_0=\left[3/(4\pi \rho _0)\right]^{1/3}`$ is the unit radius, and a Thomas-Fermi approximation up to $`\mathrm{}^2`$-order has been used to replace the kinetic energy density.

Then, the following relation involving the symmetry energy $`ϵ_I`$ can be used to relate $`x_0`$ to the already known parameter combinations and to the free parameters $`(x,y,z)`$:

$$t_0x_0=\left(\frac{3}{2}x+\frac{9}{10}y+\frac{2}{5}(T_02T_S)\right)\alpha _1\rho _0^{2/3}\frac{1}{6}z\rho _0^\sigma 2C_1(\rho _0)\frac{4ϵ_I}{\rho _0}.$$
(18)
To simplify the notation used in Section 2 we have introduced the following functions of the density:

$`C_0(\rho )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m}}+T_0\rho ,`$ (19)
$`C_1(\rho )`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{m}}+T_0\rho \right){\displaystyle \frac{\alpha _1}{\rho ^{1/3}}}+{\displaystyle \frac{1}{4}}t_0+{\displaystyle \frac{1}{24}}t_3\rho ^\sigma ,`$ (20)
$`C_2(\rho )`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{}^2}{m}}+4T_0\rho \right){\displaystyle \frac{\alpha _2}{\rho ^{1/3}}}{\displaystyle \frac{1}{2}}t_0{\displaystyle \frac{1}{24}}(\sigma +1)(\sigma +2)t_3\rho ^\sigma `$ (21)
$`C_3(\rho )`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{}^2}{m}}+T_0\rho \right){\displaystyle \frac{\alpha _2}{\rho ^{1/3}}}+{\displaystyle \frac{1}{2}}t_0+{\displaystyle \frac{1}{12}}t_3\rho ^\sigma ,`$ (22)
which contain only known parameters. We have also introduced the constants

$$\alpha _1=\left(\frac{\pi ^4}{12}\right)^{1/3},\alpha _2=\left(\frac{\pi ^4}{3}\right)^{1/3}$$
(23)
## Appendix B Landau parameters

The Landau parameters in symmetric nuclear matter are:

$`F_0`$ $`=`$ $`\left({\displaystyle \frac{3}{4}}t_0+{\displaystyle \frac{1}{16}}(\sigma +1)(\sigma +2)t_3\rho ^\sigma \right){\displaystyle \frac{2m^{}k_F}{\mathrm{}^2\pi ^2}}F_1,`$ (24)
$`G_0`$ $`=`$ $`\left({\displaystyle \frac{1}{4}}t_0(2x_01)+{\displaystyle \frac{1}{24}}t_3\rho ^\sigma (2x_31)\right){\displaystyle \frac{2m^{}k_F}{\mathrm{}^2\pi ^2}}G_1,`$ (25)
$`F_0^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{4}}t_0(2x_0+1){\displaystyle \frac{1}{24}}t_3\rho ^\sigma (2x_3+1)\right){\displaystyle \frac{2m^{}k_F}{\mathrm{}^2\pi ^2}}F_1^{},`$ (26)
$`G_0^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{4}}t_0{\displaystyle \frac{1}{24}}t_3\rho ^\sigma \right){\displaystyle \frac{2m^{}k_F}{\mathrm{}^2\pi ^2}}G_1^{},`$ (27)
$`F_1`$ $`=`$ $`3T_0{\displaystyle \frac{m^{}}{\mathrm{}^2}}\rho ,`$ (28)
$`G_1`$ $`=`$ $`3T_1{\displaystyle \frac{m^{}}{\mathrm{}^2}}\rho ,`$ (29)
$`F_1^{}`$ $`=`$ $`3T_2{\displaystyle \frac{m^{}}{\mathrm{}^2}}\rho ,`$ (30)
$`G_1^{}`$ $`=`$ $`3T_3{\displaystyle \frac{m^{}}{\mathrm{}^2}}\rho ,`$ (31)
where the $`T_i`$’s are the following parameter combinations:

$`T_0`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[3t_1+5t_2+4y\right],`$ (32)
$`T_1`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[2x+2yt_1+t_2\right],`$ (33)
$`T_2`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[2x2y+t_1t_2\right],`$ (34)
$`T_3`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[t_1t_2\right],`$ (35)
and the effective mass is:

$$\frac{\mathrm{}^2}{m^{}}=\frac{\mathrm{}^2}{m}+T_0\rho .$$
(36)
The Landau parameters in neutron matter are:

$`F_0^{(n)}`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}t_0(1x_0)+{\displaystyle \frac{1}{24}}(\sigma +1)(\sigma +2)t_3\rho ^\sigma (1x_3)\right){\displaystyle \frac{m_n^{}k_F}{\mathrm{}^2\pi ^2}}F_1^{(n)},`$ (37)
$`G_0^{(n)}`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}t_0(x_01)+{\displaystyle \frac{1}{12}}t_3\rho ^\sigma (x_31)\right){\displaystyle \frac{m_n^{}k_F}{\mathrm{}^2\pi ^2}}G_1^{(n)},`$ (38)
$`F_1^{(n)}`$ $`=`$ $`3(T_0T_2){\displaystyle \frac{m_n^{}}{\mathrm{}^2}}\rho ,`$ (39)
$`G_1^{(n)}`$ $`=`$ $`3(T_1T_3){\displaystyle \frac{m_n^{}}{\mathrm{}^2}}\rho ,`$ (40)
and the effective mass is:

$$\frac{\mathrm{}^2}{m_n^{}}=\frac{\mathrm{}^2}{m}+(T_0T_2)\rho .$$
(41)
*

*

*

*