text:
# Reduction of the sign problem using the meron-cluster approach.

## I Introduction

To stochastically study a quantum problem using a quantum Monte Carlo method it is necessary to transform it to a form that is similar to a classical statistical problem. The sign problem appears when this transformation leads to a weight function that is not positive definite. As quantum Monte Carlo methods have become increasingly efficientEvertz et al. (1993); Kawashima et al. (1994); Prokof’ev et al. (1996); Beard and Wiese (1996); Sandvik (1999); Syljuåsen and Sandvik (2002) there is a notable lack of progress in solving the sign problem. The sign problem severely limits the number of models that can be studied using quantum Monte Carlo methods, and in particular there are only very few models of interacting fermions in higher dimensions that are accessible to existing algorithms.Blankenbecler et al. (1981)

The recent development of the so-called meron-cluster solutionChandrasekharan and Wiese (1999) has extended the range of models where the sign problem can be avoided. This method uses the properties of loop quantum Monte Carlo algorithms to establish a one-to-one mapping between configurations with negative weight and corresponding configurations with positive weight. These contributions cancel each other and a fraction of the phase space with a positive definite weight function is left, which can be sampled with no sign problem.

Unfortunately the meron solution works for only a rather limited class of models.Chandrasekharan et al.  The main purpose of the present paper is to show that the meron concept can be applied also to models where the sign problem is not eliminated. We demonstrate that in a wider class of models it is possible to cancel out part of the negative configurations, and thereby reduce the sign problem. We investigate a model in an intermediate regime between a point in parameter space where the sign problem can be solved completely and a point where the meron-cluster algorithm cannot be applied. Our results show that the meron-cluster algorithm does indeed reduce the sign problem in this intermediate regime. The main focus of the study is frustrated spin models, but we also apply this method to spinless fermions.

The outline of the paper is as follows. In Sec. II the Monte Carlo algorithm is briefly explained. The sign problem and the meron-cluster algorithm are introduced in Sec. III. In Sec. IV a modified version of the stochastic series expansion is described. The origin of the sign problem for frustrated spin systems and fermions is discussed in Sec. V. In Sec. VI we demonstrate how the meron solution affects the average sign for a range of models where the sign problem cannot be eliminated. We conclude with summary and discussions in Sec. VII.

## II The quantum Monte Carlo algorithm

In order to explain the meron-cluster solution introduced in Sec. IV we here give a summary of the stochastic series expansion (SSE) methodSandvik and Kurkijärvi (1991); Sandvik (1999); Syljuåsen and Sandvik (2002).

Consider a lattice model described by a Hamiltonian $`H`$. In the SSE method the partition function $`Z`$ is Taylor expanded,

$$Z=\underset{\alpha }{}\underset{m=0}{\overset{\mathrm{}}{}}\frac{\beta ^m}{m!}\alpha |(H)^m|\alpha ,$$
(1)
where $`|\alpha `$ are states in which the above matrix elements can be calculated and $`\beta `$ denotes the inverse temperature.

For sake of clarity we will now consider a one-dimensional ferromagnetic Heisenberg model,

$$H=\underset{i=1}{\overset{N}{}}S_i^zS_{i+1}^z+\frac{1}{2}\left(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+\right),$$
(2)
where $`N`$ denotes th number of sites. The Hamiltonian is rewritten as a sum over diagonal and off-diagonal operators,

$$H=\underset{i=1}{\overset{N}{}}(H_{1,i}+H_{2,i}),$$
(3)
where

$$H_{1,i}=S_i^zS_{i+1}^z+C$$
(4)
and

$$H_{2,i}=\frac{1}{2}(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+),$$
(5)
where $`C`$ is a constant inserted to assure that the expectation value $`\alpha |H_{1,i}|\alpha `$ is positive for all states $`|\alpha `$.

To simplify the Monte Carlo update we introduce an additional unit operator $`H_{0,0}=1`$. Inserting the Hamiltonian given by Eq. (3) into Eq. (1), and truncating the sum at $`m=L`$ we obtain

$$Z=\underset{\alpha }{}\underset{S_L}{}\frac{\beta ^n(Ln)!}{L!}\alpha |\underset{k=1}{\overset{L}{}}H_{a_k,i_k}|\alpha ,$$
(6)
where $`n`$ stands for the number of non-unit operators, and $`S_L`$ denotes a sequence of operator-indices

$$S_L=(a_1,i_1),(a_2,i_2),\mathrm{},(a_L,i_L),$$
(7)
with $`a_k=1,2`$ and $`i_k=1,\mathrm{},N`$, or $`(a_k,i_k)=(0,0)`$.

The Monte Carlo procedure must sample the space of all states $`|\alpha `$, and all sequences $`S_L`$. The simulation starts with some random state $`|\alpha `$ and an operator string containing only unit operators. One Monte Carlo step consists of a diagonal and an off-diagonal update. In the diagonal update attempts are made to exchange unit and diagonal operators sequentially at each position in the operator string. The probability for inserting or deleting a diagonal operator in the operator string is given by detailed balance.Sandvik (1999)

The off-diagonal update, also called loop update, is carried out with $`n`$ fixed. Each bond operator $`H_{i_k}=H_{1,i_k}+H_{2,i_k}`$ acts only on two spins, $`S_{i_k}`$ and $`S_{i_k+1}`$. We can therefore rewrite the matrix elements in Eq. (6) as a product of $`n`$ terms, called vertices,

$$M(\alpha ,S_L)=\underset{k=1}{\overset{n}{}}W_k,$$
(8)
where the vertex weight $`W_p`$ is defined as

$`W_p=S_{i_p}^z(p)S_{i_p+1}^z(p)|H_{i_p}|S_{i_p}^z(p1),S_{i_p+1}^z(k1)`$
(9)
where $`S_i^z(p)`$ denotes the state of spin $`i`$ in a propagated state, defined by

$$|\alpha (p)\underset{k=1}{\overset{p}{}}H_{a_k,i_k}|\alpha .$$
(10)
A vertex thus consists of four spins, called the legs of the vertex, and an operator. Each term in the expansion in Eq. (6) can be viewed as a sequence of vertices. An example of one term for a four-site chain is shown in the left part of Fig. 1.

The principles of the off-diagonal update are: one of the $`n`$ vertices is chosen at random and one of its four legs is randomly selected as the entrance leg. The spin of the entrance leg is flipped. One of the legs of the operator is chosen as the exit leg, and its state is also changed. The exit leg is chosen with a probability calculated from the weight of the obtained vertex.Sandvik (1999) Thereafter the vertex list is sequentially searched for the next vertex that includes the exit spin. This spin becomes the entrance leg of the next vertex and the procedure is continued until the original entrance leg is reached. During one Monte Carlo step the loop update is repeated until on average half of the vertices have been updated.

In the method described above the spin states are altered as the loop is constructed. During one Monte Carlo step a given spin can be part of several different loops, or it may be of none. In a few special cases, such as for the isotropic Heisenberg model, the propagation of each loop through the lattice is deterministic, meaning that there is only one possible exit leg for each entrance leg.Sandvik (1999) In these special cases it is possible to divide the whole space-time lattice up into loops so that each and every spin belongs to only one loop. An example of such a configuration is shown in the right part of Fig. 1. The loop update can then be modified to identifying the unique loop structure and flipping each loop with probability one half.

## III The Sign Problem

In this section we show how the sign problem appears in quantum Monte Carlo simulations and introduce the recent meron approach to solving the sign problem. We start by considering a general form of an expectation value that can be calculated by Monte Carlo methods,

$$A=\frac{_iA(x_i)W(x_i)}{_iW(x_i)}=A(x)_W,$$
(11)
where the weight function $`W(x_i)`$ and $`A(x_i)`$ depend on the configuration $`x_i`$. When the coordinates $`x_i`$ are sampled according to relative weight, the expectation value is given by the average value of $`A(x)`$ as indicated in the last part of Eq. (11). The sign problem appears if the weight function is not positive. In this case the sampling can be done using the absolute value of the weight,

$$A=\frac{As_{|W|}}{s_{|W|}},$$
(12)
where $`s`$ denotes the sign of the weight function and equals $`\pm 1`$. However, in many cases of physical interest the average sign approaches zero exponentially as the system size is increased. The above expectation value will then suffer from very large statistical fluctuations since it becomes a ratio of two small numbers. Let us now consider how negative weight functions appear for quantum mechanical systems. The weight function $`W(\alpha ,S_L)`$ corresponding to the partition function given by Eq. (6) is

$$W(\alpha ,S_L)=\frac{\beta ^n(Ln)!}{L!}\alpha |\underset{k=1}{\overset{L}{}}H_{a_k,i_k}|\alpha ,$$
(13)
This is strictly positive, and for the ferromagnet there is no sign problem. Let us next consider an antiferromagnet. In this case the diagonal and off-diagonal operators are of the form

$$H_{1,b}=S_{i(b)}^zS_{j(b)}^z+C,$$
(14)
and

$$H_{2,b}=\frac{1}{2}(S_{i(b)}^+S_{j(b)}^{}+S_{i(b)}^{}S_{j(b)}^+).$$
(15)
By adjusting the constant $`C`$ it is still possible to have $`\alpha |H_{1,b}|\alpha 0`$. However, for the off-diagonal operator the expectation value $`\alpha |H_{2,b}|\alpha 0`$, and the sign must be taken into account. If there is an odd number of off-diagonal operators in the configuration the sign of the weight function will be negative. Due to the periodic boundary conditions in the imaginary time direction the number of off-diagonal operators on a square lattice is always even and there is no sign problem. However, if the system is frustrated, as on a triangular lattice, the sign problem appears for the antiferromagnetic spin model. For a fermionic system the anticommutator rules must be taken into account and the sign of the configuration changes sign every time two fermion world lines wrap around each other in imaginary time. We therefore see that both for frustrated spin models and fermionic models the sign problem enters into the loop update. Flipping a loop can cause the number of spin flipping operators to change parity, or it can cause two fermions to permute and thereby change the sign. Loops that cause the sign to change are called merons,Chandrasekharan and Wiese (1999) and in some cases one can, in effect, solve the sign problem by avoiding configurations that include merons.

In order to explain the meron solution we need to revisit the loop update. As was pointed out at the end of the previous section it is sometimes possible to divide the lattice up into a unique loop structure. Expectation values can then be calculated using so-called improved estimators, which are averages over all possible loop configurations. If the number of loops in the system is given by $`N_L`$ there are $`2^{N_L}`$ configurations that can be reached by flipping the loops, since all the loops can be in two states. The expectation value in Eq. (12) can therefore be rewritten as

$$A=\frac{A(x)s(x)_{|W|}}{s(x)_{|W|}},$$
(16)
where the double expectation brackets denote an average over the different loop configurations

$$s(x)A(x)_{|W|}=\frac{1}{2^{N_L}}\underset{l=1}{\overset{2^{N_L}}{}}s(x_l)A(x_l)_{|W|}.$$
(17)
If certain criteria are fulfilled, the expectation value can be expressed as

$$A=A(x)\delta _{n_M,0}_{|W|},$$
(18)
where $`n_M`$ is the number of merons. Therefore only configurations without sign changing loops give non-zero contributions to the expectation value, and the sign problem is, in effect, solved.

