text:
# Energies and collapse times of symmetric and symmetry-breaking states of finite systems with a U(1) symmetry

## Abstract

We study quantum systems of volume $`V`$, which will exhibit the breaking of a U(1) symmetry in the limit of $`V\mathrm{}`$, when $`V`$ is large but finite. We estimate the energy difference between the ‘symmetric ground state’ (SGS), which is the lowest-energy state that does not breaks the symmetry, and a ‘pure phase vacuum’ (PPV), which approaches a symmetry-breaking vacuum as $`V\mathrm{}`$. Under some natural postulates on the energy of the SGS, it is shown that PPVs always have a higher energy than the SGS, and we derive a lower bound of the excess energy. We argue that the lower bound is $`𝒪(V^0)`$, which becomes much larger than the excitation energies of low-lying excited states for a large $`V`$. We also discuss the collapse time of PPVs for interacting many bosons. It is shown that the wave function collapses in a microscopic time scale, because PPVs are not energy eigenstates. We show, however, that for PPVs the expectation value of any observable, which is a finite polynomial of boson operators and their derivatives, does not collapse for a macroscopic time scale. In this sense, the collapse time of PPVs is macroscopically long.

The symmetry breaking (SB) is a key to understand quantum systems of many degrees of freedom. Although SB is formally defined for infinite systems, the physics of SBs in finite systems have been attracting much attention  for the following reasons: (i) Progress of experimental techniques enables one to observe and examine phase transitions in small systems, such as small magnets, small superconductors , liquid Helium in a small bubble , and laser-trapped atoms , hence SBs in finite systems should be studied seriously. (ii) Although recent progress of computers enables one to obtain ground states of finite systems numerically, the relation between such ground states and SB ground states of infinite systems are non-trivial: In a mean-field approximation, ground states of a finite system of volume $`V`$ are degenerate (if it will exhibit a SB as $`V\mathrm{}`$), and each ground state approaches a SB vacuum of the infinite system as $`V\mathrm{}`$. We call such a state that has a finite expectation value of an order parameter and approaches (i.e., well approximates) a SB vacuum of the infinite system as $`V\mathrm{}`$ ‘a pure phase vacuum’ (PPV). On the other hand, if one diagonalizes the Hamiltonian of the finite system exactly, the energy spectrum is much different from that of a mean-field approximation. The ground states are not necessarily degenerate. Moreover, it often occurs that a symmetric state which does not break the symmetry is a ground state, whereas PPVs have higher energies , in the absence of a SB field, which is usually considered as an unphysical, artificial field for the breaking of a U(1) symmetry. We call such a ground state the ‘symmetric ground state’ (SGS). In contrast to PPVs, the SGS does not approach a SB vacuum of the infinite system as $`V\mathrm{}`$. Hence, for a SB to occur, the energy difference between PPVs and the SGS should be small enough. Although the magnitude of this energy difference has been attracting much attention , definite conclusions have not yet been reached for the breaking of a U(1) symmetry. For example, an exact calculation  gave only a rough estimate (see the discussion following Eq. (15)).

In this paper, we estimate much more strictly the energy difference between PPVs and the SGS for the breaking of a U(1) symmetry. Under some natural postulates on the energy of the SGS, we show that PPVs always have a higher energy than the SGS, and that the excess energy is lower-bounded by $`\mu ^{}\delta N^2/2V`$, where $`\mu ^{}`$ is the derivative of the chemical potential with respect to the density $`nN/V`$, and $`\delta N^2`$ denotes the fluctuation of $`N`$. We further argue that this lower bound is $`𝒪(V^0)`$, which becomes much larger than the excitation energies of low-lying excited states for a large $`V`$. This should be contrasted with the breaking of the Z2 symmetry, for which the energy difference between PPVs and the SGS is only $`𝒪(V^1)`$ , which becomes much smaller than the excitation energies in a three-dimensional space for a large $`V`$. We also study the collapse time $`t_{\mathrm{coll}}`$ of PPVs for the case of interacting many bosons. It is shown that $`t_{\mathrm{coll}}=𝒪(V^0)`$ for the wave function, because PPVs are not energy eigenstates. We show, however, that for PPVs $`t_{\mathrm{coll}}=𝒪(\sqrt{V})`$ for the expectation value of any observable, which is a finite polynomial of boson operators and their derivatives, if the degree of the polynomial is fixed independent of $`V`$. In this sense, the collapse time of PPVs is macroscopically long.

We consider a quantum system that has a U(1) symmetry, whose conserved charge $`\widehat{N}`$ has integral eigenvalues (in an appropriate unit). We assume that the system is uniform, with the periodic boundary conditions, in order to eliminate additional complexities caused by non-uniform potentials or surface effects. Since the system volume $`V`$ is finite, $`\widehat{N}`$ is always well-defined, hence there exist simultaneous eigenstates $`|N,\mathrm{}`$ of $`\widehat{H}`$ and $`\widehat{N}`$;

$`\widehat{H}|N,\mathrm{}`$ $`=`$ $`E_{N,\mathrm{}}|N,\mathrm{},`$ (1)
$`\widehat{N}|N,\mathrm{}`$ $`=`$ $`N|N,\mathrm{},`$ (2)
where $`\mathrm{}`$ is a (set of) quantum number(s). For each value of $`N`$, there exists the lowest-energy state $`|N,\mathrm{G}`$, which we assume is non-degenerate. In general, $`\widehat{N}`$ becomes the generator of the U(1) symmetry. Hence, $`|N,\mathrm{G}`$ is the SGS. We now make two basic postulates: postulate 1 is the extensivity of the lowest eigenenergy;

$$E_{N,\mathrm{G}}=Vϵ(N/V),$$
(3)
where $`ϵ`$ is a single-variable function of the charge density $`N/V`$. This postulate seems natural under our assumptions that the system is uniform and $`|N,\mathrm{G}`$ is non-degenerate, as long as the charge does not induce a long-range force. \[Hence, care must be taken when the present results are applied to systems for which $`N`$ is the electric charge.\] By taking the limit $`V\mathrm{}`$, we can define $`ϵ(n)`$ for every continuous value of $`n`$ ($`=N/V`$). Using this $`ϵ(n)`$ for a finite $`V`$ as well, we can regard $`N`$ in Eq. (3) as a continuous variable. Hence, we can define $`\mu (n)ϵ^{}(n)=\frac{}{N}E_{N,\mathrm{G}}.`$ Postulate 2 is

$$\mu ^{}(n)ϵ^{\prime \prime }(n)=V\frac{^2}{N^2}E_{N,\mathrm{G}}>0,$$
(4)
which also seems natural because thermodynamics requires that $`\mu `$ should be a non-decreasing function of $`n`$, for the system to be stable. Although $`\mu ^{}=0`$ for non-interacting particles (such as a photon gas whose $`\mu `$ is always zero), we assume $`\mu ^{}>0`$ because we are not interested in such a trivial case. For weakly-interacting many bosons with a repulsive interaction (with the effective coupling constant $`g>0`$), for example, these postulates are indeed satisfied; $`E_{N,\mathrm{G}}=Vϵ(N/V)`$, where $`ϵ(n)=gn^2/2+\mathrm{},`$ hence $`\mu ^{}=g+\mathrm{}>0.`$

When the system exhibits the breaking of the U(1) symmetry in the limit of $`V\mathrm{}`$ (while keeping the charge density finite), a SB state cannot be $`lim_V\mathrm{}|N,\mathrm{G}`$ or $`lim_V\mathrm{}|N,\mathrm{}`$, because they are eigenstates of $`\widehat{N}`$ . Therefore, one should take superpositions of states with different charges in order to construct a PPV, which approaches (i.e, well approximates) a SB vacuum as $`V\mathrm{}`$ . If the quantum system (of a finite $`V`$) is perfectly closed, such superpositions are forbidden for massive and/or charged particles by the charge superselection rule, which requires that any pure state must be an eigenstate of $`\widehat{N}`$ . However, if the quantum system is a subsystem of a larger system , we previously showed that one can associate a pure state, which is a coherent superposition of states with different charges, to the subsystem . We investigate energy expectation values of such states, which in general have either of the following forms;

$`|C`$ $``$ $`{\displaystyle \underset{N}{}}C_N|N,\mathrm{G},`$ (5)
$`|\stackrel{~}{C}`$ $``$ $`{\displaystyle \underset{N,\mathrm{}}{}}\stackrel{~}{C}_{N,\mathrm{}}|N,\mathrm{},`$ (6)
where $`C_N`$ and $`\stackrel{~}{C}_{N,\mathrm{}}`$ are coefficients, which are normalized as $`_N|C_N|^2=1`$ and $`_{N,\mathrm{}}|\stackrel{~}{C}_{N,\mathrm{}}|^2=1`$, respectively. We are interested in the case where the charge density $`n`$ ($`=N/V`$) approaches a finite value as $`V\mathrm{}`$. (On the other hand, the case $`n0`$ is rather trivial, while the case $`n\mathrm{}`$ is anomalous.) Namely,

$$N=𝒪(V),$$
(7)
where $`NC|\widehat{N}|C`$ or $`\stackrel{~}{C}|\widehat{N}|\stackrel{~}{C}`$, and $`n=N/V`$. Furthermore, PPVs should be consistent with thermodynamics, according to which variances of extensive variables are of $`𝒪(V)`$ or smaller. Hence,

$$\delta N^2𝒪(V),$$
(8)
where $`\delta N^2C|\delta \widehat{N}^2|C`$ or $`\stackrel{~}{C}|\delta \widehat{N}^2|\stackrel{~}{C}`$, where $`\delta \widehat{N}\widehat{N}N`$. Namely, $`C_N`$ and $`\stackrel{~}{C}_{N,\mathrm{}}`$ should localize in the $`N`$ space in such a way that Eqs. (7) and (8) are satisfied. Note that the phases of $`C_N`$ and $`\stackrel{~}{C}_{N,\mathrm{}}`$ are irrelevant to $`N`$, $`\delta N^2`$, and the energy expectation values, $`C|\widehat{H}|C=_N|C_N|^2E_{N,\mathrm{G}}`$ and $`\stackrel{~}{C}|\widehat{H}|\stackrel{~}{C}=_{N,\mathrm{}}|C_N|^2E_{N,\mathrm{}}`$.

For $`N=N+\mathrm{\Delta }N`$, where $`(\mathrm{\Delta }N)^2𝒪(V)`$, Eq. (3) can be expanded as

$$E_{N,\mathrm{G}}=V\left[ϵ\left(\frac{N}{V}\right)+\frac{\mathrm{\Delta }N}{V}\mu \left(\frac{N}{V}\right)+\frac{1}{2}\left(\frac{\mathrm{\Delta }N}{V}\right)^2\mu ^{}\left(\frac{N}{V}\right)+𝒪\left(\frac{1}{V^{3/2}}\right)\right].$$
(9)
We neglect the higher-order term, $`V𝒪(1/V^{3/2})`$, in the following analysis. Then, we can easily show that

$$C|\widehat{H}|C=E_{N,\mathrm{G}}+\frac{\delta N^2}{2V}\mu ^{}\left(\frac{N}{V}\right).$$
(10)
Since the last term is positive because of postulate 2, we conclude that $`|C`$ always has a higher energy than $`|N,\mathrm{G}`$ if they have the same value of $`N`$ . Here, it is crucial to fix the value of $`N`$ for the comparison. \[Otherwise, either state could have a higher energy because of the linear term in Eq. (9).\] Note that formula (10), although very simple, gives the energy of a general state of the form (5) very precisely, with the error being only $`𝒪(1/V^{1/2})`$. Note also that the energy expectation value is determined only by $`N`$ and $`\delta N^2`$ if the functional forms of $`E_{N,\mathrm{G}}`$ and $`\mu ^{}(n)`$ are given.

For a more general state (6), we derive an inequality. Let $`C_N^{}`$’s be some coefficients which satisfy $`|C_N^{}|^2=_{\mathrm{}}|\stackrel{~}{C}_{N,\mathrm{}}|^2`$. For $`|C^{}_NC_N^{}|N,\mathrm{G}`$,