Let us examine the necessary conditions for this to be the case:

1. The lattice can be divided up into loops so that each spin belongs to one and only one loop.
2. The weight must not change when the loops are flipped.
3. The loops must affect the sign independently.
4. The zero-meron sector must be positive definite.
5. The expectation value of the operator is unchanged when a loop is flipped.

Together these conditions place severe restrictions on which models can be studied with the meron solution. Our aim is to examine if the conditions can be relaxed to allow for a more general algorithm. Of the five conditions the last one is the least severe. Many operators for which this condition does not hold can be expressed by introducing the two-meron sector.Chandrasekharan and Wiese (1999); Henelius and Sandvik (2000) Examples of operators for which the last condition holds are the energy and heat capacity, which in SSE are given by

$`E`$ $`=`$ $`{\displaystyle \frac{n}{\beta }}`$ (19)
$`C`$ $`=`$ $`n^2n^2n.`$ (20)
where $`n`$ is the number of non-unit operators in the operator string.

The first three conditions are severe and we are not aware of a way to relax these restrictions. In this study we examine the case when the fourth condition is not met, and the zero-meron sector is not positive definite. If this is not the case Eq. (18) must be replaced by

$$A=\frac{A(x)s(x)\delta _{n_M,0}_{|W|}}{s(x)\delta _{n_M,0}_{|W|}},$$
(21)
where $`s(x)`$ is the sign of the configuration in the zero meron sector. The higher meron sectors do not contribute to the average and the sign problem is therefore reduced, but not eliminated. To find out how the average sign changes by leaving out the non-zero meron sectors is the aim of this investigation. It is not clear whether the exponential character of the sign problem will change, or how the system size will affect the average sign. This depends on the relative weight of the different meron sectors and on the average sign in the zero-meron sector.

In order to investigate systems where condition four does not hold it is necessary to satisfy the first three criteria. As was stated at the end of Sec. II there are a few special cases where a unique loop structure exists, and in these cases the first condition is automatically satisfied. For more general models, where one can choose between different exit legs, this is not true, since there no longer exists a unique way to divide the lattice into loops. It is, however, possible to cover the lattice with loops so that each spin belongs to one and only one loop. This can be done by always choosing an exit leg that is not already part of any loop. Since a unique loop structure does not exist one also has to sample all possible loop configurations, and in order to implement the meron solution efficiently this has to be done without leaving the zero and two-meron sectors. We have implemented a loop update that inherently divides the lattice up in separate loops. This can be done as long as the bounce process, where the entrance and exit legs coincide, can be neglected. With the advent of directed loopsSyljuåsen and Sandvik (2002) there now exists a way to eliminate the bounce process in many models of interest. In the next section we will introduce the idea of directed loops and intoroduce models where the first three, but not the forth, conditions are satisfied.

## IV SSE Loop Construction

Following Ref. Syljuåsen and Sandvik, 2002 we here derive the directed loops equations that enable us to neglect the bounce process, and explicitly show that the weight in the extended space of directed loops is unaltered when a loop is flipped. Then we will describe a modification to the SSE loop update and demonstrate how the meron approach can be efficiently implemented also for models where there is not a unique way of dividing the lattice into loops.

By considering that each possible loop has a “time-reversed” counterpart it was shownSyljuåsen and Sandvik (2002) that a loop move satisfies detailed balance if, for every vertex that the loop passes through,

$$W_sP(s,ij)=W_s^{}P(s^{},ji),$$
(22)
where the vertex weight $`W_s`$ for a given spin configuration $`s`$ is defined by Eq. (9) and $`P(s,ij)`$ is the probability to choose exit leg $`j`$ given entrance leg $`i`$ for a given spin configuration $`s`$.

The main idea of the method of directed loops is to attach weight also to the link that connects the entrance and exit spins at a traversed vertex. The weight of a vertex in this extended space can then be written as $`W(s,i,j)`$. We are allowed to introduce such auxiliary variables if the sum over auxiliary variables reduces to the original weight

$$\underset{j}{}W(s,i,j)=W_s.$$
(23)
We furthermore require that the weight in this extended space is not affected by flipping a loop, which translates to

$$W(s,i,j)=W(s^{},j,i).$$
(24)
Considering the criteria for detailed balance, Eq. (22) we can now relate the weights in the original and extended space as

$$W(s,e,x)=W_sP(s,e,x).$$
(25)
We have therefore shown that the total weight in our extended configuration space is unchanged when flipping a loop, and the second condition for using the meron solution is therefore fulfilled. The method of directed loops is a way to assign probabilities to the different possible exit legs given an entrance leg. The probabilities are determined by solving the system of equations that Eq. (22) and Eq. (24) generate. We want to emphasis that for many models it is possible to solve the equations so that the bounce process is forbiddenSyljuåsen and Sandvik (2002), and as we now show this makes it possible to divide the lattice up into separate loops.

Next we describe a new procedure that enables us to make this division. In our new update the physical structure of the loops is determined and updated during the diagonal move, while each loop is flipped with probability one half in the off-diagonal update. The diagonal move is quite different from the one described in section II and we here describe it in some detail.

The operator string is, as before, updated sequentially, and if a unit operator is encountered an attempt to insert a diagonal operator is made. If an operator is inserted its four legs are linked to the existing loops. If the bounce is eliminated there are three possible ways to link the entrance and exit legs of a vertex, see Fig. 2. The decision of which configuration to use is based on detailed balance equations, which states that the probability of different vertices is proportional to their weight given by Eq. 9. It is enough to treat an arbitrary leg as “entrance leg” to the vertex and determine an exit leg. With two legs connected together in this manner the other two legs are also automatically linked. If instead a diagonal operator is encountered, an attempt is made to exchange it for a unit operator, as in the standard diagonal update. If the removal is not accepted an attempt is made to alter the way in which the legs are connected. In the same way as when inserting a new diagonal operator an “exit leg” is determined for an arbitrary “entrance leg”, a procedure which may change the existing links . This change of loop structure is also done if an off-diagonal operator is encountered. We graphically demonstrate this modified update in Fig. 3. On the left side of Fig. 3 a term in the SSE expansion is divided up into loops. A possible outcome after a diagonal update is shown on the right. In the next section we will apply this update to models that suffer from the sign problem.

## V Models

### V.1 Frustrated spin

The main focus of this article is frustrated spin systems. The Hamiltonian, $`H`$, is given by,

$$H=\underset{i,j}{}\mathrm{\Delta }S_i^zS_j^z+\frac{1}{2}(S_i^+S_j^{}+S_j^+S_i^{}),$$
(26)
where the sum runs over all nearest neighbors. For $`\mathrm{\Delta }\{1,1\}`$ the directed loops equations can be solved so that the bounce is eliminated, and using the method described in the previous section the first two criteria for using the meron solution are fulfilled. The frustration is introduced by having spins positioned on a triangular lattice. The off-diagonal coupling is antiferromagnetic and as described previously the number of off-diagonal operators may be odd on a frustrated lattice and there is therefore a sign problem. In Fig. 4 an example of such a configuration with an odd number of off-diagonal operators is shown. The third criteria states that the loops need to be independent in their effect on the sign. If flipping a loop causes the total number of off-diagonal operators to change parity the sign changes and the loop is a meron. The only way this can happen is if the loop traverses an odd number of vertices. The number of transversed vertices is independent of any other loops and it is therefore clear that the loops are independent in their effect on the sign, which establishes the third criteria. We have therefore shown that the first three criteria for the application of the meron solution are fulfilled.

The fourth criteria requires that the zero-meron sector has a positive definite weight. This is not, in general, true for this model. The configuration shown in Fig. 4 has negative weight, yet the only loop in the system is not a meron. The only point where the fourth criteria hold true is for $`\mathrm{\Delta }=1`$, a point in parameter space which has been extensively studied.Henelius and Sandvik (2000) At this point it is possible to exclude all but one update, update C in Fig. 2. As we move away from $`\mathrm{\Delta }=1`$ the negative part of the zero-meron sector grows. At the other extreme, $`\mathrm{\Delta }=1`$, it is not possible to eliminate the bounces without introducing a non-ergodicity. If the bounce is eliminated the only allowed update is update A in Fig 2. With only this update there are no loops which pass through an odd number of operators, and thus there are no merons. Still the zero-meron sector, the only sector left, contains both positive and negative parts, but there is no way to switch between them and therefore the whole SSE space is not sampled. For all other values $`\mathrm{\Delta }\{1,1\}`$ the meron solution can be applied, and this model constitutes an ideal testing ground for a further study of the meron solution since, by adjusting a single parameter, we can move from a point where the meron solution eliminates the sign problem ($`\mathrm{\Delta }=1`$) to a point where the loop update becomes non-ergodic ($`\mathrm{\Delta }=+1`$). However, since much effort has been made to solve the sign problem in fermion models, we will next describe the application to a system of spinless fermions.

### V.2 Spinless fermions

Besides the frustrated spin systems we have also studied spinless fermions on a square lattice. In this context a meron is a loop which permutes an even number of fermions when it is flipped. The loops must affect the sign independently of each other for the meron method to work, as stated by the third constraint described in Sec. III. For the frustrated spin models described above this was always the case, but it is not so for fermionic models. As discussed previously there are three different ways to traverse a vertex, see Fig. 2. For fermions update C make the loops dependent of each other.Chandrasekharan et al.  To solve this problem one can give the completely empty vertex a negative weight, but in this case update B causes the loops to be dependent. We are not aware of a way solve this dependency problem and one is thus restricted to models where one of the two updates B and C can be forbidden.

In the work by Chandrasekharan et. alChandrasekharan et al. , two models are given where only one type of vertex update is allowed and where the zero meron sector, therefore, is positive. We have studied a model where it is possible to exclude update B and the bounce but where both update A and C must be allowed. The Hamiltonian for this model is

$`H`$ $`=`$ $`{\displaystyle \underset{i,j}{}}c_i^+c_j^{}+c_j^+c_i^{}+{\displaystyle \frac{4}{3}}(n_i1/2)(n_j1/2)`$ (27)
$`{\displaystyle \frac{1}{3}}(n_i+n_j),`$
where $`n_i`$ is the occupation on site $`i`$ and $`c_i^{}`$ ($`c_i^+`$) are the ordinary annihilation (creation) fermion operators. By adjusting the constant, $`C`$ added to the diagonal operator, the empty vertex is given a negative weight and the loops are independent of each other. Therefore we have a fermionic example of a model where the first three criteria for using the meron solution are fulfilled, but where the zero-meron sector is not positive definite. In the next section we analyze how the meron solution affects the average sign in these cases where the zero-meron sector is not positive definite.

## VI Results

First we consider the spin model described by Eq. (26). We have calculated the expectation value of the sign for frustrated spin models with different values of the constant $`\mathrm{\Delta }`$. The expectation value for a three-site system is shown as a function of temperature in Fig. 5. The simulation is performed in the zero and two meron sectors and expectation values are shown both including and excluding the two-meron sector. As can be seen in the figure the sign decreases exponentially in both cases. Asymptotically it appears that the sign in the zero-meron sector is increased by a constant factor as compared to the case of including the two-meron sector. The average sign is greatest close to the point $`\mathrm{\Delta }=1`$, where the sign problem is eliminated, and it decreases as the point $`\mathrm{\Delta }=1`$ is approached. This is to be expected since the meron-cluster solution cannot be applied in the present form at the Heisenberg point($`\mathrm{\Delta }=1`$).