$$\stackrel{~}{C}|\widehat{H}|\stackrel{~}{C}C^{}|\widehat{H}|C^{}=\underset{N,\mathrm{}}{}|\stackrel{~}{C}_{N,\mathrm{}}|^2(E_{N,\mathrm{}}E_{N,\mathrm{G}})0,$$
(11)
where the equality holds iff $`\stackrel{~}{C}_{N,\mathrm{}}=0`$ for all $`\mathrm{}\mathrm{G}`$. Applying Eq. (10) to $`C^{}|\widehat{H}|C^{}`$, and noting that $`\stackrel{~}{C}|\widehat{N}|\stackrel{~}{C}=C^{}|\widehat{N}|C^{}`$ and $`\stackrel{~}{C}|\delta \widehat{N}^2|\stackrel{~}{C}=C^{}|\delta \widehat{N}^2|C^{}`$, we obtain

$$\stackrel{~}{C}|\widehat{H}|\stackrel{~}{C}E_{N,\mathrm{G}}+\frac{\delta N^2}{2V}\mu ^{}\left(\frac{N}{V}\right).$$
(12)
When $`|\stackrel{~}{C}`$ and $`|C`$ have the same values of $`N`$ and $`\delta N^2`$, Eqs. (10) and (12) can be combined as the simple formula;

$$\stackrel{~}{C}|\widehat{H}|\stackrel{~}{C}C|\widehat{H}|C=E_{N,\mathrm{G}}+\frac{\delta N^2}{2V}\mu ^{}\left(\frac{N}{V}\right).$$
(13)
It is easy to show the similar result for $`\widehat{K}\widehat{H}\mu \widehat{N}`$ (where $`\mu `$ here is a constant);

$$\stackrel{~}{C}|\widehat{K}|\stackrel{~}{C}C|\widehat{K}|C=K_{N,\mathrm{G}}+\frac{\delta N^2}{2V}ϵ^{\prime \prime }\left(\frac{N}{V}\right),$$
(14)
where $`K_{N,\mathrm{G}}E_{N,\mathrm{G}}\mu N`$. Hence, the following results are applicable to the expectation values of $`\widehat{K}`$ as well, if we replace $`K_{N,\mathrm{G}}`$ and $`\mu ^{}`$ with $`E_{N,\mathrm{G}}`$ and $`ϵ^{\prime \prime }`$, respectively.

Since a PPV should take either form (5) or (6), we conclude that (i) a PPV always has a higher energy than the SGS, and (ii) the excess energy is lower-bounded by $`\mu ^{}\delta N^2/2V`$. This is the first of the main results of this paper. It should be mentioned that Ref.  tried to give an upper bound of the energy difference between the SGS and ‘low-lying states,’ a linear combination of which is a PPV, while our result gives a lower bound of the energy increase of PPVs over the SGS.

We now estimate how the lower bound $`\mu ^{}\delta N^2/2V`$ behaves with increasing $`V`$. For interacting many-bosons, we previously found the state vector of a PPV, which we called the coherent state of interacting bosons (CSIB) . This state vector, denoted by $`|\alpha ,\mathrm{G}`$, has the form of Eq. (5) with $`C_N=e^{|\alpha |^2}\alpha ^N/\sqrt{N!}.`$ Hence, $`\delta N^2=N`$, and Eq. (10) yields

$$\alpha ,\mathrm{G}|\widehat{H}|\alpha ,\mathrm{G}E_{N,\mathrm{G}}=(n/2)\mu ^{}(n)=𝒪(V^0)>0,$$
(15)
where the sign is determined by postulate 2 ($`\mu ^{}>0`$), and $`n=N/V`$. On the other hand, if we apply the inequality of Ref.  to the CSIB, we obtain $`|\alpha ,\mathrm{G}|\widehat{H}|\alpha ,\mathrm{G}E_{N,\mathrm{G}}|𝒪(\sqrt{V})`$, which diverges as $`V\mathrm{}`$. Although there is no contradiction between the two results, the present result gives a much more accurate estimate: the CSIB has a higher energy than the SGS by $`𝒪(V^0)`$, for the same value of $`N`$. Although one might expect that the energy increase would be a decreasing function of $`V`$ (as in the case of the breaking of the Z2 symmetry ), our result denies such a naive expectation.

For general systems which exhibit the breaking of a U(1) symmetry, we do not know the explicit forms of PPVs. By virtue of relation (13), however, it is sufficient to estimate $`\delta N^2`$ for the estimation of the energy increase. We argue that

$$\delta N^2N$$
(16)
for general systems which exhibit the breaking of a U(1) symmetry, for the following reasons. Macroscopic properties of PPVs must be stable against weak perturbations from environments. The environments include those which exchange charges with the system. We may apply the classical stochastic theory to estimate the stable distribution of the charges because (i) the phases of $`C_N`$ and $`\stackrel{~}{C}_{N,\mathrm{}}`$ are irrelevant to $`N`$ and $`\delta N^2`$, and (ii) the phase coherence between the environments and the system may be negligible if the dephasing times of the environments are short enough. Then, according to the classical stochastic theory, the steady-state distribution of the charges should satisfy Eq. (16) when charges are randomly exchanged with a huge environment.

From Eqs. (4), (7), (13) and (16), we obtain

$$\mathrm{PPV}|\widehat{H}|\mathrm{PPV}E_{N,\mathrm{G}}𝒪(V^0)>0,$$
(17)
for general systems which exhibit the breaking of a U(1) symmetry. Note that this result is consistent with the theory of SBs in infinite systems, according to which PPVs have the same energy density as the SGS . In fact, Eq. (17) yields $`(\mathrm{PPV}|\widehat{H}|\mathrm{PPV}E_{N,\mathrm{G}})/V𝒪(V^1)0`$ as $`V\mathrm{}`$, for the lower bound of the difference in the energy densities. On the other hand, our result denies a naive expectation that the energies of PPVs and the SGS would be ‘almost degenerate’ in the sense that the energy difference would be a decreasing function of $`V`$. Furthermore, the energy difference for a large $`V`$ becomes much larger than the excitation energies of low-lying excited states, whose wavenumber $`kV^{1/d}`$ in a $`d`$-dimensional space, because the excitation energy $`ϵ(k)`$ behaves as $`ϵ(k)|k|V^{1/d}`$ for a linear dispersion, and $`ϵ(k)|k|^2=V^{2/d}`$ for a parabolic dispersion. This should be contrasted with the breaking of the Z2 symmetry, for which the energy difference between PPVs and the SGS is only $`𝒪(V^1)`$ , which becomes much smaller than the excitation energies in a three-dimensional space for a large $`V`$. This indicates, for example, that much more care is necessary for the breaking of the U(1) symmetry than for the breaking of the Z2 symmetry, when one tries to find a PPV by numerical calculations.

We finally discuss the collapse time of PPVs, by generalizing the discussion of Ref. . In general, PPVs are not an energy eigenstate, hence their wave functions deform in finite systems as time evolves. Let $`t_{\mathrm{coll}}`$ be the collapse time, which is defined as the time scale at which this deformation becomes significant. For example, for an initial ($`t=0`$) state of the form of Eq. (5), it evolves as $`|C_NC_Ne^{iE_{N,\mathrm{G}}t}|N,\mathrm{G}|C;t,`$ where $`\mathrm{}`$ is taken unity. Since $`E_{N+\mathrm{\Delta }N,\mathrm{G}}E_{N,\mathrm{G}}=\mu \mathrm{\Delta }N+\mathrm{}`$ from Eq. (9), the difference of $`|C;t`$ from $`|C`$ becomes significant at $`t=t_{\mathrm{coll}}1/\mu \sqrt{\delta N^2}`$, because the linear term $`\mu \mathrm{\Delta }N`$ alters the relative phases among $`C_N`$’s, except when $`C_N`$’s take some special forms. As $`V`$ is increased, this time scale approaches zero if $`\sqrt{\delta N^2}`$ increases in proportion to $`V`$ as Eq. (16). On the other hand, PPVs must survive over a macroscopic time scale, i.e., $`t_{\mathrm{coll}}\mathrm{}`$ as $`V\mathrm{}`$ for PPVs. To satisfy this condition, $`C_N`$’s of PPVs must take some special forms. For interacting many-bosons, for example, $`C_N`$’s of the CSIB indeed have special forms, for which the effect of the linear term $`\mu \mathrm{\Delta }N`$ on $`|\alpha ,\mathrm{G};t`$ is completely absorbed as a time evolution of the single parameter $`\alpha `$. In fact,

$`|\alpha ,\mathrm{G};t`$ $`=`$ $`e^{|\alpha |^2}{\displaystyle \underset{N}{}}{\displaystyle \frac{\alpha ^N}{\sqrt{N!}}}e^{i[E_{N,\mathrm{G}}+\mu \mathrm{\Delta }N+\mu ^{}(\mathrm{\Delta }N)^2/2V]t}|N,\mathrm{G}`$ (18)
$`=`$ $`e^{i(E_{N,\mathrm{G}}\mu N)t}e^{|\alpha |^2}{\displaystyle \underset{N}{}}{\displaystyle \frac{(\alpha e^{i\mu t})^N}{\sqrt{N!}}}e^{i[\mu ^{}(\mathrm{\Delta }N)^2/2V]t}|N,\mathrm{G}.`$ (19)
Since the prefactor $`e^{i(E_{N,\mathrm{G}}\mu N)t}`$ has no physical meaning, we find that $`|\alpha ,\mathrm{G};t=|\alpha e^{i\mu t},\mathrm{G}`$ if $`\mu ^{}=0`$. Namely, the CSIB does not collapse at all if $`\mu ^{}=0`$; only the phase of $`\alpha `$ evolves with time as $`\alpha e^{i\mu t}`$. This result for $`\mu ^{}=0`$ is well-known. If $`\mu ^{}>0`$, on the other hand, the wave function $`|\alpha ,\mathrm{G}`$ collapses at $`tV/\mu ^{}(\mathrm{\Delta }N)^2V/\mu ^{}\delta N^2=𝒪(V^0)`$. However, this does not necessarily mean that the expectation values of observables of interest alter in this time scale. For example, if an observable is proportional to the boson operator $`\widehat{\psi }`$ or its derivative, it detects the phase relation between adjacent coefficients, $`C_{N+1}`$ and $`C_N`$. The ratio of their phases evolves, for $`N=N+\mathrm{\Delta }N`$, as

$$\frac{(e^{i\mu t})^{N+1}e^{i\mu ^{}(\mathrm{\Delta }N+1)^2t/2V}}{(e^{i\mu t})^Ne^{i\mu ^{}(\mathrm{\Delta }N)^2t/2V}}=e^{i\mu t}e^{i\mu ^{}t/2V}e^{i\mu ^{}\mathrm{\Delta }Nt/V}.$$
(20)
In the right-hand side, the first factor $`e^{i\mu t}`$ can be absorbed as the time evolution of the single parameter $`\alpha \alpha e^{i\mu t}`$, whereas the second factor $`e^{i\mu ^{}t/2V}`$ is negligible because $`\mu ^{}/2V=𝒪(1/V)`$. Hence, only the last factor $`e^{i\mu ^{}\mathrm{\Delta }Nt/V}`$ is relevant to the collapse time. We thus find

$$t_{\mathrm{coll}}V/(\mu ^{}\mathrm{\Delta }N)V/(\mu ^{}\sqrt{\delta N^2})=𝒪(\sqrt{V}).$$
(21)
For general observables which are polynomials of degree $`M`$ of $`\widehat{\psi }`$, $`\widehat{\psi }^{}`$ and their derivatives, we obtain

$$t_{\mathrm{coll}}=𝒪(\sqrt{V}/M).$$
(22)
Hence, if the degree $`M`$ of the polynomial is fixed independent of $`V`$, we again obtain $`t_{\mathrm{coll}}=𝒪(\sqrt{V})`$. Since the expectation values of observables are relevant in quantum theory, we may conclude that the collapse time of the CSIB is $`𝒪(\sqrt{V})`$, which is macroscopic in the sense that it diverges as $`V\mathrm{}`$, although the collapse time of its wave function is $`𝒪(V^0)`$ . If, on the other hand, if $`M`$ is increased in proportion to $`\sqrt{V}`$, then $`t_{\mathrm{coll}}=𝒪(V^0)`$. Except for such an abnormal case (as $`M\sqrt{V}`$), $`t_{\mathrm{coll}}`$ is macroscopic. For more general systems with the breaking of U(1) symmetry, we have not yet obtained definite conclusions on $`t_{\mathrm{coll}}`$, although we expect a situation similar to the case of interacting many bosons. This may be a subject of future studies.