We have also calculated how the expectation value of the sign changes with the system size. The results are shown in Fig. 6. In this calculation $`\mathrm{\Delta }=0.9`$ is used. The average sign appears to decrease approximately exponentially also with system size.

Besides the two already mentioned expectation values of the sign we have also calculated the expectation value of the sign without the use of merons. This value is slightly smaller than the one for the two-meron sector. The difference is small since the sector with two merons has a very much larger weight than the sector with more than two merons, at least for the temperatures and system sizes that we have studied. In Fig. 7 we show the relative weight for the zero-meron, the two-meron and the sector with more than two merons, calculated for a system with ten spins at different temperatures. The weight of the sector with more than two merons increases with a decreased temperature as a consequence of a larger configuration which results in a larger number of loops. In the figure we also indicate the relation between the number of positive and negative values in the different sectors. For the sectors with merons there is, by definition, an equal amount of positive and negative values. As can be expected from Fig. 5 the ratio of positive to negative weight in the zero-meron sector approaches one as the temperature decreases. We also note that the relative weight of the zero-meron sector decreases quite rapidly with lower temperatures. There is no sector with an odd number of merons. If there where such a sector configurations would exist where a flip of all the loops would result in a sign change. A flip of all the loops corresponds to a change of all spin states, an operation which does not change the number of off-diagonal operators.

We have also done a similar calculation for the spinless fermion model given by Eq. (27). The result for the expectation value of the sign is shown in Fig. 8 as a function of temperature. Two different system sizes are studied, two times two and four times four sites. Also here the average sign appears to decrease exponentially with inverse temperature and system size both in the zero- and two-meron sectors. Asymptotically it seems that the average sign in the zero-meron sector again is increased by a constant factor when leaving out the higher meron sectors.

We have studied the relative weight of the sectors with different number of merons also for the fermionic system. For the fermions there are configurations with only one meron. This is due to the fact that the empty vertex is given a negative weight and flipping all the loops in a configuration may change the parity of the number of empty vertices. A comparison with exact diagonalization indicates that one only needs to include either the one- or the two-meron sector in addition to the zero-meron sector. In figure 9 the relative weight for the zero-, one-, two-, and higher merons sectors are presented. At high temperatures the relative weight of the zero-meron sector again dominates, while at lower temperatures the weight in the higher meron sectors increases.

## VII Summary and discussion

We have shown that as long as it is possible to divide the system into independent loops, the meron-cluster approach can be used to decrease the sign problem even when the zero-meron sector is not positive definite. We have applied this method to both frustrated spin systems and spinless fermions. An intermediate regime between a point in parameter space where the meron-cluster algorithm eliminates the sign problem and a point where it cannot be applied is studied. In this intermediate regime the exponential character of the sign problem persists, but one can increase the average sign by a constant factor by limiting measurements to the zero-meron sector. The method is probably of most practical use in the vicinity of points where the sign problem can be eliminated using the meron solution. In a large scale application the weight of the two-meron sector should be decreased by a reweighting techniqueHenelius and Sandvik (2000) to obtain better statistics. To be able to use this algorithm we have combined stochastic series expansion with the concept of directed loops.

###### Acknowledgements.

We are grateful to A. Sandvik for stimulating discussions. The work was supported by the Swedish Research Council and the Göran Gustafsson foundation.


format: text

text:
# “I’m sorry Dave, I’m afraid I can’t do that”: Linguistics, Statistics, and Natural Language Processing circa 2001To appear in the National Research Council study on the Fundamentals of Computer Science. This is an April 2003 version.

## Ambiguity and language analysis

> At last, a computer that understands you like your mother.
> – 1985 McDonnell-Douglas ad

The snippet quoted above indicates the early confidence at least one company had in the feasibility of getting computers to understand human language. But in fact, that very sentence is illustrative of the host of difficulties that arise in trying to analyze human utterances, and so, ironically, it is quite unlikely that the system being promoted would have been up to the task. A moment’s reflection reveals that the sentence admits at least three different interpretations:

1. The computer understands you as well as your mother understands you.
2. The computer understands that you like your mother.
3. The computer understands you as well as it understands your mother.

That is, the sentence is ambiguous; and yet we humans seem to instantaneously rule out all the alternatives except the first (and presumably the intended) one. We do so based on a great deal of background knowledge, including understanding what advertisements typically try to convince us of. How are we to get such information into a computer?

A number of other types of ambiguity are also lurking here. For example, consider the speech recognition problem: how can we distinguish between this utterance, when spoken, and “… a computer that understands your lie cured mother”? We also have a word sense ambiguity problem: how do we know that here “mother” means “a female parent”, rather than the Oxford English Dictionary-approved alternative of “a cask or vat used in vinegar-making”? Again, it is our broad knowledge about the world and the context of the remark that allows us humans to make these decisions easily.

Now, one might be tempted to think that all these ambiguities arise because our example sentence is highly unusual (although the ad writers probably did not set out to craft a strange sentence). Or, one might argue that these ambiguities are somehow artificial because the alternative interpretations are so unrealistic that an NLP system could easily filter them out. But ambiguities crop up in many situations. For example, in “Copy the local patient files to disk” (which seems like a perfectly plausible command to issue to a computer), is it the patients or the files that are local?222Or, perhaps, the files themselves are patient? But our knowledge about the world rules this possibility out. Again, we need to know the specifics of the situation in order to decide. And in multilingual settings, extra ambiguities may arise. Here is a sequence of seven Japanese characters:

Since Japanese doesn’t have spaces between words, one is faced with the initial task of deciding what the component words are. In particular, this character sequence corresponds to at least two possible word sequences, “president, both, business, general-manager” (= “a president as well as a general manager of business”) and “president, subsidiary-business, Tsutomu (a name), general-manager” (= ?). It requires a fair bit of linguistic information to choose the correct alternative.333To take an analogous example in English, consider the non-word-delimited sequence of letters “theyouthevent”. This corresponds to the word sequences “the youth event”, “they out he vent”, and “the you the vent”.

To sum up, we see that the NLP task is highly daunting, for to resolve the many ambiguities that arise in trying to analyze even a single sentence requires deep knowledge not just about language but also about the world. And so when HAL says, “I’m afraid I can’t do that”, NLP researchers are tempted to respond, “I’m afraid you might be right”.

## Firth things first

But before we assume that the only viable approach to NLP is a massive knowledge engineering project, let us go back to the early approaches to the problem. In the 1940s and 1950s, one prominent trend in linguistics was explicitly empirical and in particular distributional, as exemplified by the work of Zellig Harris (who started the first linguistics program in the USA). The idea was that correlations (co-occurrences) found in language data are important sources of information, or, as the influential linguist J. R. Firth declared in 1957, “You shall know a word by the company it keeps”.

Such notions accord quite happily with ideas put forth by Claude Shannon in his landmark 1948 paper establishing the field of information theory; speaking from an engineering perspective, he identified the probability of a message’s being chosen from among several alternatives, rather than the message’s actual content, as its critical characteristic. Influenced by this work, Warren Weaver in 1949 proposed treating the problem of translating between languages as an application of cryptography (see \[other chapter\]), with one language viewed as an encrypted form of another. And, Alan Turing’s work on cracking German codes during World War II led to the development of the Good-Turing formula, an important tool for computing certain statistical properties of language.

In yet a third area, 1941 saw the statisticians Frederick Mosteller and Frederick Williams address the question of whether it was Alexander Hamilton or James Madison who wrote some of the pseudonymous Federalist Papers. Unlike previous attempts, which were based on historical data and arguments, Mosteller and Williams used the patterns of word occurrences in the texts as evidence. This work led up to the famed Mosteller and Wallace statistical study which many consider to have settled the authorship of the disputed papers.

Thus, we see arising independently from a variety of fields the idea that language can be viewed from a data-driven, empirical perspective — and a data-driven perspective leads naturally to a computational perspective.

## A “C” change

However, data-driven approaches fell out of favor in the late 1950’s. One of the commonly cited factors is a 1957 argument by linguist (and student of Harris) Noam Chomsky, who believed that language behavior should be analyzed at a much deeper level than its surface statistics. He claimed,

> It is fair to assume that neither sentence (1) \[Colorless green ideas sleep furiously\] nor (2) \[Furiously sleep ideas green colorless\] … has ever occurred …. Hence, in any \[computed\] statistical model … these sentences will be ruled out on identical grounds as equally “remote” from English.444Interestingly, this claim has become so famous as to be self-negating, as simple web searches on “Colorless green ideas sleep furiously” and its reversal will show. Yet (1), though nonsensical, is grammatical, while (2) is not.

That is, we humans know that sentence (1), which at least obeys (some) rules of grammar, is indeed more probable than (2), which is just word salad; but (the claim goes), since both sentences are so rare, they will have identical statistics — i.e., a frequency of zero — in any sample of English. Chomsky’s criticism is essentially that data-driven approaches will always suffer from a lack of data, and hence are doomed to failure.

This observation turned out to be remarkably prescient: even now, when billions of words of text are available on-line, perfectly reasonable phrases are not present. Thus, the so-called sparse data problem continues to be a serious challenge for statistical NLP even today. And so, the effect of Chomsky’s claim, together with some negative results for machine learning and a general lack of computing power at the time, was to cause researchers to turn away from empirical approaches and toward knowledge-based approaches where human experts encoded relevant information in computer-usable form.

This change in perspective led to several new lines of fundamental, interdisciplinary research. For example, Chomsky’s work viewing language as a formal, mathematically-describable object has had lasting impact on both linguistics and computer science; indeed, the Chomsky hierarchy, a sequence of increasingly more powerful classes of grammars, is a staple of the undergraduate computer science curriculum. Conversely, the highly influential work of, among others, Kazimierz Adjukiewicz, Joachim Lambek, David K. Lewis, and Richard Montague adopted the lambda calculus, a fundamental concept in the study of programming languages, to model the semantics of natural languages.

## The empiricists strike back

By the ’80s, the tide had begun to shift once again, in part because of the work done by the speech recognition group at IBM. These researchers, influenced by ideas from information theory, explored the power of probabilistic models of language combined with access to much more sophisticated algorithmic and data resources than had previously been available. In the realm of speech recognition, their ideas form the core of the design of modern systems; and given the recent successes of such software — large-vocabulary continuous-speech recognition programs are now available on the market — it behooves us to examine how these systems work.

Given some acoustic signal, which we denote by the variable $`a`$, we can think of the speech recognition problem as that of transcription: determining what sentence is most likely to have produced $`a`$. Probabilities arise because of the ever-present problem of ambiguity: as mentioned above, several word sequences, such as “your lie cured mother” versus “you like your mother”, can give rise to similar spoken output. Therefore, modern speech recognition systems incorporate information both about the acoustic signal and the language behind the signal. More specifically, they rephrase the problem as determining which sentence $`s`$ maximizes the product $`P(a|s)\times P(s)`$. The first term measures how likely the acoustic signal would be if $`s`$ were actually the sentence being uttered (again, we use probabilities because humans don’t pronounce words the same way all the time). The second term measures the probability of the sentence $`s`$ itself; for example, as Chomsky noted, “colorless green ideas sleep furiously” is intuitively more likely to be uttered than the reversal of the phrase. It is in computing this second term, $`P(s)`$, where statistical NLP techniques come into play, since accurate estimation of these sentence probabilities requires developing probabilistic models of language. These models are acquired by processing tens of millions of words or more. This is by no means a simple procedure; even linguistically naive models require the use of sophisticated computational and statistical techniques because of the sparse data problem foreseen by Chomsky. But using probabilistic models, large datasets, and powerful learning algorithms (both for $`P(s)`$ and $`P(a|s)`$) has led to our achieving the milestone of commercial-grade speech recognition products capable of handling continuous speech ranging over a large vocabulary.

But let us return to our story. Buoyed by the successes in speech recognition in the ’70s and ’80s (substantial performance gains over knowledge-based systems were posted), researchers began applying data-driven approaches to many problems in natural language processing, in a turn-around so extreme that it has been deemed a “revolution”. Indeed, now empirical methods are used at all levels of language analysis. This is not just due to increased resources: a succession of breakthroughs in machine learning algorithms has allowed us to leverage existing resources much more effectively. At the same time, evidence from psychology shows that human learning may be more statistically-based than previously thought; for instance, work by Jenny Saffran, Richard Aslin, and Elissa Newport reveals that 8-month-old infants can learn to divide continuous speech into word segments based simply on the statistics of sounds following one another. Hence, it seems that the “revolution” is here to stay.

Of course, we must not go overboard and mistakenly conclude that the successes of statistical NLP render linguistics irrelevant (rash statements to this effect have been made in the past, e.g., the notorious remark, “Every time I fire a linguist, my performance goes up”). The information and insight that linguists, psychologists, and others have gathered about language is invaluable in creating high-performance broad-domain language understanding systems; for instance, in the speech recognition setting described above, a better understanding of language structure can lead to better language models. Moreover, truly interdisciplinary research has furthered our understanding of the human language faculty. One important example of this is the development of the head-driven phrase structure grammar (HPSG) formalism — this is a way of analyzing natural language utterances that truly marries deep linguistic information with computer science mechanisms, such as unification and recursive data-types, for representing and propagating this information throughout the utterance’s structure. In sum, computational techniques and data-driven methods are now an integral part both of building systems capable of handling language in a domain-independent, flexible, and graceful way, and of improving our understanding of language itself.

### Acknowledgments

Thanks to the members of the CSTB Fundamentals of Computer Science study — and especially Alan Biermann — for their helpful feedback. Also, thanks to Alex Acero, Takako Aikawa, Mike Bailey, Regina Barzilay, Eric Brill, Chris Brockett, Claire Cardie, Joshua Goodman, Ed Hovy, Rebecca Hwa, John Lafferty, Bob Moore, Greg Morrisett, Fernando Pereira, Hisami Suzuki, and many others for stimulating discussions and very useful comments. Rie Kubota Ando provided the Japanese example. The use of the term “revolution” to describe the re-ascendance of statistical methods comes from Julia Hirschberg’s 1998 invited address to the American Association for Artificial Intelligence. I learned of the McDonnell-Douglas ad and some of its analyses from a class run by Stuart Shieber. All errors are mine alone. This paper is based upon work supported in part by the National Science Foundation under ITR/IM grant IIS-0081334 and a Sloan Research Fellowship. Any opinions, findings, and conclusions or recommendations expressed above are those of the authors and do not necessarily reflect the views of the National Science Foundation or the Sloan Foundation.

## For further reading

Charniak, Eugene. 1993. Statistical Language Learning. MIT Press.

Jurafsky, Daniel and James H. Martin. 2000. Speech and Language Processing: An Introduction to Natural
Language Processing, Computational Linguistics, and Speech Recognition. Prentice Hall. Contribut-
ing writers: Andrew Kehler, Keith Vander Linden, and Nigel Ward.

Manning, Christopher D. and Hinrich Schütze. 1999. Foundations of Statistical Natural Language Process-
ing. The MIT Press.


text:
# Neutrino condensates at center of galaxies as background for the MSW mechanism

## I Introduction

The nature of matter condensates at galaxy centers  is nowadays widely discussed since no current proposed model seems fully able to fit the data. For example, recent observations toward the Center of our Galaxy indicate that it hosts a super-massive dark compact object, Sagittarius A (Sgr A) , which is an extremely loud radio source of estimated mass $`M2.61\times 10^6M_{}`$ and radius $`R0.016`$pc$`30.6`$lds  and several evidence for the presence of similar objects have been found in other galaxies, in quasars, and active galactic nuclei in a mass-range $`10^6÷10^9M_{}`$ .

The most natural candidates to explain such peculiar objects could be either a single super-massive black hole or a very compact cluster of stellar size black holes . This last case has been ruled out by stability criteria  which give maximal lifetimes of the order of $`10^8`$ years which are much shorter than the estimated age of standard galaxies . The first hypothesis is much more supported since similar super-massive black holes allow to explain the central dynamics of several galaxies as M87 , or NGC4258 . However, if Sgr A were a super-massive black hole, its luminosity should be more than $`10^{40}`$erg/s, in contrast to observations which give a bolometric luminosity of $`10^{37}`$erg/s (this is the so called blackness problem or the black hole on starvation). Besides, the most recent observations probe the gravitational potential at a radius larger than $`4\times 10^4`$ Schwarzschild radii of a black hole of mass $`2.6\times 10^6M_{}`$ , so that unambiguous proof that the super-massive object at the center of Galaxies is a black hole is still lacking.

Several alternative models have recently appeared in literature: for example, in , the hypothesis that the Galactic Center could consist of a super-massive boson star is investigated and some proposals have been done in order to seek a signature like, for example, the Cerenkov effect .

Viollier et al.  proposed that the dark matter at the center of galaxies could be made by non-baryonic matter (massive neutrinos) which interacts gravitationally forming super-massive balls in which the degeneracy pressure of fermions balances their self–gravity. Such neutrino balls could be formed in the early epochs during a first–order gravitational phase transition and their dynamics could be reconciled with some adjustments to the Standard Model of Cosmology . The neutrino mass required for fitting observational data by mean of super-massive degenerate neutrino stars is 10keV$`m_{\nu _\tau }25`$keV, $`m_{\nu _\tau }`$ being the mass of $`\tau `$ neutrinos making up the stars. The gravitational effects of such a structure on the stars orbiting around it can be investigated by gravitational lensing , while astrophysical constraints and signatures can be investigated by next generation of $`X`$-ray satellites as XEUS or Constellation-X . In , a more physical model is discussed taking into account a neutrino ball in gravitational equilibrium of a semi-degenerate fermion gas. Density and pressure within the ball are defined by adopting a formalism based on a distribution function in phase space, which allows to consider neutrinos with a degeneracy degree varying from the center to the border of the system. Limiting cases are fully degenerate fermion systems and the classical isothermal spheres well-known in literature.

In this paper, we take into account the fact that a condensate of massive neutrino could naturally act as the background for the Mikheyev-Smirnov-Wolfenstein (MSW) mechanism which gives rise to neutrino oscillations. In Sect.II, we outline the simplest neutrino star model in order to define the environment where such oscillations should happen. The MSW effect and the active-sterile oscillations are discussed in Sect.III. Conclusions are drawn in Sect.IV.

## II The Neutrino Stars Model

In this section, we shortly recall the simplest model of heavy neutrino condensates, bound by gravity (for details, see ). In the Thomas–Fermi model for fermions, the Fermi energy $`E_F`$ and the gravitational potential binding the system are related by (in natural units)

$$\frac{k_F^2(r)}{2m_\nu }m_\nu \mathrm{\Phi }(r)=E_F=m_\nu \mathrm{\Phi }(r_0),$$
(1)
where $`\mathrm{\Phi }(r)`$ is the gravitational potential, $`k_F`$ is the Fermi wave number and $`\mathrm{\Phi }(r_0)`$ is a constant chosen to cancel the gravitational potential for vanishing neutrino density. The constant $`r_0`$ is the estimated size of the halo of the ball. If one takes into account a degenerate Fermi gas, one gets $`k_F(r)=\left(6\pi ^2n_\nu (r)/g_\nu \right)^{1/3}`$, where $`n_\nu (r)`$ is the neutrino number density, assumed being the same for neutrinos and antineutrinos within the halo. The number $`g_\nu `$ is the spin degeneracy factor. If in the center of the neutrino condensate there is a baryonic star (which is approximated as a point source), the gravitational potential will obey a Poisson equation where neutrinos (and antineutrinos) are the source terms,

$$\mathrm{}\mathrm{\Phi }=4\pi Gm_\nu n_B,$$
(2)
where $`n_B=N/V`$ is the neutrino background number density with $`N`$, the number of neutrinos making up the condensate, $`N=M/m_\nu `$, and $`V=4\pi R^3/3`$ its volume. We have

$$n_B=\frac{N}{V}=\frac{3Mm_\nu }{4\pi R^3}.$$
(3)
Eq.(2) is valid everywhere except at the origin. In the case of spherical symmetry, the Poisson equation reduces to a radial Lané–Emden differential equation .

The general solution of (2), has scaling properties and it is able to reproduce the observations . In fact, a degenerate neutrino star of mass $`M=2.6\times 10^6M_{}`$, consisting of neutrinos with mass $`m12.0`$keV for $`g_\nu =4`$, or $`m14.3`$keV for $`g_\nu =2`$, does not contradict the observations. Considering a standard accretion disk, the data are in agreement with the model if Sgr A is a neutrino star with radius $`R=30.3`$ ld ($`10^5`$ Schwarzschild radii) and mass $`M=2.6\times 10^6M_{}`$ with a luminosity $`L10^{37}`$erg/s. Similar results hold also for the dark object ($`M3\times 10^9M_{}`$) inside the center of M87.

For our purposes, we assume a simple model where the neutrino density is constant around a point-like baryonic mass (e.g. a baryonic star). Such a model has been also considered in  where lensing experiments have been proposed for probing the existence of neutrino balls, without invoking the presence of super-massive black holes at the center of our Galaxy. If $`n_B=constant`$, it follows, from Eq. (2), that $`\mathrm{\Phi }r^2`$, in agreement with the Newton theorem applied to a spherical mass distribution. Such a situation is achieved considering a Fermi gas at temperature $`T=0`$ .

## III MSW Effect and $`\nu _\tau \nu _s`$ Oscillation

In the framework of Violler et al. model, it is natural to investigate the possibility that neutrinos crossing the neutrino background making up the supermassive condensate, could oscillate owing to the mechanism equivalent to MSW effect . To be more specific, assuming that $`\tau `$ neutrinos are the constituent of neutrino balls, we explore the possibility that a beam of $`\tau `$ neutrinos ($`\nu _\tau `$) propagating in the $`\nu _\tau `$-background, could oscillate in sterile neutrinos $`\nu _s`$ ($`\nu _\tau \nu _s`$) through a neutral current interaction. The neutrino $`\nu _s`$ has the property that does not undergo electroweak interactions.

According to the standard model, the $`\nu \nu `$ scattering is described by the effective Lagrangian density