format: text

text:
# Exploiting VSIPL and OpenMP for Parallel Image Processing

## 1. Introduction

The Vector, Signal and Image Processing Library (VSIPL)  is an open standard C language Application Programmer Interface (API) that allows portable and optimized single processor programs. OpenMP  is an open standard C/Fortran API that allows portable thread based parallelism on shared memory computers. Both of these standards have enormous potential to allow users to realize the goal of portable applications that are both parallel and optimized.

Exploiting these new open standards requires integrating them into existing applications as well as using them in new efforts. Image processing is one of the key areas where VSIPL and OpenMP can make a large impact. Currently, a large fraction of image processing applications are written in the Interpreted Data Language (IDL) environment . The goal of this work is to show that it is possible to bring the performance benefits of these new standards to the image processing community in a high level manner that is transparent to users.

## 2. Approach

Wide area 2D convolution is a staple of digital image processing (see Figure 1). The advent of large format CCDs makes it possible to literally “pave” with silicon the focal plane of an optical sensor. Processing of the large images obtained from these systems is complicated by the non-uniform Point Response Function (PRF) that is common in wide field of view instruments. This paper presents a fast, FFT based algorithm for convolving such images. This algorithm has been transparently implemented within IDL environment using VSIPL (for optimized single processor performance) with added OpenMP directives (for parallelism).

The inputs of image convolution with variable PRFs consists of a source image, a set of PRF images, and a grid which locates the center of each PRF on the source image. The output image is the convolution of the input image with each PRF linearly weighted by its distance from its grid center. The computational basis of this convolution are 2D overlap and add FFTs with interpolation (see Figure 2). Today, typical images sizes are in the millions (2K x 2K) to billions (40K x 40K) of pixels. A single PRF is typically thousands of pixels (100 x 100) pixels, but can be as small 10 x 10 or as large as the entire image. Over a single image a PRF will be sampled as few as once but as many as hundreds of times depending on the optical system.