$$_{eff}=\frac{G_F}{\sqrt{2}}[\overline{\nu }_\tau \gamma ^\mu (1\gamma ^5)\nu _\tau ][\overline{\nu }_B\gamma ^\mu (1\gamma ^5)\nu _B]$$
(4)
where $`\nu _B`$ indicates background neutrino, i.e. $`\nu _B\nu _\tau `$. For simplicity in what follows we write $`\nu `$ instead of $`\nu _L(1\gamma ^5)\nu /2`$. Neutrinos forming the condensate are non relativistic so that the axial current contribution is negligible, as well as the spatial components of the 4-current since it is related to the average velocity of $`\nu _B`$. Then, by considering the average over the background, the 0-component of the 4-current gives, in agreement with standard procedure ,

$$<\overline{\nu }_B\gamma ^0\nu _B>=<\nu _B^{}\nu _B>n_B.$$
(5)
Eq. (5) adds an amounts of $`\sqrt{2}G_Fn_B`$ to the enegy of neutrinos propagating within the neutrino star. Equation of evolution of neutrinos is (for simplicity we consider only the radial motion)

$$i\frac{d}{dr}\left(\begin{array}{c}\nu _\tau \\ \nu _s\end{array}\right)=\mathrm{\Lambda }\left(\begin{array}{c}\nu _\tau \\ \nu _s\end{array}\right),$$
(6)
where the mixing matrix is defined as

$$\mathrm{\Lambda }\left(\begin{array}{cc}\frac{\mathrm{\Delta }m^2}{2E}\mathrm{cos}2\theta +a& \frac{\mathrm{\Delta }m^2}{2E}\mathrm{sin}2\theta \\ \frac{\mathrm{\Delta }m^2}{2E}\mathrm{sin}2\theta & \frac{\mathrm{\Delta }m^2}{2E}\mathrm{cos}2\theta \end{array}\right),$$
(7)
where $`a\sqrt{2}G_Fn_B`$. The mass eigenstates $`|\nu _1>`$ and $`|\nu _2>`$ are determined by diagonalizing the mixing matrix in (7): One writes the mass eigenstates as a superposition of flavor eigenstates

$`|\nu _1>`$ $`=`$ $`\mathrm{cos}\stackrel{~}{\theta }|\nu _\tau >\mathrm{sin}\stackrel{~}{\theta }|\nu _s>,`$ (8)
$`|\nu _2>`$ $`=`$ $`\mathrm{sin}\stackrel{~}{\theta }|\nu _\tau >+\mathrm{cos}\stackrel{~}{\theta }|\nu _s>,`$ (9)
where the mixing angle $`\stackrel{~}{\theta }`$ is defined as

$$\mathrm{tan}2\stackrel{~}{\theta }=\frac{\mathrm{\Delta }m^2\mathrm{sin}2\theta }{\mathrm{\Delta }m^2\mathrm{cos}2\theta 2aE}.$$
(10)
The corresponding eigenvalues are

$$\lambda _{1,2}=\pm \frac{\mathrm{\Delta }m^2}{4E}\sqrt{(\mathrm{cos}2\theta b)\mathrm{cos}2\theta +\mathrm{sin}^22\theta },$$
(11)
being $`b=2aE/\mathrm{\Delta }m^2`$. The resonance condition occurs when

$$\mathrm{cos}2\theta =b.$$
(12)
We remember that the neutrino density $`n_B`$ is assumed constant. From Eq. (12), and using Eq.(3) one infers the relation

$$\mathrm{\Delta }m^27\times 10^{14}\left(\frac{30.6\text{ld}}{R}\right)^3\frac{E}{\text{GeV}}\frac{M}{M_{}}\frac{\text{eV}}{m_{\nu _\tau }}\text{eV}^2.$$
(13)
Neutrino stars are characterized by the parameters $`R30.6`$ld and $`M2.6\times 10^6M_{}`$, thus Eq. (13) becomes

$$\mathrm{\Delta }m^21.4\times 10^7\frac{E}{\text{GeV}}\frac{\text{eV}}{m_{\nu _\tau }}\text{eV}^2.$$
(14)
The result does depend on neutrino energy and the mass of the background neutrinos ($`m_{\nu _\tau }10`$keV). From (14) we derive the allowed solutions for $`\mathrm{\Delta }m^2`$:

$$E\text{MeV},\mathrm{\Delta }m^210^{14}\text{eV}^2,$$

$$E\text{GeV},\mathrm{\Delta }m^210^{11}\text{eV}^2,$$

$$E\text{TeV},\mathrm{\Delta }m^210^8\text{eV}^2,$$

which agree with the best fit for active-sterile neutrinos of $`\tau `$ type, $`\mathrm{\Delta }m^210^6`$ obtained from the big bang nucleosynthesis predictions and in the absence of lepton asymmetry. The case of neutrinos with very high energy

$$E\text{PeV},\mathrm{\Delta }m^210^5\text{eV}^2,$$

is excluded.

In presence of lepton symmetry, the above bound $`\mathrm{\Delta }m^210^6`$ becomes $`\mathrm{\Delta }m^21`$, and the range of energy $`1÷10^6`$GeV is permissible in order that neutrino stars might induce an appreciable MSW effect

These results have non trivial consequences in relation to neutrino asymmetries in the early Universe, widely discussed in Refs. .

## IV Discussion and Conclusions

In this paper we have investigated the active-sterile neutrino oscillations for neutrinos propagating in a neutrino star which acts as background matter. Taking into account the neutral current weak interaction, we have studied the resonance condition, in analogy to MSW mechanism for solar neutrinos.

The basic results here obtained are two-fold. From one side, our results show that the present data on neutrino oscillations do not contradict a possible existence of such exotic super-massive objects. Actually much more investigation on this results is necessary. First of all, we should take into account more realistic models where the density of the condensate changes with the radius . Furthermore, one should investigate the $`\nu _e:\nu _\mu :\nu _\tau `$ ratio which becomes $`\nu _e:\nu _\mu :\nu _\tau ^{}`$, where $`\nu _\tau ^{}`$ are the $`\tau `$-neutrinos number reduced (with respect to the expected ones) as a consequence of the $`\nu _\tau \nu _s`$ oscillation induced by the neutrino star environment (for this topic, see also ). On the other side, our results show that oscillations affected by matter, i.e. induced by neutrino stars, could have effects on the generation of neutrino asymmetries in the Universe . The latter occur as a consequence of active-sterile neutrino oscillations, which generate a discrepancy of neutrino and antineutrino number density (such a subject is widely discussed in Refs. ). The lepton number of a neutrino flavor $`f`$ is defined as

$$L_f=\frac{n_{\nu _f}n_{\overline{\nu }_f}}{n_\gamma (T)},$$
(15)
where $`n_{\nu _f}`$ ($`n_{\overline{\nu }_f}`$) is the number density of neutrinos (antineutrinos), and $`n_\gamma `$ one of photons at temperature $`T`$. Constraints on $`L_f`$ are derived from the cosmic microwaves background and from big bang nucleosynthesis. In the case of neutrinos $`\nu _\tau `$, which are of interest in our case, bounds on the lepton number are $`|L_{\nu _\tau }|6.0`$. The creation of the lepton number determines a strong modification on the upper limit of the mass squared difference value. In fact, it goes from $`\mathrm{\Delta }m^210^6`$eV2, coming from the big bang nucleosynthesis predictions, to $`\mathrm{\Delta }m^21`$eV2 .

As final comment, it is worth to point out that if at the center of galaxies a black hole is present, then no oscillation phenomena could occur, unless a violation of the equivalence principle is invoked . This is opposite to results obtained in this paper, where $`\nu _\tau `$-background matter affects oscillations giving a new scenario and a strong signature for future investigations aimed to probe the (possible) existence of such exotic objects, the neutrino balls, at the center of galaxies.


text:
# Fracture of a viscous liquid

## Abstract

When a viscous liquid hits a pool of liquid of same nature, the impact region is hollowed by the shock. Its bottom becomes extremely sharp if increasing the impact velocity, and we report that the curvature at that place increases exponentially with the flow velocity, in agreement with a theory by Jeong and Moffatt. Such a law defines a characteristic velocity for the collapse of the tip, which explains both the cusp-like shape of this region, and the instability of the cusp if increasing (slightly) the impact velocity. Then, a film of the upper phase is entrained inside the pool. We characterize the critical velocity of entrainment of this phase and compare our results with recent predictions by Eggers.

We all observed that diving in a swimming-pool produces air bubbles, while entering slowly in the same pool does not: the motion of a solid penetrating a liquid at a high speed largely deforms the surface of this liquid, and induces air entrainment. This phenomenon often defines the maximum speed of coating of a solid, which enters a bath before leaving it coated. Generally air entrainment must be avoided because of the resulting bubbles (which are all the more harmful since the coating solution generally contains surfactants, and thus is likely to transform into a foam). Different experiments using either fibers or plates as solids indeed prove the existence of a threshold (in velocity) for air entrainment Burley (1992), but there is today no theoretical picture for understanding quantitatively this phenomenon.

Quite similarly, when a viscous liquid enters a pool of same nature (as it occurs for example for the stem wave of a moving ship), air entrainment may occur (producing bubbles, which ultimately burst – a serious cause of noise for a ship). Such a bubble production was observed for a liquid jet impinging a bath of same nature, as reported by Lin and Donelly, who studied experimentally the minimum velocity of air entrainment as a function of different characteristics of the jet (diameter, viscosity or surface tension) Lin and Donelly (1981). This question is important in many industrial processes where a viscous liquid (typically molten glass, metal or polymer) is poured inside a mould – and there again, the formation of bubbles is detrimental. A remarkable geometry for studying this phenomenon, and more generally the phenomenology of the regimes of deformation of the liquid bath, was proposed by Joseph Joseph et al. (1991); Joseph (1992). Two counter-rotative cylinders are partially immersed in a bath, whose upper surface is observed. At a low speed of rotation, the interface between the cylinders is hardly deformed, but increasing this speed leads to the formation of a cusp-like shape. This device thus produces a singularity at a free interface, and more generally the observation and description of similar singularities has been a very active field for the past years, either theoretically Jeong and Moffat (1992) or experimentally - let us quote for example selective withdrawal Cohen and Nagel (2002), the pinch-off of a jet Cohen et al. (1999), drops running-down on a tilted solid Podgorski et al. (2001) or jet eruptions Zeff et al. (2000).

Surface tension naturally opposes the formation of such singularities at a liquid/fluid interface, and for the problem of the cusp between two counter-rotative cylinders, Jeong and Moffatt calculated in a classical paper the (non-zero) radius $`r`$ of the tip, in the limit of small Reynolds numbers Jeong and Moffat (1992). They derived an analytical expression, which can be summarized using a physical argument proposed by Hinch Jeong and Moffat (1992). In this two-dimensional geometry, the tip can be seen as a thin (hemi-)cylinder, of radius $`r`$. The viscous drag (per unit length) on a cylinder scales as $`\eta V/\mathrm{ln}(d/r)`$, denoting $`d`$ as an external length, $`\eta `$ as the liquid viscosity and $`V`$ as the flow velocity around the cylinder. On the other hand, the surface tension $`\gamma `$ draws this cylinder (as it does on the edge of a free sheet of liquid) with a force $`2\gamma `$ (there again, per unit length). Balancing both these forces yields as a tip radius:

$$rd\mathrm{exp}(Ca)$$
(1)
where $`Ca`$ is the capillary number ($`Ca=\eta V/\gamma `$).

This law first helps to understand why a tip forms around a well-defined velocity: the cusp becomes very sharp when $`V`$ is of the order of $`\gamma /\eta `$, a logical quantity since both these parameters play antagonist roles in the formation of the tip. But it also predicts that $`r`$ essentially vanishes above this velocity, and Eggers pointed out that this should lead to the crack of the tip, because of the lubrication pressure of the upper fluid, forced to flow in this very confined region Eggers (2001). He showed that beyond a critical radius for the tip, the pressure in the cusp breaks its stationary shape. Thus, this critical radius $`r_c`$ must increase with the viscosity $`\eta _0`$ of the upper fluid, and Eggers found that $`r_c`$ should scale as $`(\eta _0/\eta )^{4/3}`$ Eggers (2001). Together with Moffatt’s result (Eq. 1), this yields as a threshold velocity $`V_c`$ of entrainment of the upper phase:

$$V_c\frac{\gamma }{\eta }\mathrm{ln}\left(\eta /\eta _0\right)$$
(2)
We could have expected (naively) that entrainment of the upper fluid would occur when the viscous forces overcome capillary ones (i.e. $`V_c\gamma /\eta `$). Eggers’ theory (Eq. 2) shows that taking into account the flow of the upper phase in the cusp should induce a logarithmic correction to this simple statement. As a consequence, the threshold velocity should depend on the viscosity of the upper phase, i.e. on its nature. This might sound a posteriori quite logical, but remarkably, this dependence is expected to be extremely weak.

Here, we investigate experimentally the different situations related to this problem. We first describe the macroscopic deformation of the interface, due to the incoming flow of liquid. Then, we focus on the smaller scale of the tip. Finally, we characterize the condition of entrainment of the upper phase.

We used a horizontal cylinder of radius $`R=2`$ cm, half-immersed in a bath of glycerol (or in some cases silicone oil), of density $`\rho `$ and viscosity $`\eta `$ around 1 Pa.s (Figure 1). The whole device was set in a transparent tank, and the upper fluid was either air or light oil (of density $`\rho _0<\rho `$ and viscosity $`\eta _0\eta `$). The interfacial tension between the two fluids, denoted as $`\gamma `$, was measured for each experiment. The cylinder could be rotated at a controlled velocity $`\mathrm{\Omega }`$, and the motion filmed with a digital camera.

A viscous liquid film of thickness $`h`$ (of typically 5 mm, thus thick enough to be independant of the nature of the roller surface) is dragged out of the bath. This layer impinges the bath at a velocity $`V\mathrm{\Omega }(R+h)`$, and deforms the free surface at this place by a depth $`L`$. We denote $`r`$ as the radius of curvature of the bottom of this hollow region. Figure 2 shows the different shapes obtained with a silicon oil ($`\eta =0.97`$ Pa.s), and increasing the speed of the cylinder.

It is observed that both the characteristic lengths $`L`$ and $`r`$ are dynamic in essence: the depth of immersion $`L`$ increases with the speed of cylinder, while the tip radius $`r`$ decreases. Moreover, these lengths are very different from each other: while the depth of immersion is around 1 cm, the tip has a radius of curvature much smaller than 1 mm. This radius can become very small when the speed is large, as first reported by Joseph et al. Joseph et al. (1991) in a similar situation. But we do not observe as they did a clear threshold velocity for the formation of a cusp – both the lengths $`L`$ and $`r`$ varying continuously as a function of the cylinder velocity.

We first focus on the macroscopic length $`L`$, namely the depth of immersion of the tip into the bath, and report in Fig. 3 the variation of $`L`$ as a function of the impacting velocity $`V`$, for a bath of silicon oil ($`\eta =0.97`$ Pa.s) .

The depth $`L`$ of immersion of the tip is found to increase with $`V`$. For speeds larger than 10 cm/s, $`L`$ is several millimetres, which is significantly larger than the capillary length $`\kappa ^1=\sqrt{\gamma /\mathrm{\Delta }\rho g}1.5`$ mm (denoting $`\mathrm{\Delta }\rho `$ as the difference of density between both fluids). We conclude that capillary forces can be neglected in this regime of large velocities. The viscous friction which pulls downward the tip can be written dimensionally $`\eta V/L^2`$, while gravitational forces per unit volume (which pull it upwards) are $`\mathrm{\Delta }\rho g`$. Balancing both these forces yields:

$$L\left(\eta V/\mathrm{\Delta }\rho g\right)^{1/2}$$
(3)
Equation 3 defines some kind of viscous length (on the model of $`\kappa ^1`$ defined just above, replacing capillary effects by viscous ones). It describes fairly well the data in Figure 3, where $`L`$ is observed to increase as the square root of $`V`$. Some deviations appear at small $`V`$ : the deformation then is small ($`L\kappa ^1`$), and capillary effects must be considered. Since they oppose the formation of the hollow region, they indeed make $`L`$ smaller than predicted by Eq. 3.