There are many opportunities for parallelism in this algorithm. The simplest is to convolve each PRF separately on a different processor and then combine all the results on a single processor. This approach works well with VSIPL, OpenMP and IDL (see Figure 3. At the top level a user passes the inputs into an IDL routine which passes pointers to an external C function. Within the C function OpenMP forks off multiple threads. Each thread executes its convolution using VSIPL functions. The OpenMP threads are then rejoined and the results are added. Finally a pointer to the output image is returned to the IDL environment in the same manner as any other IDL routine.

## 3. Results

This algorithm was implemented on an SGI Origin 2000 at Boston University . This machine consists of 64 300 MHz MIPS 10000 processors with an aggregate memory of 16 GBytes. IDL version 5.3 from Research Systems, Inc. was used along with SGI’s native OpenMP compiler (version 7.3.1) and the TASP VSIPL implementation. Implementing the components of the system was the same as if each were done separately. Integrating the pieces (IDL/OpenMP/VSIPL) was done quickly, although care must be taken to use the latest versions of the compilers and libraries. Once implemented the software can be quickly ported via Makefile modifications to any system that has IDL, OpenMP, and VSIPL (currently these are SGI, HP, Sun, IBM, and Red Hat Linux). We have conducted a variety of experiments which show linear speedups using different numbers of processors and different image sizes (see Figure 4). Thus, it possible to achieve good performance using open standards underneath existing high level languages.

## References

 http://www.vsipl.org/

 http://www.openmp.org/

 http://www.rsi.com/

 http://scv.bu.edu/


text:
# On the nature of spinor Bose-Einstein condensates in rubidium

## Abstract

We perform detailed close-coupling calculations for the rubidium isotopes 85Rb and 87Rb to ascertain the nature of their spinor Bose-Einstein condensates. These calculations predict that the spinor condensate for the spin-1 boson 87Rb has a ferromagnetic nature. The spinor condensates for the spin-2 bosons 85Rb and 87Rb, however, are both predicted to be polar. The nature of a spin-1 condensate hinges critically on the sign of the difference between the s-wave scattering lengths for total spin 0 and 2 while the nature of a spin-2 condensate depends on the values of the differences between s-wave scattering lengths for the total spin 0, 2 and 4. These scattering lengths were extracted previously and found to have overlapping uncertainties for all three cases, thus leaving the nature of the spinor condensates ambiguous. The present study exploits a refined uncertainty analysis of the scattering lengths based on recently improved result from experimental work by Roberts et al., which permits us to extract an unambiguous result for the nature of the ground state spinor condensates.

In a conventional magnetic trap for ultra-cold alkali atoms the spin degrees of freedom are “frozen out” since the atom must be in a weak-field seeking Zeeman state to be trapped. In an optical trap, however, the spins of the alkali atoms are essentially free, and all magnetic substates $`|f,m`$ for a given spin $`f`$ can be populated. Since the atom-atom interaction depends on spin, these magnetic substates can be changed in a scattering event. Accordingly, it is of interest to see how the spins are organized in the ground state and to explore the nature of the spin-mixing dynamics in an optically trapped Bose-Einstein condensate.

Multi-component condensates have been formed in magnetic traps. For instance Ref. Myatt:97  used a double magneto-optical trap and a magnetic trap to create condensates in either the $`|f=2,m=2`$ or the $`|f=1,m=1`$ spin state of 87Rb, and in a mixture of both by cooling $`|1,1`$ evaporatively and $`|2,2`$ via thermal contact with the $`|1,1`$ atoms. In this case the spin projections are approximately frozen out because the spin-flip cross sections in 87Rb are anomalously small Burke:97 ; Julienne:97 ; Kokkelmans:97 . By contrast, Ref. Stamper-Kurn:98  made a sodium condensate consisting simultaneously of all three magnetic substates of the $`f=1`$ atomic state, by cooling the atoms in a magnetic trap and then transferring them into a optical trap. This experimental technique produces what is referred to as a spinor condensate, because it can explore its full range of spin degrees of freedom. In the theoretical description of Ref. Ho:98 ; Ho:00 , the spinor condensates are classified according to the relative values of certain characteristic scattering lengths. Note that alternative theoretical treatments Law:98 ; Pu:00  differ in their detailed predictions concerning the nature of the spinor BEC ground state. Nevertheless, in this paper we determine the interaction parameters for spinor condensates of 85Rb and 87Rb which based on Ref. Ho:98 ; Ho:00 ; Ciobanu:00  fall into the two following categories:

Spin-1 atoms (87Rb) Let $`F`$ be the total spin of two bosonic spin $`f=1`$ atoms, and let $`a_F`$ be the s-wave scattering length for the total spin $`F`$ symmetry. Since $`f=1`$, only $`F=0,2`$ are allowed by Bose symmetry for an s-wave collision. The nature of the spin-1 BEC ground state depends critically on the relative values of $`a_0`$ and $`a_2`$. According to Ho Ho:98  a spinor Bose condensate composed of spin-1 bosons in an optical trap can be either “ferromagnetic” or “antiferromagnetic” in nature Ho:98 ; Ho:00 . The antiferromagnetic state has alternatively been termed “polar”, and we use this terminology here. The difference between the scattering lengths $`a_0`$ and $`a_2`$ determines the nature of the spin-1 condensate: the ferromagnetic state emerges when $`a_0>a_2`$, whereas the polar state emerges when $`a_0<a_2`$ Ho:98 . In the ferromagnetic state virtually all atoms reside in the same spin substate (either $`m=1`$ or $`m=1`$); in the polar state the spin projections are mixed.

Spin-2 atoms (85Rb, 87Rb) Two bosonic spin $`f=2`$ atoms possess $`F=0`$, $`2`$, $`4`$ total spin states exhibiting the appropriate Bose symmetry for an s-wave collision. For spin-2 87Rb the scattering lengths $`a_0`$, $`a_2`$ and $`a_4`$ are determined by the real part of the phase-shift since the inelastic scattering processes are also allowed. According to Ciobanu et al. Ciobanu:00 , a spinor condensate of spin-2 bosons in an optical trap can be one of the three types “ferromagnetic”, “polar”, or “cyclic” in nature, which we abbreviate as F, P or C respectively. Ferromagnetic and polar condensates are similar to those above. The name “cyclic” arises from a close analogy with d-wave BCS superfluids. The nature of the spin-2 BEC ground state depends critically on the relative values of $`a_0a_2`$ and $`a_2a_4`$ Ciobanu:00 .

The three states emerge under following conditions:

P: $`a_0a_4<0`$, $`\frac{2}{7}a_2a_4<\frac{1}{5}a_0a_4`$,

F: $`a_2a_4>0`$, $`\frac{1}{5}a_0a_4+\frac{2}{7}(a_2a_4)>0`$,

C: $`a_2a_4<0`$, $`\frac{1}{5}a_0a_4\frac{2}{7}(a_2a_4)>0`$.

For spin-1 87Rb the total spin $`F=0`$ and $`F=2`$ scattering lengths $`a_0`$ and $`a_2`$ are almost equal. They have been calculated before Ho:98  based on the analysis of Ref. Roberts:98 , but the uncertainties determined still overlap for $`a_0`$ and $`a_2`$, so that the sign of the difference has remained uncertain. In particular, the scattering lengths have been interpreted rather conservatively in Ref. Roberts:98 . For spin-2 85Rb and 87Rb the uncertainties for the total spin $`F`$ scattering lengths $`a_0`$, $`a_2`$ and $`a_4`$ have been too large to uniquely identify the nature of the spinor condensates Ho:00 . The uncertainty region for 85Rb was large enough to overlap all three regions P, F and C, while the uncertainty region for 87Rb overlapped both the polar and the cyclic region. In the present study we determine the scattering lengths $`a_0`$ and $`a_2`$, and their uncertainties for spin-1 87Rb, and $`a_0`$, $`a_2`$ and $`a_4`$ and their uncertainties for spin-2 85Rb and 87Rb. We concentrate on an accurate determination of the difference $`a_0a_2`$ for spin-1 87Rb and the pair ($`a_0a_4`$, $`a_2a_4`$) for spin-2 85Rb and 87Rb. If one accepts the spinor condensate treatment of Ref. Law:98 ; Pu:00  this analysis gives an unambiguous determination of the nature of the BEC ground states.

Uncertainties in the scattering lengths arise primarily from imperfect knowledge of three parameters: the long-range van der Waals coefficient $`C_6`$, and the singlet and the triplet s-wave scattering lengths, $`a_s`$ and $`a_t`$ respectively. In addition, when using potential curves determined for one isotope to predict scattering for another isotope, the results can depend on the precise number of bound states in the triplet potential, $`N_b`$ as well as the precise number of bound states in the singlet potential. Roberts et al. analyzed a magnetic-field Feshbach resonance to determine “state of the art” potentials for 85Rb Roberts:98 . Recently they have revisited some of the rethermalization measurements in Ref. Roberts:98  and improved the uncertainties for the long-range van der Waals coefficient and the singlet end triplet s-wave scattering lengths for 85Rb Roberts:01 .

Using these new values of $`C_6`$, $`a_s`$ and $`a_t`$ we show below unambiguously that $`a_0>a_2`$ for spin-1 87Rb. This result in turn implies that the spinor condensate is definitely ferromagnetic, as was previously suspected Ho:98 . By contrast the spin-1 23Na scattering lengths, recently determined in Ref. Crubellier:99 , imply that a 23Na $`f=1`$ spinor BEC is polar, as has been suggested before Ho:98 . By extracting the scattering length from a spectroscopic experiment, Crubellier et al. found that, for 23Na, $`a_0=50.0\pm 1.6`$ $`a.u.`$ and $`a_2=55.0\pm 1.7`$ $`a.u.`$ Crubellier:99 . They calculated the scattering lengths for two values of the $`C_6`$ coefficient for 23Na and found that the influence of the $`C_6`$ value is very small (a 4% change in $`C_6`$ results in a variation in the scattering length of the order of 0.1%). Consequently, the analysis for 23Na Crubellier:99 , in conjunction with the present analysis for 87Rb, implies that both types of spin-1 condensates can be realized with the atoms used most frequently in BEC experiments (23Na and 87Rb).

The improved results for $`C_6`$, $`a_s`$ and $`a_t`$ also predict that $`a_0a_4<0`$, $`\frac{2}{7}a_2a_4<\frac{1}{5}a_0a_4`$, for both spin-2 85Rb and 87Rb. This result implies that both spin-2 85Rb and 87Rb will be polar. Previously, it was estimated that 87Rb would be polar, but that 87Rb would be cyclic Ciobanu:00 . This implies that the ground state for spin-2 85Rb and 87Rb will have the same nature as spin-2 23Na. Spin-2 23Na was already unambiguously classified since the uncertainties on differences between the relevant scattering lengths place $`a_0a_4`$ and $`a_2a_4`$ within the polar region Ciobanu:00 . The results for spin-1 87Rb and for spin-2 85Rb and 87Rb are summarized in Figs. 1, 2 and 3 respectively.

Our calculations start from the singlet and triplet Born-Oppenheimer potentials between two rubidium atoms that were calculated in Ref. Krauss:90 , where the sing-let potential is adjusted to have 125 bound states private . These potentials are matched smoothly at $`r=22.0`$ $`a.u.`$ to the standard long-range van der Waals potentials using the new value of the long-range coefficient $`C_6`$ inferred from the experiment in Ref. Roberts:98  and reanalyzed according to Ref. Roberts:01 , and using the $`C_8`$ and $`C_{10}`$ coefficients from the calculations of Ref. Marinescu:94 . The potentials are adjusted to match the scattering length by including short-range inner-wall corrections that are parameterized for each spin by: $`c\mathrm{arctan}((rr_{min})^2/(c_r))`$ for $`r<r_{min}`$. $`c_r`$ is a constant (the same order of magnitude as $`r_{min}`$; slightly different for the singlet and the triplet), the inner-wall parameters $`c`$ are of the order of $`10^5`$ to $`10^4`$ $`a.u.`$, $`r`$ is the separation between the two Rb atoms and $`r_{min}`$ is the separation for which the potential is minimal. The inner-wall parameters $`c`$ are varied over a range that reproduces the recently improved values of $`a_s`$ and $`a_t`$. The improved values of $`C_6`$ for rubidium, $`a_s`$ and $`a_t`$ for 85Rb are: $`C_6=4660\pm 20`$ $`a.u.`$, $`a_s=3650_{670}^{+1500}`$ $`a.u.`$ and $`a_t=332\pm 18`$ $`a.u.`$ Roberts:01 , while the calculations of Ref. Marinescu:94  determined that $`\overline{C}_8=550600`$ $`a.u.`$. These are the values we adopt in the present calculations. Our calculations here do not allow for variance in $`C_8`$. This is reasonable because the dependence of $`C_8`$ is one order of magnitude smaller than the dependence of $`C_6`$. Furthermore the number of bound states in the triplet potential was previously believed to be $`39\pm 1`$ private ; Tsai:97 , but more refined experimental analysis suggests that it is instead $`40N_b42`$ private ; Leo . The present calculations are done for $`N_b=39,40,41,42`$. The number of bound states in the singlet potential is not changed.

Our calculations have been carried out for three values of $`C_6`$ that span the empirical range ($`46404680`$ $`a.u.`$). These values are adequate since the quantities of interest vary smoothly with $`C_6`$ over the range of interest. We also tested the triplet potential for each one of the four relevant $`N_b`$ ($`39`$, $`40`$, $`41`$, $`42`$). For each value of $`N_b`$ we determine the values of the inner-wall corrections that correspond to the uncertainty range of 85Rb $`a_s`$ and $`a_t`$ for each of the three values of $`C_6`$. These calculations are carried out at zero magnetic field and 130 $`nK`$ since the given value of $`C_6`$, $`a_s`$ and $`a_t`$ are determined from collisions at this temperature Roberts:01 . The same potentials optimized for 85Rb have been used in our 87Rb calculations, except for an appropriate change in the reduced mass. Since $`N_b`$ in rubidium is unknown at present, and since we utilize the same potentials determined by the 85Rb singlet and triplet scattering length in our 87Rb calculations, we incorporate $`N_b`$ in our analysis of the uncertainties for 87Rb. The singlet and triplet potentials are used in multichannel calculations to compute the $`a_0`$ and $`a_2`$ for spin-1 87Rb and $`a_0`$, $`a_2`$ and $`a_4`$ for spin-2 85Rb and 87Rb (again at zero magnetic field and at $`130`$ $`nK`$), with fixed $`C_6`$ and $`N_b`$. The calculations are then repeated for each value of $`C_6`$ and $`N_b`$ with the corresponding new values of the inner-wall corrections. These calculations span the empirical $`C_6`$, $`N_b`$, $`a_s`$ and $`a_t`$ range, which permits us to extract the over-all uncertainty in the difference between the relevant scattering lengths, $`a_0a_2`$ for spin-1 and the two relevant differences ($`a_0a_4`$, $`a_2a_4`$) for spin-2.

As a confirmation, the same rubidium potentials have also been used to calculate the singlet and triplet s-wave scattering lengths $`a_s`$ and $`a_t`$ for 87Rb (at 130 $`nK`$ energy). These single-channel calculations have been repeated for the three values of $`C_6`$ and four values of $`C_6`$ that span the empirical range. This permits us to check whether $`a_s`$ and $`a_t`$ fall within the range of previous measured values. The single-channel triplet scattering lengths are found to be:

$`\begin{array}{cc}N_b& a_t\\ 39& 107\pm 1a.u.\\ 40& 103\pm 1a.u.\\ 41& 100\pm 1a.u.\\ 42& 97\pm 1a.u.\end{array}`$

and the single-channel singlet scattering length is found to be: $`a_s=91\pm 1`$ $`a.u.`$. The $`a_s`$ and $`a_t`$ values for $`N_b=39`$ are in good agreement with previous work Roberts:98 ; Burke:99 . As another confirmation, and since $`N_b`$ is unknown we have also calculated the 85Rb scattering length for $`f=2`$, $`m=2`$ at various magnetic fields to compare the obtained values with the values from Roberts et al.; Roberts:01 . This comparison shows good agreement. For each of the $`N_b`$, the scattering lengths obtained for the specific magnetic fields exhibits an uncertainty greater than the one given by Crubellier:99 .

The values for $`a_0`$ and $`a_2`$ for spin-1 87Rb,along with their uncertainties, are shown in Fig. 1. $`a_0`$ is always greater than $`a_2`$ in the multichannel calculations, which unambiguously determine the nature of spin-1 87Rb to be ferromagnetic. The global difference lies between $`0.3`$ and $`2.7`$ $`a.u.`$ over the uncertainty range. The difference is an increasing function of $`N_b`$, while $`a_0`$ and $`a_2`$ themselves are decreasing functions of $`N_b`$. For a given $`N_b`$, the range of possible values of $`a_0`$ and $`a_2`$ varies only weakly with $`C_6`$, as was the case for 23Na Crubellier:99 .

The results of $`a_0a_4`$ and $`a_2a_4`$ with uncertainties for 87Rb spin-2 are shown in Fig. 2. For all four values of $`N_b`$ and all three values of $`C_6`$ the uncertainty region is within the “polar” region, making the nature of 87Rb spin-2 condensate unambiguously determined. The pair ($`a_0a_4`$, $`a_2a_4`$) moves closer to the boundary between the polar and the cyclic regions as $`N_b`$ is increased but never reaches the boundary, within the present uncertainties. For a fixed value of $`N_b`$ $`a_2a_4`$ is increasing as a function of $`C_6`$, while $`a_0a_4`$ is almost independent of $`C_6`$. The uncertainty region for a fixed value of $`N_b`$ is very narrow (especially for higher $`N_b`$). The long axis of this region corresponds to the difference $`a_sa_t`$, whereas the narrow axis corresponds to the sum $`a_s+a_t`$.

The results for 85Rb spin-2 are shown on Fig. 3. In contrast to the case of 87Rb, here $`a_2a_4`$ and $`a_0a_4`$ are more dependent on the value of $`C_6`$ than on $`N_b`$, but only very little. $`a_2a_4`$ and $`a_0a_4`$ are slowly increasing functions of $`C_6`$ as well as of $`N_b`$. The uncertainties for all $`N_b`$ and values of $`C_6`$ unambiguously determine the nature of 85Rb spin-2 to be polar. The uncertainty region is again very narrow. The long axis of this region corresponds to $`a_t`$, whereas the narrow axis corresponds to $`a_s`$.

Since the graphs for spin-2 85Rb and 87Rb show only scattering lengths differences rather than scattering lengths, we summarize $`a_0`$, $`a_2`$ and $`a_4`$ in the following table. The scattering lengths for 85Rb show only very little dependence of $`N_b`$. Over the entire range of $`a_s`$, $`a_t`$, $`C_6`$ and $`N_b`$ the estimated scattering lengths for 85Rb are (in $`a.u.`$):

$`a_0=680\pm 80`$, $`a_2=540\pm 50`$, $`a_4=390\pm 30`$

The scattering lengths for spin-2 87Rb over the entire range of $`a_s`$, $`a_t`$ and $`C_6`$ are estimated as (in $`a.u.`$):

$`\begin{array}{ccccc}N_b& a_0& a_2& a_4& \\ 39& 90.5\pm 1& 97.5\pm 1& 106.8\pm 1& \\ 40& 89.0\pm 1& 94.8\pm 1& 103.6\pm 1& \\ 41& 87.7\pm 1& 92.4\pm 1& 100.5\pm 1& \\ 42& 86.4\pm 1& 90.2\pm 1& 97.4\pm 1& \end{array}`$

These numbers conservatively give the global uncertainties for each $`N_b`$. In the context of spinor condensates it is necessary to consider the actual allowed regions of the parameters, as we have done above and which permit us to draw meaningful conclusions.

To see how the results change when the multichannel and single-channel energy changes, we calculated the scattering lengths at various energies (with 1 pico Kelvin in the single and multichannel as the lowest value) to cover the relevant temperature for some experiments. This did not change our conclusions about the nature of the spinor Bose-Einstein condensates in rubidium.

Since the shape of the inner-wall potential is not known exactly and since we change it to have potentials with the four different values of $`N_b`$ we also performed the calculations with a quadratic inner-wall correction ($`c(rr_{min})^2`$ for $`r<r_{min}`$) instead of the $`\mathrm{arctan}`$ form. This did not change the conclusions and only changed the calculated scattering lengths by about $`0.1`$ %.

The present values of $`a_0`$, $`a_2`$ for spin-1 87Rb and $`a_0`$, $`a_2`$, $`a_4`$ for spin-2 85Rb and 87Rb are consistent with values obtained from Ref. Roberts:98 . Note that to carry out the calculations based on Ref. Roberts:98 , the correlations among $`a_s`$, $`a_t`$ and $`C_6`$ must be taken into account. We have separately calculated the values of $`a_0`$, $`a_2`$ for spin-1 87Rb and $`a_0`$, $`a_2`$, $`a_4`$ for spin-2 85Rb and 87Rb form $`a_s`$, $`a_t`$ for 85Rb and $`C_6`$ as given in Ref. Roberts:98 , and find that they support our classifications of the spinor condensates as presented in this paper. The new values from Ref. Roberts:01  allow us to determine a smaller uncertainty on the calculated scattering lengths, but they do not change our conclusions.

In summary our analysis based on the new results for the values of $`C_6`$, $`a_s`$, $`a_t`$ and the number of bound states in the triplet potential demonstrate that the nature of the ground states of 85Rb and 87Rb spin-2 condensates should be polar. In addition, the ground state of the 87Rb spin-1 condensate should be ferromagnetic. Therefore, in view of the known scattering parameters for 23Na, both ferromagnetic and polar spin-1 condensates are experimental accessible whereas no cyclic or ferromagnetic spin-2 condensate appears to exist for the most common rubidium isotopes.

We gratefully acknowledge D. M. Stamper-Kurn for suggesting that we perform a refined uncertainty analysis for spin-1 87Rb. we also thank E. Snyder for permitting the use of his computer programs, J. P. Burke for assistance and valuable discussions and J. L. Roberts for communication the results of Ref. Roberts:01  prior to publication. This work is supported by the National Science Foundation. N.N.K. acknowledges support from the Danish Fulbright Commission.


text:
# Acknowledgments

## Acknowledgments

This work was supported in part by the National Natural Science Foundation of China.

## Figure captions

The variations of the effective masses for the nucleon and the $`\mathrm{\Lambda }`$ hyperon, $`\delta M_i^{}=M_i^{}M_i`$ ($`i=N,\mathrm{\Lambda }`$), as functions of the quark mass correction, $`\delta m_q`$. The results in the QMF model with $`\chi _c=\frac{1}{2}kr^2`$ are shown by solid curves, while those with $`\chi _c=\frac{1}{2}kr^2(1+\gamma ^0)/2`$ are shown by dashed curves.

$`\mathrm{\Lambda }`$ Single particle energies in $`{}_{\mathrm{\Lambda }}{}^{41}\mathrm{Ca}`$, $`{}_{\mathrm{\Lambda }}{}^{91}\mathrm{Zr}`$, and $`{}_{\mathrm{\Lambda }}{}^{209}\mathrm{Pb}`$. $`\mathrm{QMF}^\mathrm{I}`$ and $`\mathrm{QMF}^{\mathrm{II}}`$ denote the models with $`\chi _c=\frac{1}{2}kr^2(1+\gamma ^0)/2`$ and $`\chi _c=\frac{1}{2}kr^2`$, respectively. The results in the QMC model  are also shown for comparison. The experimental data are taken from Refs. .

The scalar and vector potentials, $`U_S`$ and $`U_V`$, for the $`1s_{1/2}`$ $`\mathrm{\Lambda }`$ state in $`{}_{\mathrm{\Lambda }}{}^{41}\mathrm{Ca}`$ and $`{}_{\mathrm{\Lambda }}{}^{209}\mathrm{Pb}`$.


text:
# Noise parametric identification and whitening for LIGO 40-meter interferometer data

## I Introduction

The building of large interferometers is going to reach the phase of data taking: TAMA (Japanese)  is already working; GEO (British/German)  will begin to take data next year; LIGO (U.S.A.)   in 2002; VIRGO (French/Italian)  in 2003 (see e.g.    for general and exhaustive description of interferometric Gravitational Waves detectors).

The large amount of data produced by gravitational wave detectors will be essentially noise and, hopefully, buried in noise there will be the signal we are looking for. Ground-based interferometric detectors are sensible to a broad band frequencies (2-3Hz to more than 1kHz) in revealing relative displacement of test masses at the near and far extremities of interferometer arms due to GW signal, but unfortunately a lot of other factors can cause a similar displacement. The test masses are suspended to pendular structures in order to isolate them from seismic noise , but thermal noise of the suspension chain will cause a displacement of the mass . Also shot noise and radiation pressure of the laser will move the mirrors . The physicists are working in modeling all possible causes of noise in the interferometer giving out a sensitivity curve of the apparatus .

This curve is limited at low frequencies by seismic noise; in the middle band by thermal noise and at high frequencies (higher than $`0.71`$kHz) by shot-noise. The sensitivity curve of these detector is a broad-band noise plus several very narrow peaks due to the violin modes of the suspensions wires, that will make the detection of a gravitational signal in this frequency band very difficult. For this reason efforts have been made in the preparation of the analysis of data for cutting  out these resonances.

The analysis of data to detect the gravitational signal requires an accurate knowledge of the noise, which means a statistical characterization of the stochastic process, and in the case of local stationarity and Gaussian nature an accurate estimation of the Power Spectral Density (PSD).

Moreover if the output of the interferometer will be non stationary over a long period of time, we must be able in following the changes in the PSD. A way to achieve this is to estimate the PSD on a chunk of data at different interval of time, using classical techniques .

We proposed the use of adaptive methods to follow on line the change in the feature of the spectrum in such a way to have at any desired instant the correct curve for the PSD.

If we are able in identifying the noise of our detector we can also apply the whitening procedure of the data.

The goal of a whitening procedure is to make the sequence of data delta-correlated, removing all the correlation of the noise. Most of the theory of detection is in the frame of a wide-sense stationary Gaussian white noise, but in our problem the noise is surely a colored one and, in principle, there could be present non stationary and non Gaussian features. If we whiten the data, supposing hence to be in the frame of a stationary and Gaussian noise, we can apply the optimal algorithm detection.

Moreover when we are searching a transient signal of unknown form it is very important to have a whitened noise . The importance of whitening data is also linked to the possibility of reducing their dynamical range .

In this paper we show how to identify and how to whiten the noise data produced by LIGO-40 meter interferometric detector before applying any algorithm detection.

The data taken from Caltech 40-meter prototype interferometer in November 1994  have been given from the LIGO collaboration at disposal of the data analysis groups of other experiments to perform some test of data analysis algorithms. A typical PSD of the output of LIGO-40 meter data is plotted in figure 1. It is characterized by a huge number of spectral lines. Part of them are due to the violin modes, part to the laser noise and mostly of them to the 60Hz and its harmonics power supply.

We use these data to check the algorithms of parametric noise identification and the algorithms which perform the whitening on realistic data taken by an interferometric gravitational waves detector, having already analyzed and tuned these techniques on Virgo-like simulated data (see reference  ).

The concepts and the notation in this paper are described in detail in our earlier work . Sections II and III give a brief summary of the methods, but the reader is directed to Ref.  for details.

## II Parametric Modeling and Whitening

The advantages of power spectral parametric estimation with respect to the classical spectral methods are described in an exhaustive way in reference . We want to underline that one of its advantage in linked to the better spectral resolution we obtain if we suppose that the process we are analyzing is governed by a dynamical law, because we can use the estimation of autocorrelation function until a certain lag and then extrapolate its value at successive lags, under the dynamical hypothesis we made.

An other advantage is that we may compress the information of the PSD in a restricted number of parameters (the parameters of our model) and not in the full autocorrelation function. This can help us, for example, if we want to create a data base of noise sources.

In this paper we analyze only the Auto-Regressive (AR) parametric modeling because, as a byproduct, it offers the possibility to write down a linear whitening filter.

### A AR models

An Auto-Regressive process $`x[n]`$ of order $`P`$ ($`AR(P)`$ ) is governed by the relation

$$x[n]=\underset{k=1}{\overset{p}{}}a[k]x[nk]+w[n],$$
(1)
which links the value of the data $`x[n]`$ to the previous $`P`$ ones by the AR parameters $`a[k]`$, being $`w[n]`$ a Gaussian white noise process. The theoretical form of the PSD for such a process is given by

$$P_{AR}(f)=\frac{\sigma ^2}{|1+_{k=1}^Pa_k\mathrm{exp}(i2\pi kf)|^2},$$
(2)
being $`\sigma `$ one of the parameters to be estimated linked to the RMS of the process.

The relationship between the parameters $`a[k]`$ of the AR model and the autocorrelation function $`r_{xx}(n)`$ of the process $`x[n]`$ is given by the Yule–Walker equations

$$r_{xx}[k]=\{\begin{array}{cc}_{l=1}^pa_lr_{xx}[kl]\hfill & \text{for}k1\hfill \\ _{l=1}^pa_lr_{xx}[l]+\sigma ^2\hfill & \text{for}k=0.\hfill \end{array}$$
(3)
The problem of determining the AR parameters is the same of that of finding the optimal “weights vector” $`𝐰=w_k`$, for $`k=1,\mathrm{}P`$ for the problem of linear prediction . In the linear prediction we would predict the sample $`x[n]`$ using the $`P`$ previous observed data $`𝐱[n]=\{x[n1],x[n2]\mathrm{}x[nP]\}`$ building the estimate as a transversal filter:

$$\widehat{x}[n]=\underset{k=1}{\overset{P}{}}w_kx[nk].$$
(4)
We choose the coefficients of the linear predictor by minimizing a cost function that is the mean squares error $`ϵ=[e[n]^2]`$, being

$$e[n]=x[n]\widehat{x}[n]$$
(5)
the error we make in this prediction, obtaining the so called Normal or Wiener-Hopf equations

$$ϵ_{min}=r_{xx}[0]\underset{k=1}{\overset{P}{}}w_kr_{xx}[k],$$
(6)
which are identical to the Yule–Walker ones with

$`w_k`$ $`=`$ $`a_k`$ (7)
$`ϵ_{min}`$ $`=`$ $`\sigma ^2`$ (8)
This relationship between AR model and linear prediction assures us to obtain a filter which is stable and causal : this is the fundamental relation in building our whitening filter.

A method of solving the Yule–Walker equation is the Durbin algorithm .

The strategy of this method is the knowledge of the optimal $`(P1)`$th order filter to calculate from it the optimal $`P`$th order.

### B Link between AR model and whitening filter

When we find the AR $`P`$ parameters that fit a PSD of a noise process, what we are doing is to find the optimal vector of weights that let us reproduce the process at the time $`n`$ knowing the process at the $`P`$ previous time. All the methods that involves this estimation try to make the error signal (see equation (5) ) a white process in such a way to throw out all the correlation between the data (which we use for the estimation of the parameters). So one of the output of the filter which identifies the noise will be the whitened part of the data.

The Durbin algorithm, used to estimated the parameters $`a[k]`$, introduces in a natural way the Lattice structure for the whitening filter. The equations which represent our lattice filter in the time domain could be written

$`e_p^f[n]`$ $`=`$ $`e_{p1}^f[n]+k_pe_{p1}^b[n1],`$ (9)
$`e_p^b[n]`$ $`=`$ $`e_{p1}^b[n1]+k_pe_{p1}^f[n],`$ (10)
being $`e_p^b`$ the backward error, that is the error we make, in a backward way, in the prediction of the data $`x[np+1]`$ using $`p1`$ successive data, and $`e_p^f`$ the forward error, that is the error we made in estimating the data $`x[n]`$ using the $`P`$ previous data. The coefficients $`k_p`$ are the so-called reflection or PARCOR coefficients, being linked to the correlation that the data $`x[n]`$ has with the $`x[np]`$, neglecting the $`p1`$ values of $`x`$ in between. In figure 2 is showed how the lattice structure is used to estimate the forward and backward errors (see reference ).

Using a lattice structure we can implement the whitening filter following these steps:

* estimate the values of the autocorrelation function $`\widehat{r}_{xx}[k],0kP`$ of our process $`x[n]`$;
* use the Durbin algorithm to find the reflection coefficients $`k_p,1pP`$;
* implementation of the lattice filter with these coefficients $`k_p`$ initiating the filter $`e_0^f[n]=e_0^b[n]=x[n]`$.

In this way the forward error at the stage $`P`$-th is equivalent to the forward error of a transversal filter and represents the output of the whitening filter.

The procedure of whitening will be accomplished before applying the algorithms for the detection of gravitational signal of different wave forms. The level of whiteness of the data needed for the various algorithms could be different. We want now to introduce a parameter that let us quantify the level of whiteness of data at the output of whitening filter.

### C The “whiteness” of data: measure of flatness of PSD

The spectral flatness measure for a PSD is defined as

$$\xi =\frac{\mathrm{exp}(\frac{1}{N_s}_{N_s/2}^{N_s/2}\mathrm{ln}(P(f))𝑑f)}{\frac{1}{N_s}_{N_s/2}^{N_s/2}P(f)𝑑f}$$
(11)
where the integral is extended in the bandwidth of Nyquist frequency; this parameter satisfies

$$0\xi 1$$
(12)
If $`P(f)`$ is very peaky, then $`\xi 0`$, if $`P(f)`$ is flat than $`\xi =1`$.

Using the definition (11) the flatness for a process at the output of a whitening filter built with a minimum phase filter (as the AR filter is) is

$$\xi _e=\xi \frac{r_{xx}[0]}{r_{ee}[0]},$$
(13)
where $`r_{xx}[0]`$ and $`r_{ee}[0]`$ are the values of the autocorrelation function of the process before and after the whitening procedure, and $`\xi `$ is the value of flatness for the initial sequence .

### D Order Selection

The idea of the whitening filter is that the process we analyze is an autoregressive one and that once we have the AR parameters we can use them in the whitening filter.

In general we don’t know the order of our process, even if we suppose that it is an AR one. If it is an AR of order P, and we use an order $`p<P`$, the fitted spectrum will be smoother than the original one; if we choose an order $`p>P`$, there may be spurious peaks in the spectrum. In both cases the whitening will be not good.

If our process is not AR, the number of parameters could be in principle infinite. We must then fix a criterion that let us select the right order of the process, or at least the best one.

We used the classical order selection criteria , that is the Akaike information criterion (AIC), the forward prediction error (FPE) the Parzen’s criterion (CAT) and the minimum description length (MDL) one

$`AIC(P)`$ $`=`$ $`N\mathrm{log}ϵ(P)+2P,`$ (14)
$`FPE(P)`$ $`=`$ $`ϵ(P){\displaystyle \frac{N+P+1}{NP1}},`$ (16)
$`CAT(P)`$ $`=`$ $`\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{j=1}{\overset{P}{}}}{\displaystyle \frac{Nj}{Nϵ_j}}\right){\displaystyle \frac{NP}{Nϵ_P}},`$ (18)
$`MDL(P)`$ $`=`$ $`N\mathrm{log}ϵ(P)+P\mathrm{log}N,`$ (19)
where $`ϵ(P)`$ is the mean square error at the order $`P`$ and $`N`$ is the length of data. In literature the MDL criterion is considered the best among them, because it is robust with respect to the length of the sequence, while the others depend a lot on N.

## III Adaptive methods:Least Squares Lattice

The implementation of an adaptive filter follows two steps: the filtering of the input data and the adjustment of the filter parameters with which we process the data to the next iteration. The filters parameters are updated by minimizing a cost function. The way in which we build this cost function distinguishes the adaptive methods  The Least Squares based methods build their cost function using all the information contained in the error function at each step, writing it as the sum of the error at each step up to the iteration $`n`$ :

$$ϵ[n]=\underset{i=1}{\overset{n}{}}\lambda ^{ni}e^2(i|n),$$
(20)
being

$$e(i|n)=d[i]\underset{k=1}{\overset{N}{}}x_{ik}w_k[n],$$
(21)
where $`d`$ is the signal to be estimated, $`x`$ are the data of the process and $`w`$ the weights of the filter. We introduced the forgetting factor $`\lambda `$ that let us tune the learning rate of the algorithm. This coefficient can help when there are non stationary data in the sequence and we want that the algorithm have a short memory. If we have stationary data we fix $`\lambda =1`$.

There are two ways to implement the Least Squares methods for the spectral estimation: in a recursive way (Recursive Least Squares or Kalman Filters) or in a Lattice Filters using fast techniques . The first kind of algorithm, examined in , has a computational cost proportional to the square of the order of filter, while the cost of the second one is linear in the order $`P.`$ We will report only the application of Least Squares Lattice (LSL) algorithm on LIGO 40-meters data, because in our previous work the authors showed that this is the best choice among the adaptive algorithms the authors analyzed for the problem we have to face.

## IV Results on 40-meter LIGO data

We use two frames of the LIGO 40-meter interferometer data: the 18nov.2.frame and the 19nov.2.frame. We use the GRASP interface to read the data using only the frames taken with the interferometer locked. We are not interested to the absolute sensitivity of such interferometer, but only in the features of the Power Spectral Density, so we don’t use calibrated data, and we analyzed the IFO output of the data in ADC counts.This is a fast channel which was sampled at about $`10`$kHz.

In figure 1 a typical power spectral density of these data is displayed: it is characterized by a high number of lines, which could be mainly harmonics of the electromagnetic interference at 60Hz, by two broad peaks around $`180`$Hz and $`300`$Hz which could be due to the laser power supply, by the violin resonances around $`500`$-$`600`$Hz, and a continuous background, due to electronic noise, shot noise, seismic noise etc.

It is evident that we would need a high number of parameter to fit the spectrum and to build the respective whitening filter.

### A 18nov94.2.frame

We try to obtain an estimate of the order we need for whitening the noise PSD of this frame. The MDL criterion gives an order of 515. As first test we use this value to whiten the 18nov94.2.frame data.

In figure 3 we plotted the PSD of 18nov94.2.frame averaged on $`30`$ sequence of data of $`6`$s in input to the Durbin and LSL whitening filters. In the same figure we can see the AR(515) fit to the PSD and the PSD at the output of the whitening filters.

In the first figure it is evident that the filters succeed in whitening the overall behavior of the noise power spectrum, but not in eliminating all the spectral lines.

Besides the several harmonics spectral lines, the main problem in the whitening of these data is linked to the peaks in the range of frequencies from $`0`$ to $`1000`$Hz. These peaks are broad, and it is difficult to fit them with an all-poles models as it is evident in the zoom we made in this band of frequencies in figure 3.

Notwithstanding these problems, the values of flatness at the output of whitening filters is quite good as we see in figure 4 and in table I.

This is reflected in the time-domain output of the filter where it is evident that the rms of the data is remarkably reduced (see figure 5 ).

It is to be noted that the results of the adaptive algorithm reported in figure 4 and in table I are a little better that those of the ’static’ Durbin algorithm. These can be explained by the fact that the estimation of parameter of Durbin whitening needs a previous estimation of autocorrelation function, while the adaptive LSL filter estimate the parameters directly by the time data, avoiding the problem of errors propagation. Moreover in those runs the adaptive algorithm was tuned in such a way to have a long memory to let the filter continue to adjust its parameters with the input of new data.

#### 1 Whitening filter of higher order

Until now we used the order of whitening filter as the one estimate on the averaged PSD or autocorrelation function. Now we run the order selection criteria on a single sequence of data long $`6`$ s. The results are reported in table II.

We try the whitening procedure using firstly a filter of $`P=662`$ coefficients and then a filter of P=$`2682.`$

In figure 6 we report the PSD at the output of whitening filter of order $`P=662`$ as averaged on $`30`$ samples of $`6`$s of data. As it is evident there exists always the problem in the low frequency part of the spectrum due to the two broad peaks. A typical value of the flatness on one sample of data is $``$$`0.5`$ for the Durbin filter and $``$$`0.5`$$`5`$ for the LSL one.

In figure 7 we plotted the outputs of whitening filters with $`P=2682`$ as suggested from the AIC, FPE and CAT criteria. The higher order filters succeeded in eliminating the narrow spectral lines at higher frequencies, but the problems in the band under $`1000`$Hz still remain. The flatness of the output PSD is remarkably better in this case. A typical value for the Durbin whitening filter is $`\xi 0.6`$, while for the LSL one is $`\xi 0.7`$

#### 2 Whitening the data on a restricted band of frequencies.

By considering the fact that the signal detection algorithm will act on restricted band of frequencies, where it is most probable to find a gravitational signals and that the whitening procedure will be limited to an interval of frequencies and not on all the Nyquist interval, we try to whiten the 18nov94.2.frame on the band of frequencies from $`0`$ to $`1234`$Hz, where the algorithms have met with major difficulties in whitening the noise spectrum. To decimate the set of data at our disposal we let the data pass through a Butterworth low pass filter of 4 order that cuts off all the frequencies over $`1234`$Hz and then we re-sample the data with a sampling factor $`4.`$

The output of this pre-filter procedure has been sent to the whitening filter. We estimate the best order for the whitening filter using a single realization of $`6`$s of data and applying to it the order selection criteria (see table III).

We find again that the MDL gives the minimum value for $`P`$ . We whiten the data both with the value found by the MDL criterium, that is $`P=329`$, and with the value $`P=1154`$ which is the value found by the AIC and FPE criteria.

The PSD at the outputs of whitening filters with $`P=329`$ are plotted in figure 8.

Even if we have restricted the band of interest the difficulties of whitening the broad peaks remain. Moreover the band of frequencies around $`600`$Hz is characterized by several spectral lines plus some broad peaks, due to violin resonances . Most of the parameters of the whitening fit try to cancel out the broad peaks, so the whitening is worse also on narrow lines. Perhaps it could be useful separate the whitening in two steps, eliminate the narrow lines and whiten the broad peaks or augmenting the number of parameters. In fact if we use $`P=1154`$ the whitening we obtain in this band of frequencies is good as you can see in figure 9.

The typical value of the flatness for the whitening with $`P=1154`$ is between $`0.6`$ and $`0.8`$ and the Durbin filter and LSL filter have similar performance.

Moreover the choice of order of the filter is not tied up by the order selection criteria.

As we have already outlined  this number is an indicative one; the selection criteria depend also on the length of sample to analyze. If we use for example frame of data long $`26`$s, we find $`P=779`$ for the MDL and $`P=2694`$ for the FPE, AIC and CAT.

Just as a trial we make a whitening with $`P=3000`$ on a single set of these data. The result is plotted in figure 10. The whitening is very good: in fact the parameter of flatness is $`0.77`$ for the Durbin whitening algorithm and $`0.83`$ for the LSL whitening algorithm.

This is a proof that we can succeed in whitening the data if we choose a sufficiently high number of parameters to compensate the fact that we have all-poles models and that we need a high number of poles to whiten broad peaks.

### B 19nov94.2.frame

We perform the same kind of analysis even on the data set 19nov94.2.frame.

We find the order of the filter using the results of the order selection criteria on a set of data long $`6`$s, choosing $`P=`$ 663 that is the best order find by the MDL criterion (see table IV).

In figure 11 we showed the results of whitening filters on 19nov94.2.frame data at an order $`P=663`$.

The kind of problem we met in the whitening results is the same as in 18nov94.2.frame, that is besides the several lines at high frequencies, it is difficult to whiten the low frequency broad peaks. A typical value of the flatness at the output of Durbin filter is $`0.53`$, whilst for the LSL filter is $`0.55`$.

In figure 12 we report the whitening results obtained with filter of $`P=2351`$ order. The whitening is remarkably better .The values of the flatness are $``$$`0.6`$ for the Durbin filter and $``$$`0.63`$ for the LSL one.

#### 1 Whitening on a restricted band of frequencies

We perform the whitening procedure also on this frame using a restricted band of frequencies. The selection order criteria give in this case $`P=412`$ for the MDL and $`P=912`$ for the others (see table V).

The whitening filter with these order have been implemented on $`30`$ samples of length $`6`$s. The results are plotted in figures 13 14.

The whitening with $`P=912`$ is good; in this case the typical value for the flatness parameter is $`0.65`$ for Durbin filter and $`0.7`$ for LSL one. We make the whitening, just as a trial, on this set of data using filter of $`P=4096.`$ In figure 15 we plotted the results. The whitening procedure is good and the values of flatness is $`0.66`$ for Durbin filter and $`0.77`$ for LSL one. These values are not very different from the ones obtained with filters of order $`P=912`$.

## V Conclusion

We showed that we can obtain a whitened PSD on realistic interferometric data. In the two frames 18nov94.2 and 19nov94.2 of LIGO 40-meter interferometer we processed we found that, restricting our analysis to a band of frequencies and using a sufficiently high number of parameters, we can obtain whitened data in time domain on which we can apply our detection algorithms. It could be useful to statistically characterize the LIGO data to check out the presence of non stationary or non linearities in order to study the behavior of the whitening filters. Indeed we verified that the adaptive algorithm has a better performance with respect to the static algorithm so it could be useful if we have to face with non stationary data.

In whitening the LIGO 40-meter data we applied the filters directly to the data but, in principle, it could be useful before applying the whitening filters to cancel out the harmonics of $`60`$Hz  . We want to underline that we apply these techniques to check their capability in whitening the power spectral density, but we could apply them in a more efficient way if we know in advance the kind of analysis we want to perform after the whitening, because we can choose to have a smoother level of whiteness or to analyze a narrow band of frequencies. The filters we use belong to the family of Kalman filter, but the dynamical law we suppose to underline the process is a general one, so in principle we can whiten all the spectral lines of the PSD. This could be advantageous if we want to perform the whitening in one step, because we know that the ’risk’ of finding a periodic gravitational signal with an amplitude comparable to the ones of the violin modes or other spectral lines in the PSD is very low.

## Acknowledgments

The authors wish to thank the LIGO Laboratory for having made its prototype interferometer data from November 1994 available to the Virgo collaboration.They especially wish to acknowledge the LIGO authors in reference  who had a significant role in the acquisition of the data. Moreover the authors wish to thanks Albert Lazzarini for his courtesy and Andrea Vicere’ for having made at their disposal these data and for very useful discussions.


text:
# Chemical freeze-out parameters at RHIC from microscopic model calculations

## Abstract

The relaxation of hot nuclear matter to an equilibrated state in the central zone of heavy-ion collisions at energies from AGS to RHIC is studied within the microscopic UrQMD model. It is found that the system reaches the (quasi)equilibrium stage for the period of 10-15 fm/$`c`$. Within this time the matter in the cell expands nearly isentropically with the entropy to baryon ratio $`S/A=150170`$. Thermodynamic characteristics of the system at AGS and at SPS energies at the endpoints of this stage are very close to the parameters of chemical and thermal freeze-out extracted from the thermal fit to experimental data. Predictions are made for the full RHIC energy $`\sqrt{s}=200`$ AGeV. The formation of a resonance-rich state at RHIC energies is discussed.

Thermalization and chemical equilibration of hot and dense nuclear matter produced in ultrarelativistic heavy-ion collisions is a topic of great importance for the interpretation of current SPS and RHIC results. In our investigation we study the relaxation process in central Au+Au collisions at $`\sqrt{s}=200`$ AGeV within the microscopic transport UrQMD model . Earlier studies at AGS and SPS energies revealed that the central reaction volume defined by a cubic cell of volume $`V=125`$ fm3 is well suited for this kind of study . It contains enough particles to be treated as a statistical system, and its macroscopic characteristics become isotropic after some time. Figure 2 depicts the velocity distributions of hadrons in the cell in transverse ($`x`$ and $`y`$) and in longitudinal ($`z`$) directions. At $`t=3`$ fm/$`c`$ the longitudinal velocity distribution differs considerably from the distributions in the transverse plane, while at $`t=5`$ fm/$`c`$ the magnitudes and widths of all three distributions become very close to each other. Isotropy of the velocity distributions results in the isotropy of pressure in the cell. Pressure in longitudinal and in transverse direction is shown in Fig. 2(a) for AGS, SPS, and RHIC energies. It is widely believed that the thermalization at RHIC sets in quite early. Indeed, at RHIC the pressure in the cell becomes isotropic at $`t5`$ fm/$`c`$ compared with $`t8`$ fm/$`c`$ (SPS) and $`t10`$ fm/$`c`$ (AGS). Starting from $`t=5`$ fm/$`c`$ the results of the microscopic calculations are compared with the predictions of the statistical model (SM) of an ideal hadron gas . The values of the energy density $`\epsilon `$, baryon density $`\rho _B`$, and strangeness density $`\rho _S`$, determined microscopically, are used as an input to obtain particle yields, partial energy densities, pressure, and entropy density via the temperature $`T`$, baryochemical potential $`\mu _B`$, and strangeness chemical potential $`\mu _S`$. As seen in Fig. 2(a), the microscopic pressure is very close to the grand canonical pressure after the onset of the (quasi)equilibrium stage.

The entropy density per baryon $`S/A`$, defined both microscopically and macroscopically (via the Gibbs thermodynamic identity $`Ts=\epsilon +P\mu _B\rho _B\mu _S\rho _S`$) is presented in Fig. 2(b). The expansion in the cell proceeds nearly isentropically with $`S/A=12`$ (AGS), 32 (SPS), and 150 (RHIC) (cf. $`(S/A)^{\mathrm{therm}.\mathrm{fit}}14`$ (AGS) and $`36`$ (SPS) ). Therefore, it would be very interesting to compare the UrQMD estimate $`(S/A)^{\mathrm{RHIC}}=150170`$ with the value extracted from the thermal model fit to RHIC experimental data.

The equation of state (EOS) in the $`(P,\epsilon )`$-plane is shown in Fig. 4. For all three energies it can be well approximated by a simple linear dependence $`P/\epsilon =0.12`$ (AGS), and 0.15 (SPS and RHIC). Note, that this version of the model does not imply the formation of the quark-gluon plasma, therefore, there are no kinks in this plot that can be attributed to quark-hadron phase transition. The evolution of the EOS in the $`(T,\mu _B)`$-plane is depicted in Fig. 4 together with the chemical freeze-out and the thermal freeze-out parameters obtained from a thermal fit to experimental data at AGS and SPS energies. One can see that temperatures and chemical potentials in the cells at the beginning and at the end of the equilibrated stage are close to the thermodynamic parameters of chemical and thermal freeze-out, respectively. The UrQMD predicts that at chemical freeze-out $`T_{\mathrm{chem}.\mathrm{FO}}^{\mathrm{RHIC}}=195\pm 5`$ MeV and $`\mu _{\mathrm{chem}.\mathrm{FO}}^{\mathrm{RHIC}}=43\pm 2`$ MeV, while at thermal freeze-out $`T_{\mathrm{therm}.\mathrm{FO}}^{\mathrm{RHIC}}=130\pm 5`$ MeV and $`\mu _{\mathrm{therm}.\mathrm{FO}}^{\mathrm{RHIC}}=43\pm 3`$ MeV, i.e., the evolution of nuclear matter in the cell proceeds at constant baryon chemical potential. Calculations show very weak (within few MeV) difference between the results for Au+Au at $`\sqrt{s}=130`$ AGeV and at $`\sqrt{s}=200`$ AGeV. A thermal fit to particle ratios measured by the STAR Collaboration (preliminary data) yields $`T_{\mathrm{chem}.\mathrm{FO}}=190\pm 20`$ MeV and $`\mu _q=15\pm 5`$ MeV , i.e. $`\mu _B=45\pm 15`$ MeV, which is in remarkable agreement with the UrQMD calculations.

It is interesting to check the correspondence of the cell conditions at chemical freeze-out to the criterion $`E/N_hM_h+3/2T1`$ GeV introduced in . The mean energy

This means that not only the temperature, but also the mean mass of a particle is increased in the cell at RHIC energies, i.e., there should be more heavy resonances compared with the SPS cell. Therefore, the ratios of hadronic abundances are studied (see Fig. 5). Here the results are presented separately for non-strange and strange baryons and mesons. In the baryon sector the resonances dominate over the strange and non-strange baryons until the end of the simulations. This can be taken as an indication of the creation of long-lived resonance-rich matter. The fraction of baryon resonances is almost 70% of all baryons in the cell at RHIC at $`5t19`$ fm/$`c`$, while at SPS and AGS the number of baryon resonances decreases from 70% to 35%, and from 60% to 25%, respectively. The meson fractions of resonances shrink within the time interval $`5t19`$ fm/$`c`$ from 60% to 30% (RHIC), 50% to 20% (SPS), and 40% to 15% (AGS). But at RHIC energies the

The results of our study can be summarised as follows. The formation of long lived resonance-abundant matter is found. UrQMD predicts that $`T=195\pm 5`$ MeV, $`\mu _B=43\pm 2`$ MeV, $`150S/A170`$, and $`E/N_h1`$ GeV at chemical freeze-out in central Au+Au collisions at RHIC. The equation of state has a linear dependence $`P=0.15\epsilon `$. The UrQMD cell calculations show that strangeness to entropy ratio monotonically increases with rising $`\sqrt{s}`$ as $`R_S0.025`$ (AGS), 0.04 (SPS), and 0.05 (RHIC).


text:
# Black Hole Excision for Dynamic Black Holes

\[

## Abstract

We extend previous work on 3D black hole excision to the case of distorted black holes, with a variety of dynamic gauge conditions that respond naturally to the spacetime dynamics. We show that in evolutions of highly distorted, rotating black holes, the combination of excision and the gauge conditions we use is able to drive the coordinates to a frame in which the system looks almost static at late times. Further, we show for the first time that one can extract accurate waveforms from these simulations, with the full machinery of excision and dynamic gauge conditions. The evolutions can be carried out for long times, far exceeding the longevity and accuracy of better resolved 2D codes.

\]

Introduction. The long term simulation of black hole (BH) systems is one of the most challenging and important problems in numerical relativity. For BHs, the difficulties of accuracy and stability in solving Einstein’s equations numerically are exacerbated by the special problems posed by spacetimes containing singularities. At a singularity, geometric quantities become infinite and cannot be handled easily by a computer.

Traditionally, in the 3+1 approach the freedom in choosing the slicing has been used to slow down the approach of the time slices towards the singularity (“singularity avoidance”), while allowing them to proceed outside the BH . Singularity avoiding slicings are able to provide accurate evolutions, allowing one to study BH collisions and extract waveforms , but only for limited evolution times. Combining short full numerical evolutions with perturbation methods, one can even study the plunge from the last stable orbit of two BHs . But a breakthrough is required to push numerical simulations far enough to study orbiting BHs, requiring accurate evolutions exceeding time scales of $`t100M`$. In 3D, traditional approaches have not been able to reach such time scales, even in the case of Schwarzschild BHs. Characteristic evolution codes, on the other hand, are well-adapted to the long-term evolution of single black holes , but here we concentrate on non-characteristic methods for their ease of generalization to multiple black holes.

A more promising approach involves cutting away the singularity from the calculation (“singularity excision”), assuming it is hidden inside an apparent horizon (AH) . This work has been progressing, from early spherical proof of principle  to recent 3D developments . However, beyond a few spherical test cases , excision has yet to be used in conjunction with live gauge conditions designed to respond to both the dynamics of the BH and the coordinate motion through the spacetime.

In this paper we extend recent excision work  to the case of distorted, dynamic BHs in 3D, using a new class of gauge conditions. These gauge conditions not only respond naturally to the true spacetime dynamics, but also drive the coordinates to a frame where the system looks almost static at late times. We show that not only are the evolutions accurate as indicated by the mass associated with the apparent horizon, but also that very accurate waveforms can be extracted, even when the waves carry only a tiny fraction of the energy of the spacetime. We also show that the 3D evolutions of dynamic BHs we are now able to perform, are superior, in terms of accuracy, stability, and longevity, to previous 3+1 BH simulations, whether carried out in full 3D or even when restricted to 2D. These results indicate that BH excision can be made to work under rather general circumstances, and can significantly improve the length of the evolutions, and the accuracy of the waveforms extracted, which will be crucial for gravitational wave astronomy.

Initial Data. For this paper we consider single distorted BH spacetimes  that have been used to model the late stages of BH coalescence . Following , the initial 3-metric $`\gamma _{ab}`$ is chosen to be

$$ds^2=\mathrm{\Psi }^4\left[e^{2q}\left(d\eta ^2+d\theta ^2\right)+\mathrm{sin}^2\theta d\varphi ^2\right],$$
(1)
where the “Brill wave” function $`q`$ is a general function of the spatial coordinates, subject to certain regularity and fall off restrictions, that can be tailored to produce very distorted 3D BHs interacting with nonlinear waves. The radial coordinate $`\eta `$ is logarithmic in the cartesian radius $`r`$. There are two classes of data sets used here corresponding to even- and odd-parity distortions. The even-parity data have vanishing extrinsic curvature, while the cases containing an odd-parity component have nontrivial extrinsic curvature $`K_{ij}`$. As shown in , these distorted BH data sets can include rotation as well, corresponding to spinning, distorted BHs that mimic the early merger of two orbiting BHs. Hence they make an ideal test case for the development of our techniques. We leave the details of the construction of these BH initial data sets to Refs. . An important point that we wish to emphasize is that such data are not of the Kerr-Schild form with ingoing coordinates at the horizon. That form of initial data sets has been recently advocated since it is not conformally flat  and is well adapted to inward propagation of quantities at the horizon.

Evolution and Excision Procedures. Our simulations have been performed using what we refer to as the “BSSN” version of the 3+1 evolution equations , which we have found to have superior stability properties when compared to standard formulations. As detailed in , we actively force the trace of the conformal-traceless extrinsic curvature $`\stackrel{~}{A}_{ij}`$ to remain zero, and we use the independently evolved “conformal connection functions” $`\stackrel{~}{\mathrm{\Gamma }}^i`$ only in terms where derivatives of these functions appear. All the simulations described below have been performed using a 3-step iterative Crank-Nicholson scheme and a radiative (Sommerfeld) outer boundary condition. We refer the reader to Ref.  for the details of the numerical implementation.

We use the simple excision approach described in . Our algorithm is based on the following ideas: (a) Excise a cube contained inside the AH that is well adapted to cartesian coordinates; (b) Use a simple boundary condition at the sides of the excised cube: copying of time derivatives from their values one grid point out along the normal directions; (c) Use centered (non-causal) differences in all terms except for advection terms on the shift (terms of the form $`\beta ^i_i`$). For these terms we use second order upwind along the shift direction. These simplifications in excision reduce the complexity in the algorithm, avoid delicate interpolation issues near the excision boundary, and have allowed us to make rapid progress. Currently, the method is implemented for non-moving excision regions, although they are allowed to grow. One can hope that colliding black holes can be treated even with this restriction through the use of co-moving coordinates. A more detailed description of our excision algorithm can be found in Ref. .

Gauge Conditions. For the lapse we use a hyperbolic slicing condition motivated by the Bona-Massó family of slicing conditions ,

$$_t\alpha =\alpha ^2f\left(\alpha \right)(KK_0),$$
(2)
where $`K`$ is the trace of the extrinsic curvature, $`K_0`$ is its value in the initial data, and $`f`$ is a (positive) function of $`\alpha `$ which we specify below. With this condition, the lapse will evolve as long as $`\alpha ^2f(\alpha )`$ and $`KK_0`$ are non-vanishing.

For the shift $`\beta ^i`$ we have considered families of elliptic, parabolic, and hyperbolic conditions that relate the shift with the evolution of the conformal connection functions $`\stackrel{~}{\mathrm{\Gamma }}^i`$. We obtain parabolic and hyperbolic shift conditions by making either $`_t\beta ^i`$ or $`_t^2\beta ^i`$ proportional to the elliptic operator for $`\beta ^i`$ contained in the “Gamma freezing” condition $`_t\stackrel{~}{\mathrm{\Gamma }}^k=0`$ (see Ref. ), itself closely related to the well known minimal distortion family . Elliptic conditions have the disadvantage of requiring boundary data at the excision region where it is difficult to know what to impose, while parabolic conditions force a strong restriction on the stability of the differencing scheme: $`\mathrm{\Delta }t(\mathrm{\Delta }x)^2`$ (this is true for explicit schemes, implicit schemes have no such restriction). We have then concentrated on hyperbolic conditions of the form

$$_t^2\beta ^i=\zeta _t\stackrel{~}{\mathrm{\Gamma }}^i\xi _t\beta ^i,$$
(3)
where $`\zeta `$ and $`\xi `$ are positive functions. We call such evolution conditions for the shift hyperbolic “driver” conditions (see ).

In the spirit of the puncture method for evolutions , we use a BSSN scheme with the usual time-dependent conformal factor $`e^{4\varphi }`$ and an additional time-independent conformal factor $`\mathrm{\Psi }^4`$ that comes from the initial data. In our examples we use $`\zeta =k/\mathrm{\Psi }^4`$, where $`k`$ is a positive constant. The division by $`\mathrm{\Psi }^4`$ helps to slow down the evolution of the shift in the vicinity of the black hole. We have found it important to add a dissipation term with a constant coefficient $`\xi `$ in order to reduce some initial oscillations in the shift. Notice that in contrast with $`k`$, the coefficient $`\xi `$ is not dimensionless (it has dimensions of inverse length), so in practice we rescale it using the total mass of the system. Experience has shown that by tuning the value of $`\xi `$ we can almost freeze the evolution of the system at late times.

The parameters used for all simulations described below are: $`\alpha `$ is given by Eq. (2), with $`f=2/\alpha `$; $`\beta ^i`$ is given by Eq. (3) with $`\zeta =0.75/\mathrm{\Psi }^4,\xi =3/M`$ ($`M`$ is the initial ADM mass of the system). As initial conditions we take $`\alpha =1`$, $`\beta ^i=0`$, $`_t\alpha =_t\beta ^i=0`$, except in one case mentioned below where we perform a single maximal slicing solve to obtain a more appropriate initial lapse. Given these initial conditions, we let the gauge conditions take care of the rest. We use the same gauge parameters for all results in this paper, whether they are applied to Schwarzschild, distorted, or rotating BHs, showing the strength and generic nature of these conditions.

Results. The first example we show is Schwarzschild, written in the standard isotropic coordinates used in many BH evolutions. Note that with this initial data and our starting gauge conditions, the BH should evolve rapidly. If $`\alpha `$ and $`\beta ^i`$ were held fixed at their initial values, the slices would hit the singularity at $`t=\pi M`$. Instead, $`\alpha `$ and $`\beta ^i`$ work together with excision to rapidly drive the coordinates to a frame where the system looks essentially static, corresponding to the true physical situation.

In Fig. 1 we show the radial metric function $`g_{rr}/\mathrm{\Psi }^4`$ vs. time. The grid covers an octant with $`128^3`$ points ($`\mathrm{\Delta }x=0.2`$, $`M=2`$). Appropriate symmetry conditions are applied on the faces of the octant for the different dynamical variables. We have checked that removing the octant symmetry (while using a lower resolution) does not change the results for the evolution times reported here (in particular no instabilities were encountered, cmp. ). Notice that the metric begins to grow, as it does without a shift, but as the shift builds up the growth slows down significantly. At this stage, the system is effectively static, even though we started in the highly dynamic isotropic coordinates. We also show the time development of $`\alpha `$ and $`\beta ^r`$ in Fig. 2, which evolve rapidly at first but then effectively freeze, bringing the whole system to an almost static state by $`t=10M`$. The evolution of the metric and gauge variables then proceeds only very slowly with time until the simulation is stopped well after $`t=200M`$. The decision to stop the code at this time is simply a CPU time consideration, but we notice that in this and all following examples the code is stopped once all the interesting dynamics have finished.

Figure 3 shows the AH mass $`M_{AH}=\sqrt{\mathrm{Area}_{AH}/16\pi }`$, determined with a 3D AH finder . For comparison we also show the value of $`M_{AH}`$ for the 3D run without shift, and for a highly resolved 2D simulation with no shift and no excision . The 2D code uses maximal slicing, so the coordinate time $`t`$ refers to different slices, but the slicings turn out to be even more similar than is to be expected from Eq. (2). While the 3D simulation with shift and excision continues well beyond $`t=200M`$, the 2D result becomes inaccurate and the code crashes due to axis instabilities by $`t=150M`$, and the 3D run without shift crashes already by $`t=50M`$. Notice that in the 2D case, after around $`t=35M`$, $`M_{AH}`$ grows rapidly due to numerical errors associated with grid stretching. With excision and our new gauge conditions, the 3D run has less than a few percent error by $`t=200M`$, while the 2D case has more than 100% error before it crashes at $`t150M`$. For the excision run, notice also that while there is some initial evolution in the metric and the coordinate size of the AH (see Figs. 1 and 2) the AH mass changes very little.

Next, we turn to a truly dynamic, even-parity distorted BH. This system contains a strong gravitational wave that distorts the BH, causing it to evolve, first nonlinearly, and then oscillating at its quasi-normal frequency, finally settling down to a static Schwarzschild BH. This provides a test case for our techniques with dynamic, evolving BH spacetimes, and allows us to test our ability to extract gravitational waves with excision for the first time. In this case, in the language of , we choose the Brill wave parameters to be $`Q_0=0.5`$, $`\eta _0=0`$, $`\sigma =1`$, corresponding to a highly distorted BH with $`M=1.83`$. Just as before, we use a grid that covers one octant, with $`128^3`$ points and $`\mathrm{\Delta }x=0.2`$

In Fig. 4 we show the AH mass $`M_{AH}`$ as a function of time for the distorted BH simulations carried out in both 2D and 3D. $`M_{AH}`$ grows initially as a nonlinear burst of gravitational waves is absorbed by the BH, but then levels off as the BH goes into a ring-down phase towards Schwarzschild. Fig. 5 shows the proper polar circumference of the AH divided by its proper equatorial circumference. This ratio allows an estimate of the size of the local dynamics during the run. Notice how the horizon starts far from spherical (with a ratio close to 2), it later oscillates from prolate to oblate and back again, and finally settles on a sphere (with a ratio of 1).

In the 3D case, the gauge conditions and excision quickly drive the metric to an almost static configuration, as the system itself settles towards a static Schwarzschild BH. The evolution is terminated at around $`t=160M`$. To our knowledge, distorted BHs of this type have never been evolved for so long, nor with such accuracy, in either 2D or 3D. By comparison, in the more highly resolved 2D case with zero shift and no excision, the familiar grid stretching effects allowed by the gauge choice lead to highly inaccurate evolutions after some time with the error in $`M_{AH}`$ again approaching 100% when the code finally crashes at $`t100M`$.

In Fig. 6, we show the results of extracting waves from the evolution of this highly distorted, excised BH. Using the standard gauge-invariant waveform extraction technique, the Zerilli function is shown for both the 2D and 3D simulations discussed above. There is a slight but physically irrelevant phase difference in the two results due to differences in the slicing; otherwise the results are remarkably similar (the waves are extracted at the same Schwarzschild radius in both cases). This shows conclusively that the excision and live gauge conditions do not adversely affect the waveforms, even if they carry a small amount of energy (around $`10^3M_{ADM}`$ in this case).

We now turn to a rather different type of distorted BH, including rotation and general even and odd-parity distortions. In the language of Ref. , the parameters for this simulation are $`Q_0=0.5`$, $`\eta _0=0`$, $`\sigma =1`$, $`J=35`$, corresponding to a rotating distorted BH with $`M=7.54`$ and an effective rotation parameter $`a/M=0.62`$. Previously, such data sets could be evolved only to about $`40M`$ . For the purposes of this paper we have chosen an axisymmetric case so that we can compare the results to those obtained with a 2D code. Since this example is much more demanding, we have found it important in order to increase the accuracy of our runs to perform a single initial maximal solve to reduce the initial gauge dynamics. The symmetries of this example are now not consistent with the evolution of just one octant. However, we still have reflection symmetry on the $`z=0`$ plane, so we evolve only the positive $`z`$ half of the domain. The grid used in this case has $`199^2\times 100`$ points and $`\mathrm{\Delta }x=0.4`$. The gauge conditions work well even in the presence of rotation: the shift drives the evolution to an almost static state as the system itself settles down to a Kerr BH. The metric functions (not shown) evolve in a similar way to those shown before, essentially freezing at late times. In Fig. 7, we show the extracted waveforms, now computed using the imaginary part of the Newman-Penrose quantity $`\mathrm{\Psi }_4`$ (e.g. ), which includes contributions from all $`\mathrm{}`$modes at the same time. The results from the 2D and 3D codes agree very closely, except for a slight phase shift due to slicing differences, until the 2D code becomes inaccurate and crashes. The 3D simulation continues well beyond this point, and is terminated at $`t=120M`$.

Conclusions. We have extended recently developed 3D BH excision techniques, using a new class of live gauge conditions that dynamically drive the coordinates to a frame where the metric looks essentially static at late times, when the system itself settles to a stationary Kerr BH. Our techniques have been tested on highly distorted, rotating BHs, and are shown to be very robust. For the first time, excision is tested with wave extraction, and waveforms are presented and verified. The results are shown to be more accurate, and much longer lived, than previous 3D simulations and even better resolved 2D simulations of the same data. Such improvements in BH excision are badly needed for more astrophysically realistic BH collision simulations, which are in progress and will be reported elsewhere.

Acknowledgements. This work was supported by AEI. Calculations were performed using the Cactus code at AEI, NCSA, PSC, and RZG.