We then focused on the small scale of the deformation, namely the tip itself, and took macrophotographs of this region. A numerical camera was placed behind an optical microscope lighted by a classical white lamp, and the whole device first calibrated using micrometric fibers. The tip was observed to be slightly deeper close to the container walls (which allowed us to make conveniently our measurements), but remains far from it (about 15 mm because of the lubrification of the wall. As reported in Figure 4, the tip has a typical radius smaller than 100 µm, which decreases strongly with the flow velocity.

Analyzing the contour of the interface close to the tip allows a determination of the curvature: we can fit this contour by a parabolic function $`y=ax^2`$, from which we deduce the curvature $`r=1/2a`$. We display in Figure 5 the measured tip radii as a function of the capillary number $`Ca=\eta V/\gamma `$ for air and oil tips formed in a bath of glycerol.

The tip radius is found to decrease in a continuous way, as a function of the flow velocity $`V`$. More precisely, its logarithm varies linearly with $`V`$ – in agreement with Moffatt’s law (Eq. 1), on nearly two orders of magnitude for the tip radius. This exponential behavior helps understanding why the cusp-like shape apparently sets for a well-defined velocity, of the order of $`\gamma /\eta `$, although the actual variation of the radius is continuous with $`V`$. But it also suggests that slightly above this velocity, the cusp should become extremely sharp. This is not the case: a small increase of the roller velocity leads to the destruction of the tip, and the upper phase then is entrained, as reported in Figure 6. Note that the entrainment occurs simultaneously all along the tip, showing the 2D nature of the phenomenon.

We denote $`V_c`$ as the threshold velocity and $`Ca_c`$ as the corresponding capillary number. Our experiments consisted in determining $`Ca_c`$. Firstly, we varied the bath viscosity $`\eta `$, using mixtures of water and glycerol and air as the upper fluid. $`\eta `$ was measured for each solution. The corresponding data are displayed in black in Figure 7. Secondly, we used pure glycerol for the bath, and various light oils (such as hexane, hexadecane and silicone oils) as the upper phase. We measured $`\gamma `$ for each oil (it was found to be quite constant, around $`26\pm 1`$ mN/m). The corresponding critical capillary numbers are represented by white points, in Figure 7, where all our results are plotted as a function of the viscosity ratio $`\eta _0/\eta `$, in a semi-logarithmic scale.

We clearly observe an effect of the viscosity $`\eta _0`$ of the upper fluid (open points): the larger $`\eta _0`$, the smaller the critical capillary number. Nevertheless, this dependency is weak, since it appears to be well described by a logarithmic law. It is also limited: if the upper fluid is too viscous (typically $`\eta _0/\eta >0.1`$), the lower liquid is not entrained on the roller. For a given upper fluid, varying the bath viscosity $`\eta `$ also leads to a (weak) variation of the threshold in capillary number – which thus is not simply given by a value of order unity. All the data lie on the same semi-logarithmic plot, in agreement with Eggers’ model (Eq. 2), for a variation on five orders of magnitude of the viscosity ratio $`\eta _0/\eta `$. However, the slope deduced from the fit is $`0.22`$, larger than $`1`$ the value predicted by Eq. 2.

On the whole, these experiments may help clarifying the question of the stationary pattern which sets as a viscous liquid hits itself, in a two-dimensional geometry. If the impinging velocity is large enough (i.e. for capillary numbers typically larger than 0.1), two very different dynamic lengths can be defined: a global one, which is the depth of liquid hollowed by the impact, and which results from a balance between viscous effects and gravity; and a local one, at the bottom of this hollow region which can be expressed as the radius of curvature at this place. As predicted by Moffatt Jeong and Moffat (1992), this radius decreases exponentially with the capillary number. This law defines a characteristic velocity for which the hollow region becomes very sharp, i.e. takes a cusp-like shape, as first reported by Joseph Joseph et al. (1991). At the same time, it predicts a collapse of the tip radius for finite velocities, which implies the crack of the tip (entrainment of the upper phase), as postulated by Moffatt Jeong and Moffat (1992) and shown theoretically by Eggers Eggers (2001). We indeed report here such a crack, and show that the critical velocity for which it occurs is mainly governed by the bath viscosity, but also depends (in a much weaker manner) on the viscosity of the upper phase – in agreement with Eggers’ theory.

Different studies could naturally complement this one. Firstly, it would be worth understanding (both experimentally and theoretically) the thickness of the film of air, in the regime of entrainment. The selection of this thickness could remind the Landau-Levich problem (selection of the thickness of the film adhering to a solid withdrawn out of a bath), but the problem here is purely dynamical, which raises an interesting (and open) question (as in selective withdrawal, for the radius of the emitted filament). The problem of the stability of this film, which eventually breaks in bubbles, also remains to be studied. Another issue concerns the practically important question of water, for which the threshold velocity is much higher, which would require to study the influence of inertia in all the different processes described all along this study. Secondly, we observed a deviation towards the law displayed in Figure 7 using a very viscous silicone oil as a bath (instead of glycerol), and air above. Then, for $`\eta >10^5\eta _0`$, we measured a critical capillary number much higher, typically of the order of 5 to 10 (and increasing rapidly as decreasing the ratio $`\eta _0/\eta `$). A possible cause of this discrepancy could be the polymeric nature of the liquid (for which the chains are all the longer since the viscosity is high), which could lead to non-Newtonian effects close to the tip - and thus to changes in the threshold of air entrainment. Using newtonian silicone oils of smaller viscosity (with either air or ethanol above) indeed yields data in very close agreement with the ones displayed in Figure 7. Finally, it would be useful to connect this study with similar ones in other dimensionality. In particular, Cohen stressed that for selective withdrawal (i.e. one-dimensional tip forming when sucking locally a liquid above an interface), the viscosity of the fluid forming the tip does not affect the value of the cut-off curvature Cohen (2002). This could be due to the fact that the Stokes law on a (hemi-)sphere does not contain any logarithmic dependency, as it does for a (hemi-)cylinder, which finally emphasizes that the effects reported here are two-dimensional in essence.

We thank Jens Eggers for stimulating discussions, Ghislaine Ducouret who helped us in the viscosity measurements and Gilles Jandeau and Daniel Détry for their technical help.


text:
# The Fe XXII I(11.92 Å)/I(11.77 Å) Density Diagnostic Applied to the Chandra High Energy Transmission Grating Spectrum of EX Hydrae

## 1 Introduction

Magnetic cataclysmic variables (mCVs; polars and intermediate polars) are bright X-ray sources because the highly supersonic \[free-fall velocity $`v_{\mathrm{ff}}=(2GM_{\mathrm{wd}}/R_{\mathrm{wd}})^{1/2}3600\mathrm{km}\mathrm{s}^1`$\] material raining down on the white dwarf primary passes through a strong shock where most of its kinetic energy is converted into thermal energy and the plasma is heated to a temperature $`T_\mathrm{s}=3GM_{\mathrm{wd}}\mu m_\mathrm{H}/8kR_{\mathrm{wd}}200`$ MK. Because of the magnetic funneling of the mass lost by the secondary, the factor-of-four density jump across the accretion shock, and the settling nature of the post-shock flow (wherein the pressure is roughly constant and the density scales inversely with the temperature; Frank, King, & Raine 1992), the hot plasma in mCVs is expected to be dense. For a mass-accretion rate $`\dot{M}=10^{15}\mathrm{g}\mathrm{s}^1`$ ($`L_\mathrm{X}=GM_{\mathrm{wd}}\dot{M}/2R_{\mathrm{wd}}3\times 10^{31}\mathrm{erg}\mathrm{s}^1`$) and a relative spot size $`f=0.1`$, the density of the flow immediately behind the shock is $`n=\dot{M}/4\pi fR_{\mathrm{wd}}^2\mu m_\mathrm{H}(v_{\mathrm{ff}}/4)10^{13}\mathrm{cm}^3`$.

The standard density diagnostic of high-temperature plasmas is the intensity ratio $`Rf/i`$ of the forbidden to intercombination lines of He-like ions (Gabriel & Jordan, 1969; Blumenthal, Drake, & Tucker, 1972; Porquet et al., 2001), but this diagnostic is compromised in mCVs for two reasons. First, the critical density of this ratio increases with $`Z`$ and hence temperature—the opposite of the trend in the accretion column—so the $`R`$ line ratio is effective over only a very narrow slice of the accretion column. Second, in mCVs and other ultraviolet-bright stars, photoexcitation competes with collisional excitation to depopulate the upper level of the forbidden line, so the $`R`$ line ratio can appear to be in the “high-density limit” if the radiation field is sufficiently strong at the appropriate wavelengths: for a plasma illuminated by a 30 kK blackbody, the $`R`$ line ratios of all elements through Mg lie in the high-density limit, regardless of the density (Mauche, 2002).

To circumvent the density/photoexcitation ambiguity of the $`R`$ line ratio of He-like ions, we undertook an investigation of potential density diagnostics of lines of Fe L-shell ions that are observable in the high quality Chandra High Energy Transmission Grating (HETG) spectrum of the intermediate polar EX Hya (Mauche, 2000, 2002). Toward this end, Mauche, Liedahl, & Fournier (2001) investigated the density, temperature, and photoexcitation dependence of the $`I(17.10\mathrm{\AA })/I(17.05\mathrm{\AA })`$ line ratio of Fe XVII, and found that the anomalous ratio observed in EX Hya can be explained if the electron density of its $`T_\mathrm{e}4`$ MK plasma is $`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}2\times 10^{14}\mathrm{cm}^3`$ or if the photoexcitation temperature $`T_{\mathrm{bb}}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}55`$ kK. In this Letter, we investigate the density, temperature, and photoexcitation dependence of the $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio of Fe XXII, which also is anomalous in EX Hya. We find that this line ratio is very insensitive to temperature and photoexcitation and that the observed ratio can be explained if the electron density of EX Hya’s $`T_\mathrm{e}12`$ MK plasma is $`n_\mathrm{e}1\times 10^{14}\mathrm{cm}^3`$.

## 2 Model Spectra

The density sensitivity of the X-ray lines of Fe XXII have been discussed previously by Doschek, Meekins, & Cowan (1973), Mason & Storey (1980), Doschek, Feldman, & Cowan (1981), Phillips et al. (1982), Fawcett et al. (1987), Phillips et al. (1996), and Wargelin et al. (1998) in studies of the X-ray spectra of solar and laboratory high-temperature plasmas. These diagnostics have not received more attention because of the relative weakness of Fe XXII emission lines (the Fe XXII ionization fraction peaks at just 22% and is greater than 1% over the relatively narrow temperature range $`T_\mathrm{e}7`$–26 MK; Mazzotta et al. 1998), the relative weakness of $`n2`$ lines for $`n4`$, and the fact that the critical densities of these transitions are of order $`10^{14}\mathrm{cm}^3`$, larger than the densities typically observed in the Sun, tokamaks, and electron beam ion traps. We became interested in the density sensitivity of the $`32`$ lines of Fe XXII because the critical density of this transition is comparable to the lower limit of the electron density derived from the Fe XVII $`I(17.10\mathrm{\AA })/I(17.05\mathrm{\AA })`$ line ratio in the Chandra HETG spectrum of EX Hya, and because the Fe XXII $`32`$ $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio is anomalous (the 11.92 Å line is unusually strong) in the Chandra HETG spectrum of EX Hya compared to that of the bright late-type binaries HR 1099, Capella, and $`\sigma ^2`$ CrB (Ayres et al. 2001, Phillips et al. 2001, and Osten et al. 2003, respectively).

Before wading into the details, we refer the reader to Figure 1, which illustrates the origin of the density dependence of the Fe XXII spectrum. At low densities ($`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }<\text{ }}{\text{ }\text{ }}10^{12}\mathrm{cm}^3`$), the Fe XXII electron population is primarily in the $`2s^22p_{1/2}`$ ground state, and collisional excitations are predominantly from the ground state into the $`2s2p^2`$ manifold (8 levels total) and the $`2s^23d_{3/2}`$ level, both of which decay primarily to ground, producing lines in the extreme ultraviolet (EUV) ($`\lambda \lambda =101`$–217 Å) and at 11.77 Å, respectively. However, these levels also decay to the $`2s^22p_{3/2}`$ first-excited level with approximately 15% probability, producing lines in the EUV ($`\lambda \lambda =114`$–349 Å) and at 11.93 Å; when that level decays to ground, it produces a line in the ultraviolet (UV) ($`\lambda =846`$ Å). As the density increases ($`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}10^{13}\mathrm{cm}^3`$), electron population builds up in the $`2s^22p_{3/2}`$ first-excited level because the M1 transition to ground is slow ($`A_{21}=1.47\times 10^4\mathrm{s}^1`$). At high densities, the $`2s^22p_{3/2}`$ first-excited level is fed primarily by radiative decays from the $`2s2p^2`$ manifold. Collisional excitations out of the first-excited level are primarily into the $`2s^23d_{5/2}`$ level, which decays primarily back to the first-excited level, producing a line at 11.92 Å. Consequently, the Fe XXII 11.92 Å line is relatively strong in the X-ray spectra of high-density plasmas.

To calculate quantitative models of the X-ray spectrum of Fe XXII, we employed the Livermore X-ray Spectral Synthesizer (LXSS), a suite of IDL codes that calculates spectral models of highly charged ions based on Hebrew University/Lawrence Livermore Atomic Code (HULLAC) atomic data. HULLAC calculates atomic wavefunctions, level energies, radiative transition rates $`A_{\mathrm{ji}}`$, and oscillator strengths $`f_{\mathrm{ij}}`$ according to the fully relativistic, multiconfiguration, parametric potential method (Klapisch, 1971; Klapisch et al., 1977). Electron impact excitation rate coefficients are computed quasi-relativistically in the distorted wave approximation (Bar-Shalom, Klapisch, & Oreg, 1988) assuming a Maxwellian velocity distribution. Our Fe XXII model includes electron impact excitation rate coefficients and radiative transition rates for E1, E2, M1, and M2 decays for levels with principal quantum number $`n5`$ and azimuthal quantum number $`l4`$ for a total of 228 levels. Using these data, LXSS calculates the level populations for a given temperature and density assuming collisional-radiative equilibrium; the line intensities are then simply the product of the level populations and the radiative decay rates.

In a preliminary report, we (Mauche, Liedahl, & Fournier, 2003) used LXSS to calculate the X-ray spectrum of Fe XXII for electron densities $`n_\mathrm{e}=10^{11}`$ to $`10^{17}\mathrm{cm}^3`$ and electron temperatures $`T_\mathrm{e}=6.3`$, 12.8, and 25.5 MK ($`\frac{1}{2}`$, 1, and 2 times the temperature at which the Fe XXII ionization fraction peaks; Mazzotta et al. 1998). At the resolution of the Chandra HETG, the $`2s^23d_{3/2}`$$`2s^22p_{1/2}`$ line at 11.77 Å is blended with the (very much weaker) $`2s^23d_{5/2}`$$`2s^22p_{1/2}`$ line at 11.76 Å, and the $`2s^23d_{5/2}`$$`2s^22p_{3/2}`$ line at 11.92 Å is blended with the (much weaker) $`2s^23d_{3/2}`$$`2s^22p_{3/2}`$ line at 11.93 Å, so we report as $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ the ratio of the sum of these blends. Figure 3 of Mauche, Liedahl, & Fournier (2003) shows that this line ratio has a critical density $`n_\mathrm{c}10^{14}\mathrm{cm}^3`$, is approximately 0.3 at low densities and 1.4 at high densities, and has only a weak temperature dependence.

### 2.1 Modifications to LXSS: Collisional Excitation

Motivated by the detailed study of Fournier et al. (2001) of the density-sensitive Fe lines in the EUV spectrum of the high-density Frascati Tokamak Upgrade plasma, for this Letter we investigated the effect on the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio of two modifications to the collisional excitation data used in LXSS. First, we investigated the effect of replacing, for all transitions between and among the $`2s^22p`$ and $`2s2p^2`$ levels of Fe XXII, the electron impact excitation rate coefficients computed with HULLAC with those of Zhang & Pradhan (1997), computed with the relativistic R-matrix method. At $`T_\mathrm{e}=12.8`$ MK, the Zhang & Pradhan electron impact excitation rate coefficients are larger than those of HULLAC by 1.95 for the $`2s^22p_{1/2}`$$`2s^22p_{3/2}`$ transition, 0.92–1.70 (average 1.27) for the eight $`2s^22p_{1/2}`$$`2s2p^2`$ transitions, 1.05–3.23 (average 1.55) for the eight $`2s^22p_{3/2}`$$`2s2p^2`$ transitions, and 0.45–8.99 (average 2.20) for transitions between levels within the $`2s2p^2`$ manifold (but only 1.00–1.02 for the four $`2s^22p`$$`2s^23d`$ transitions). Despite the significant differences between these excitation rate coefficients, we find that for $`T_\mathrm{e}=12.8`$ MK the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio increases by 7% at $`n_\mathrm{e}=10^{14}\mathrm{cm}^3`$, but the low- and high-density limits are unaffected. Second, we investigated the effect of adding proton excitations for transitions among the levels of the $`2s^22p`$ and $`2s2p^2`$ configurations. This was accomplished by fitting the proton impact excitation rate coefficients of Foster, Keenan, & Reid (1997) to an expression of the form $`Q(T_\mathrm{e})=10^{12}T_\mathrm{e}^{1/2}_{i=1}^3A_i\mathrm{exp}(B_i/T_\mathrm{e})`$, where $`Q(T_\mathrm{e})`$ is the proton impact excitation rate coefficient in units of $`\mathrm{cm}^3\mathrm{s}^1`$, and we have assumed that the proton and electron temperatures are equal. Adding proton excitations to the LXSS population kinetics calculation under the assumption that the proton number density $`n_\mathrm{p}=0.85n_\mathrm{e}`$ (as is appropriate for a high-temperature plasma with solar abundances), we find that for $`T_\mathrm{e}=12.8`$ MK the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio increases by 19% at $`n_\mathrm{e}=10^{14}\mathrm{cm}^3`$, and the high-density limit increases by 7% to approximately 1.5.

With both of these modifications to LXSS, we recalculated the X-ray spectrum of Fe XXII for electron densities $`n_\mathrm{e}=10^{11}`$ to $`10^{17}\mathrm{cm}^3`$ and electron temperatures $`T_\mathrm{e}=6.3`$, 12.8, and 25.5 MK. The upper panel of Figure 2 shows that the model Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio has a critical density $`n_\mathrm{c}5\times 10^{13}\mathrm{cm}^3`$, is approximately 0.3 at low densities and 1.5 at high densities, and has only a very weak temperature dependence.

### 2.2 Modifications to LXSS: Photoexcitation

To investigate the effects of photoexcitation on the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio, we added to the LXSS population kinetics calculation the photoexcitation rates $`(\pi e^2/m_\mathrm{e}c)f_{\mathrm{ij}}F_\nu (T)`$, where $`F_\nu (T)`$ is the spectrum of the photoexcitation radiation field. For simplicity, we assume that $`F_\nu (T)=(4\pi /h\nu )B_\nu (T_{\mathrm{bb}})`$ (i.e., the radiation field is that of a blackbody of temperature $`T_{\mathrm{bb}}`$) and the dilution factor of the radiation field is equal to $`\frac{1}{2}`$ (i.e., the X-ray emitting plasma is in close proximity to the source of the photoexcitation continuum). With these assumptions, we used LXSS to calculate the X-ray spectrum of Fe XXII for an electron temperature $`T_\mathrm{e}=12.8`$ MK, electron densities $`n_\mathrm{e}=10^{11}`$ to $`10^{17}\mathrm{cm}^3`$, and photoexcitation temperatures $`T_{\mathrm{bb}}=20`$, 40, 60, 80, and 100 kK. The lower panel of Figure 2 shows that below $`T_{\mathrm{bb}}=60`$ kK the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio has essentially no photoexcitation dependence, but as the photoexcitation temperature increases the line ratio increases at low electron densities and the critical density shifts to higher densities, while the high-density limit is unaffected. Compared to the Fe XVII $`I(17.10\mathrm{\AA })/I(17.05\mathrm{\AA })`$ line ratio, the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio is insensitive to photoexcitation because the $`2s^22p_{1/2}`$$`2s^22p_{3/2}`$ transition in the UV is not optically allowed ($`f_{12}=3.15\times 10^6`$) and the $`2s^22p_{1/2}`$$`2s2p^2`$ transitions lie at such short wavelengths in the EUV.

## 3 Observations

The Chandra HETG/ACIS-S observation of EX Hya was performed between 2000 May 18 $`9^\mathrm{h}41^\mathrm{m}`$ and May 19 $`2^\mathrm{h}54^\mathrm{m}`$ UT for a total exposure of 59 ks. Extraction of the grating spectra and calculation of the effective area files was accomplished with the CIAO 2.1 suite of software using the reprocessed data products and new calibration data files (version R4CU5UPD8) for sequence 300041. Various aspects of this spectrum have been discussed by Mauche (2000, 2002) and Mauche, Liedahl, & Fournier (2001); here we discuss only a tiny portion ($`\lambda \lambda 11.60`$–12.05 Å) of the medium-energy grating (MEG) spectrum as it bears on the relative strengths of the Fe XXII $`n=32`$ emission lines. This spectrum is shown in Figure 3, where the ordinate is in units of $`\mathrm{counts}\mathrm{s}^1\mathrm{\AA }^1`$, the wavelength bin width $`\mathrm{\Delta }\lambda =0.005`$ Å, and we have combined $`\pm `$ first order counts. This portion of the spectrum contains the 11.77 and 11.92 Å emission lines of Fe XXII as well as the 11.74 Å emission line of Fe XXIII.

To determine the flux in these emission lines, we fit the MEG spectrum over the $`\lambda \lambda =11.60`$–12.05 Å wavelength interval with a model consisting of a linear background (to account for the thermal bremsstrahlung continuum obvious in the broadband spectrum) and three Gaussians. Assuming a common Gaussian width, the fitting function $`f(\lambda ;\stackrel{}{a})=a_1+a_2(\lambda \lambda _0)+_{i=0}^2a_{5+2i}\mathrm{exp}[(\lambda a_{4+2i})^2/2a_3^2]/\sqrt{2\pi }a_3`$. To perform the fit, we separately accounted for $`\pm `$ first orders and used $`\mathrm{\Delta }\lambda =0.005`$ Å bins, so there were 180 data points and 171 degrees of freedom (dof). The resulting fit, combining $`\pm `$ first orders, is shown by the thick histogram in Figure 3. It gives $`\chi ^2/\mathrm{dof}=104/171=0.60`$ and fit parameters as follows: Gaussian width $`a_3=1.4_{0.9}^{+1.6}`$ mÅ and line wavelengths and fluxes $`I(11.740\mathrm{\AA })=7.74(\pm 0.80)\times 10^5`$, $`I(11.776\mathrm{\AA })=4.27(\pm 0.67)\times 10^5`$, and $`I(11.923\mathrm{\AA })=4.52(\pm 0.68)\times 10^5\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1`$, where the errors on the wavelengths are $`\pm `$ 1–2 mÅ and all errors are $`1\sigma `$ (68% confidence) for one interesting degree of freedom. Formally (e.g., Bevington 1960), the $`I(11.92\mathrm{\AA })/I(11.78\mathrm{\AA })`$ line ratio equals $`1.06\pm 0.23`$ (68% confidence), but this propagation of errors ignores a subtle interplay between the error bars of the numerator and denominator of this ratio. To generate more rigorous values for the error bars of this ratio, we used a Monte Carlo technique with $`10^6`$ realizations to determine that the observed $`I(11.92\mathrm{\AA })/I(11.78\mathrm{\AA })`$ line ratio equals 1.06 and the errors are $`\genfrac{}{}{0pt}{}{+0.26}{0.20}`$ (68% confidence), $`\genfrac{}{}{0pt}{}{+0.48}{0.32}`$ (90% confidence), and $`\genfrac{}{}{0pt}{}{+0.89}{0.47}`$ (99% confidence). This best-fit line ratio and these error envelopes are shown superposed on both panels of Figure 2.

## 4 Analysis

Figure 2 shows that the combination of our LXSS models and Chandra HETG spectrum places significant constraints on the density of the $`T_\mathrm{e}12`$ MK plasma in EX Hya. Specifically, allowing for a factor-of-two uncertainty in the electron temperature, the upper panel of Figure 2 shows that the electron density $`n_\mathrm{e}=1.0_{0.5}^{+2.0}\times 10^{14}\mathrm{cm}^3`$ at 68% confidence, $`n_\mathrm{e}3.2\times 10^{13}\mathrm{cm}^3`$ at 90% confidence, and $`n_\mathrm{e}1.7\times 10^{13}\mathrm{cm}^3`$ at 99% confidence. Conversely, the lower panel of Figure 2 shows that the electron density can be arbitrarily low if the photoexcitation temperature $`T_{\mathrm{bb}}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}100`$ kK, but that alternative is excluded by UV observations of EX Hya. First, assuming that the distance to EX Hya is 100 pc, that the radius of its white dwarf is $`R_{\mathrm{wd}}=10^9`$ cm, that its 1010 Å flux density equals $`2.5\times 10^{13}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1\mathrm{\AA }^1`$ (Mauche, 1999), and that 100% of this emission is due to a 100 kK hotspot (i.e., assuming that the white dwarf, accretion disk, and accretion curtains do not radiate in the UV), the fractional emitting area of the spot is approximately 0.5%, far smaller than is inferred from optical and UV light curves of mCVs. Second, fits to EX Hya’s UV spectrum yield effective temperatures $`T_{\mathrm{eff}}25`$ kK (Greeley et al., 1997; Eisenbart et al., 2002), far lower than that required to produce the anomalous Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio via photoexcitation.

## 5 Discussion

To date, four different spectroscopic diagnostics have provided evidence of high densities in the X-ray–emitting plasma of the intermediate polar EX Hya. Hurwitz et al. (1997) used the line ratio of the Fe XX/Fe XXIII 133 Å blend to the Fe XXI 129 Å line observed in the 1994 Extreme Ultraviolet Explorer spectrum of EX Hya to infer $`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}10^{13}\mathrm{cm}^3`$ under the assumption that the plasma temperature $`T_\mathrm{e}=10`$ MK. Using the Chandra HETG spectrum of EX Hya, Mauche (2002) showed that the He-like ion $`R`$ line ratios of O, Ne, Mg, Si, and S are all in their high-density limit; Mauche, Liedahl, & Fournier (2001) used the Fe XVII $`I(17.10\mathrm{\AA })/I(17.05\mathrm{\AA })`$ line ratio to infer $`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}2\times 10^{14}\mathrm{cm}^3`$ for its $`T_\mathrm{e}4`$ MK plasma; and we here used the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio to infer $`n_\mathrm{e}1\times 10^{14}\mathrm{cm}^3`$ for its $`T_\mathrm{e}12`$ MK plasma. Of these diagnostics, the Fe XXII line ratio is the most reliable because it has the highest critical density and is the least sensitive to temperature and photoexcitation. We conclude first that the Chandra HETG spectrum of EX Hya requires plasma densities that are orders of magnitude greater than those typically observed in the Sun or other late-type stars, and second that the Fe XXII $`I(11.92\mathrm{\AA })/I(11.77\mathrm{\AA })`$ line ratio \[like the Fe XVII $`I(17.10\mathrm{\AA })/I(17.05\mathrm{\AA })`$ line ratio\] is a useful density diagnostic for sources like mCVs in which the plasma densities are high and the efficacy of the He-like ion density diagnostic is compromised by the presence of a bright UV continuum. Finally, we note that our density determinations for EX Hya are consistent with those expected for plasma in the accretion column, where $`n_\mathrm{e}\genfrac{}{}{0pt}{}{\text{ }\text{ }}{\text{ }\text{ }}T_\mathrm{e}^1`$. Additional density diagnostics are needed to begin to map the temperature/density profile of the accretion column of this and other mCVs.

We thank H. Tananbaum for the generous grant of Director’s Discretionary Time that made possible the Chandra observations of EX Hya. Support for this work was provided in part by NASA through Chandra Award Number DD0-1004B issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.