text:
# Bounds on the number of time steps for simulating arbitrary interaction graphs

## 1 Introduction

Simulating Hamiltonian evolutions of arbitrary quantum systems on a quantum computer is an idea that goes back already to Feynman . This would be an attractive application of future quantum computers since there are no known efficient classical algorithms for simulating generic dynamics of many particle systems. Here a quantum computer can be any quantum system provided that its time evolution can be controlled in a universal way. In particular, the problem of controllability includes the question which Hamiltonian evolutions can be simulated efficiently by the considered Hamiltonian system \[2–13\] Assume the natural Hamiltonian $`H`$ to act on an $`n`$-fold tensor product Hilbert space

$$_n:=\mathrm{}$$

where each $``$ denotes the Hilbert space of a qudit, i.e. a $`d`$-dimensional system. The Lie algebra $`su(d)`$ of traceless Hermitian operators on $``$ is a $`m:=(d^21)`$-dimensional real vector space. Let $`B:=\{\sigma _\alpha \alpha =1,\mathrm{},m\}`$ be an orthogonal basis of $`su(d)`$ with respect to the trace inner product $`A|B:=tr(A^{}B)/d`$ for $`A,Bsu(d)`$.

Then the most general qudit-qudit interaction on $`n`$ coupled qudits is given by

$$H_J:=\underset{k<l}{}\underset{\alpha \beta }{}J_{kl;\alpha \beta }\sigma _\alpha ^{(k)}\sigma _\beta ^{(l)},$$

where $`J`$ is chosen to be a real symmetric $`mn\times mn`$-matrix with zeros for $`k=l`$. Note that the symmetry of the coupling matrix $`J`$ does not imply any physical symmetry of the interaction. It is a consequence of our redundant representation that turns out to be very useful. The coupling matrix $`J`$ consists of $`m\times m`$-blocks. The $`m\times m`$-matrix $`J_{kl}`$ given by the block at position $`(k,l)`$ describes the coupling between the qudits $`k`$ and $`l`$. We have $`J_{lk}=J_{kl}^T`$, i. e. the matrix describing the coupling between the qudits $`l`$ and $`k`$ is just the transpose of the matrix describing the coupling between $`k`$ and $`l`$. The blocks on the diagonal are zero matrices.

In the setting discussed here and in most other articles on simulation of Hamiltonians the only possibilities of external control are given by local unitaries on each qudit. We assume that one is able to implement them independently. Formally, all control operations are elements of the group $`𝒞:=SU(d)SU(d)\mathrm{}SU(d)`$. A common approximation is to assume that all operations in $`𝒞`$ can be implemented arbitrarily fast (“fast control limit”). The simulation of Hamiltonians is based on the following “average Hamiltonian” approach that has successfully been used for describing Nuclear Magnetic Resonance experiments since many years (e.g. ).

Let $`t_1,t_2,\mathrm{},t_n`$ be positive numbers and $`u_1,u_2,\mathrm{},u_N𝒞`$ be control operations. Then the algorithm

> perform $`u_1`$, wait $`t_1`$, perform $`u_1^{}`$, perform $`u_2`$, wait $`t_2`$, perform $`u_2^{}`$, …, perform $`u_N`$, wait $`t_N`$, perform $`u_N^{}`$

implements approximatively the unitary evolution

$$\underset{j=1}{\overset{N}{}}\mathrm{exp}(u_j^{}Hu_jt_j).$$

If the times $`t_j`$ are small compared to the time scale of the natural evolution according to the natural Hamiltonian $`H`$ this is approximatively the evolution according to the average Hamiltonian

$$\underset{j}{}t_ju_jHu_j^{}/\tau ,$$
(1)
where $`\tau :=_jt_j`$ is the slow down factor of teh evolution, i.e., the time overhead of the simulation. For investigating the time overhead and the number of time steps of a simulation it turns out to be useful to work with the coupling matrices instead of considering the Hamiltonians themselves. To express the effect of the control operations on the coupling matrix $`J`$ note that any unitary operation $`uSU(d)`$ corresponds to a rotation on the $`m`$-dimensional sphere via the relation

$$u^{}\left(\underset{\alpha }{}c_\alpha \sigma _\alpha \right)u=\underset{\alpha }{}\stackrel{~}{c}_\alpha \sigma _\alpha ,$$

where the vector $`\stackrel{~}{c}=(\stackrel{~}{c}_1,\stackrel{~}{c}_2,\mathrm{},\stackrel{~}{c}_m)`$ is obtained by applying a rotation $`USO(m)`$ on the vector $`c=(c_1,c_2,\mathrm{},c_m)`$. It is straightforward to verify that conjugation of $`H_J`$ by $`v:=u^{(1)}u^{(2)}\mathrm{}u^{(n)}`$ corresponds to conjugation of $`J`$ by a block diagonal matrix of the form

$$V:=U^{(1)}U^{(2)}\mathrm{}U^{(n)}\underset{k=1}{\overset{n}{}}SO(m).$$

The condition for correct simulation is hence given by

$$\stackrel{~}{J}=\underset{j}{}t_jV_jJV_j^T,$$
(2)
where the orthogonal matrix $`V_j`$ corresponds to the unitary $`v_j`$ for $`j=1,\mathrm{},N`$.

The question of optimal simulation has been completely solved so far only for the case of two-qubit Hamiltonians . Optimal simulation protocols are constructed in . The number of time steps is between $`1`$ and $`3`$.

For $`n`$ qudits we have shown in  that the eigenvalues of the coupling matrix $`J`$ and $`\stackrel{~}{J}`$ provide lower bounds on the simulation time overhead. Graph theoretical notions can also provide upper bounds. In certain cases the bounds are known to be tight .

Here we do not focus on the time overhead but on the number of time steps in the $`n`$ partite case. Upper bounds for decoupling (switching off the Hamiltonian) and time-reversal (simulating of $`H`$ by $`H`$) are directly given by the parameters of known schemes . Lower bounds for time-reversal have been derived in . These results shown that time overhead and the number of time steps are not connected in any obvious way. Here we address the general problem of mutual simulation of pair-interactions on $`n`$ qudits.

In Section 2 we derive lower bounds on the number of time steps for the general problem of mutual simulation of pair-interactions. The type of coupling is assumed to be the same between all nodes, but the strengths and the signs may vary. Like the bounds on the time overhead , the bound derived here make us of the spectrum of the corresponding coupling matrices. In the special case that one wants to cancel some interactions while keeping others (“the simulation of a certain interaction graph”) the bounds refer directly to the spectra of the corresponding adjacency matrices. Surprisingly, the greatest lower bound can be given if the smallest eigenvalue of the adjacency matrix is irrational. In Section 3 the bounds shall be applied to three different cases where long range interactions are present between all spins. The first two cases are the simulation of nearest neighborhood coupling in a square lattice and a cyclic chain. The third case is a certain interaction graph that has been proposed to prepare the states of a graph code .

## 2 Lower bounds on the number of time steps

We restrict our attention to interactions with an additional symmetry, namely Hamiltonians of the following form

$$H:=\underset{k<l}{}w_{kl}\underset{\alpha \beta }{}c_{\alpha \beta }\sigma _\alpha ^{(k)}\sigma _\beta ^{(l)}.$$
(3)
The matrix $`W:=(w_{kl})`$ is a real symmetric $`n\times n`$-matrix with zeros on the diagonal. It describes the coupling strengths and the signs of the interactions between all qudits. The matrix $`C=(c_{\alpha \beta })`$ is a real symmetric $`m\times m`$-matrix characterizing the type of the coupling. This means that all qudits interact with each other via the same interaction and that only the coupling strengths and the signs vary. It is important that in this special case the coupling matrix $`J`$ can be expressed as a tensor product of $`W`$ and $`C`$, i. e., $`J=WC`$.

To derive a general lower bound for simulating arbitrary interactions $`\stackrel{~}{J}`$ by a tensor product interaction $`J=WC`$ it is useful to observe that the decisive condition in eq. (2) is invariant with respect to the following rescaling of interactions: Multiply each $`m\times m`$-block $`k,l`$ of $`J`$ and $`\stackrel{~}{J}`$ with the same factor $`r_{kl}`$. In the case $`WC`$ we can therefore assume w.l.o.g. that $`W`$ is a matrix with only $`1`$ as non-diagonal entries as long as we restrict our attention to Hamiltonians with complete interaction graphs, i.e., $`w_{kl}0`$ for $`kl`$. Formally, rescaling is denoted as follows. Let $`A/B`$ be the entry-wise quotient of the matrices $`A`$ and $`B`$ provided that $`A`$ has only zero entries at those positions where $`A`$ has also a zero. Then we consider the simulation of $`\stackrel{~}{J}/(WI)`$ by $`KC`$ where $`K`$ is the matrix with only one as non-diagonal entries and $`0`$ in the diagonal and $`I`$ is the matrix with only $`1`$ as entries.

In order to emphasize that the assumption $`w_{kl}=1`$ is not an assumption on the real physical coupling strength we formulate the following theorem for general $`W`$ and use rescaling only in the proof.

###### Theorem 1 (Lower bound)

Let $`J:=WC`$ be the coupling matrix of the system Hamiltonian, $`\stackrel{~}{J}`$ an arbitrary coupling matrix of the interaction that is simulated, and $`\mu `$ the time overhead. Denote the minimal and maximal eigenvalues of $`C`$ by $`\lambda _{\mathrm{min}}`$ and $`\lambda _{\mathrm{max}}`$, respectively and its rank by $`r(C)`$. Let $`I`$ be the $`m\times m`$-matrix whose all entries are $`1`$. Let $`s`$ be the number of eigenvalues of $`\stackrel{~}{J}/(WI)`$ that are not contained in the interval

$$:=[\mu \lambda _{\mathrm{max}},\mu \lambda _{\mathrm{min}}].$$

Then the number of time steps required to simulate $`H_{\stackrel{~}{J}}`$ by $`H_J`$ is at least $`s/r(C)`$.

Proof: The condition for a scheme $`t_1,V_1,t_2,V_2,\mathrm{},t_N,V_N`$ to be a simulation of $`\stackrel{~}{J}`$ reads

$$\underset{j=1}{\overset{N}{}}t_jV_j(WC)V_j^T=\stackrel{~}{J}.$$
(4)
Since the matrices $`V_j`$ are block-diagonal we can rescale each $`m\times m`$-block such that we obtain

$$\underset{j=1}{\overset{N}{}}t_jV_j(KC)V_j^T=\stackrel{~}{J}/(WI).$$

In the following we denote the rescaled coupling matrix $`\stackrel{~}{J}/(WI)`$ by $`J^{}`$.

Set $`R:=_{j=1}^Nt_jV_j(\mathrm{𝟏}C)V_j^T`$. By adding $`R`$ on both sides we obtain

$$\underset{j=1}{\overset{N}{}}t_jV_j\left((K+\mathrm{𝟏})C\right)V_j^T=J^{}+R.$$
(5)
The rank of the matrix $`(K+\mathrm{𝟏})`$ is $`1`$ since all its entries are $`1`$. Consequently, the rank of the left hand side of eq. (5) is at most $`Nr(C)`$.

Set $`a:=\mu \lambda _{\mathrm{min}}`$ and $`b:=\mu \lambda _{\mathrm{max}}`$. Note that the eigenvalues of $`R`$ are contained in the interval $`[a,b]`$. The rank of the right hand side is at least the number of eigenvalues of $`J^{}`$ outside the interval $`=[b,a]`$. This is seen as follows.

Let $`P`$ and $`Q`$ be the projections onto the sums of all eigenspaces of $`J^{}`$ with eigenvalues smaller than $`b`$ and greater than $`a`$, respectively. Clearly, $`s`$ is equal to the dimension of the image of $`PQ`$. Denote by $`\widehat{J}`$ and $`\widehat{R}`$ the $`s\times s`$-submatrices of $`J^{}`$ and $`R`$ defined by $`(PQ)J^{}(PQ)`$ and $`(PQ)R(PQ)`$, respectively.

Due to the choice of $`P`$ and $`Q`$ the spectrum of $`\widehat{J}`$ is contained in the interval $`(\mathrm{},b)(a,\mathrm{})`$. The spectrum of $`\stackrel{~}{R}`$ is contained in the interval $`[a,b]`$ since the minimal (maximal) eigenvalues of a matrix cannot decrease (increase) when projecting the matrix.

We prove the theorem by showing that $`\widehat{J}+\widehat{R}`$ has full rank. Set $`f:=\frac{1}{2}(a+b)`$. From the triangle inequality it follows that for every unit vector $`|\mathrm{\Psi }^s`$ one has

$`(\widehat{J}+\widehat{R})|\mathrm{\Psi }`$ $`=`$ $`(\widehat{J}+f\mathrm{𝟏}+\widehat{R}f\mathrm{𝟏})|\mathrm{\Psi }`$
$``$ $`(\widehat{J}+f\mathrm{𝟏})|\mathrm{\Psi }(\widehat{R}f\mathrm{𝟏})|\mathrm{\Psi }.`$
The eigenvalues of the shifted operators $`\widehat{J}+f\mathrm{𝟏}`$ and $`\widehat{R}f\mathrm{𝟏}`$ are contained in the intervals $`(\mathrm{},\frac{1}{2}(ba))(\frac{1}{2}(ba),\mathrm{})`$ and $`[\frac{1}{2}(ba),\frac{1}{2}(ba)]`$, respectively. This implies that

$$(\widehat{J}+f\mathrm{𝟏})|\mathrm{\Psi }(\widehat{R}f\mathrm{𝟏})|\mathrm{\Psi }>0$$

because the norms can be bounded by the absolute values of the eigenvalues: $`(\widehat{J}+f\mathrm{𝟏})|\mathrm{\Psi }>|\frac{1}{2}(ba)|`$ and $`(\widehat{R}f\mathrm{𝟏})|\mathrm{\Psi }|\frac{1}{2}(ba)|`$. This completes the proof. $`\mathrm{}`$

The following upper bound can be easily derived from Caratheodory’s theorem .

###### Theorem 2 (Upper bound)

Every simulation that is possible can be achieved within

$$\frac{n(n1)}{2}m^2+1$$

time steps.

Proof If $`\stackrel{~}{J}`$ can be simulated by $`J`$ with time overhead $`\tau `$ then $`\stackrel{~}{J}/\tau `$ is in the convex span of the matrices $`V_jJV_j^T`$ with notation as above. The dimension of this convex set is at most $`m^2n(n1)/2`$ since the diagonal blocks are empty and each matrix $`V_jJV_j^T`$ is symmetric. Caratheodory’s theorem states that each point in an $`M`$ dimensional convex set can be written as a convex sum of at most $`M+1`$ extreme points. $`\mathrm{}`$

There are interesting cases where the bound of Theorem 1 can be tightened. Assume $`\stackrel{~}{J}=\stackrel{~}{W}C`$ with the same matrix $`C`$ as the interaction that is used for the simulation. In other words, only the strengths and the signs of some interactions should be changed. By the same rescaling trick as above we can consider the problem to simulate $`(\stackrel{~}{W}/W)C`$ by $`KC`$ where $`K`$ has only $`1`$ as non-diagonal entries.

###### Theorem 3 (Lower bound)

Let $`WC`$ be the coupling matrix of the natural Hamiltonian and $`\stackrel{~}{W}C`$ the coupling that we want to simulate. Assume all non-diagonal entries of $`W`$ to be non-zero.

1. Let $`C`$ be a positive semidefinite matrix. Then the number of time steps is at least the number of positive eigenvalues of $`\stackrel{~}{W}/W`$.
2. Let $`C=diag(1,1,\mathrm{},1)`$ be the $`m\times m`$-identity matrix. Then the number of time steps is at least $`nk`$, where $`k`$ is the multiplicity of the smallest eigenvalue $`\mu _{\mathrm{min}}`$ of $`\stackrel{~}{W}/W`$.
3. Let the natural coupling be $`C:=diag(0,0,1)`$, i.e., we have $`\sigma _z\sigma _z`$ interactions between all spin-1/2-particles. Let the set of local control operations be restricted to $`i\sigma _x`$–transformations. Then one requires at least $`nk`$ time steps with $`k`$ as in Case 2. If $`\mu _{\mathrm{min}}`$ is irrational then at least $`n`$ steps are necessary. In any case, $`n(n1)/2+1`$ time steps are always sufficient.

Note that Case (3) is of special interest since it deals with simulation procedures that do not rely on any first order approximation. In this case all summands in eq. (1) commute and the unitary transformation implemented by the simulation scheme coincides exactly with the exponent of the average Hamiltonian.

Proof (of Theorem 3)

Case 1 This statement is a corollary of Theorem 1 since $`\stackrel{~}{J}/(WI)=(\stackrel{~}{W}/W)C`$. The number of positive eigenvalues of this tensor product matrix is $`r(C)`$ times the number of positive eigenvalues of $`\stackrel{~}{W}/W`$ if $`r(C)`$ is the rank of $`C`$.

Case 2 Consider the rescaled problem to simulate $`AC`$ by $`KC`$. In this case the right hand side of eq. (5) reduces to

$$AC+\tau \mathrm{𝟏},$$
(6)
where $`\tau :=_jt_j`$ is the time overhead and $`A:=\stackrel{~}{W}/W`$. In  we have shown that $`\tau `$ is at least $`\mu _{\mathrm{min}}`$. The reason is that the spectrum of $`\tau K`$ has to majorize the spectrum of $`A`$ since $`A`$ is a convex sum of conjugates of $`\tau K`$. Therefore the rank of the matrix in expression (6) is at least $`m(nk)`$. Since the left hand side of eq. (5) has at most the rank $`mN`$ the number of time steps is at least $`nk`$.

Case 3 In this case one can use a more convenient formalism. We drop the matrix $`C`$ and characterize the interactions by $`W`$ and $`\stackrel{~}{W}`$. If a $`zz`$-interaction is conjugated by $`i\sigma _x`$-transformations only two possibilities occur. The interaction term between spin $`k`$ and spin $`l`$ acquires a minus sign if exactly one spin of both is subjected to a conjugation. The interaction is unchanged if either no spin or both spins are conjugated. Instead of representing the time step $`j`$ by the $`3n\times 3n`$ block diagonal matrix $`V_j`$ we can represent it by an $`n\times n`$ diagonal matrix $`X_j`$. The diagonal entries are $`\pm 1`$ and indicate which spins are subjected to conjugation. Then eq. (4) reduces to

$$\underset{j=1}{\overset{N}{}}t_jX_jKX_j=A.$$

We add the identity matrix on both sides. Due to $`X_j\mathrm{𝟏}X_j=\mathrm{𝟏}`$ we obtain

$$\underset{j}{}t_jX_j(K+\mathrm{𝟏})X_j=A+\underset{j}{}t_j\mathrm{𝟏}.$$

The rank of the left hand side is at most the number $`N`$ of time steps. To estimate the rank of the right hand side note that the time overhead $`\tau :=_jt_j`$ is at least $`\mu _{\mathrm{min}}`$. Hence the dimension of the kernel of $`A+\tau \mathrm{𝟏}`$ can at most be the multiplicity of the eigenvalue $`\mu _{\mathrm{min}}`$. This shows that the number of time steps is at least $`nk`$.

Assume $`\mu _{\mathrm{min}}`$ to be irrational. Then the time overhead $`\tau `$ is necessarily greater than $`\mu _{\mathrm{min}}`$. This can be seen as follows. The optimization with respect to the time overhead reduces to the following convex problem. Consider the matrix $`X_jKX_j`$ for an arbitrary time step $`j`$. Its non-diagonal entries are $`\pm 1`$ and indicate which interactions acquire a minus sign in $`j`$th step. In graph-theoretical language, the set of matrices that can occur as $`X_jKX_j`$ are exactly the Seidel matrices of complete bipartite graphs (see ). Then the optimal $`\tau `$ is the minimal positive number such that $`A/\tau `$ is in the convex span of the set of Seidel matrices of complete bipartite graphs. Geometrically, the convex span is a polytope having the Seidel matrices as its extreme points. It is embedded in the $`n(n1)/2`$ dimensional vector space of real symmetric matrices with zeros on the diagonal. Let $`O`$ be the origin. Consider the semi-line $`\nu A`$ for $`\nu 0`$. Then the optimal simulation is the unique intersection point $`P`$ of the semi-line with the boundary of the polytope. The quotient of the distance between $`0`$ and $`A`$ and between $`0`$ and $`P`$ is the optimal time overhead. This quotient can never be irrational. The reason is that $`P`$ has rational entries since it is the solution of a linear equation over the field $``$ of rational numbers. Hence $`\tau `$ is greater than $`\mu _{\mathrm{min}}`$ and $`A+\tau \mathrm{𝟏}`$ has necessarily full rank. This proves that we need at least $`n`$ time steps.

To see that $`n(n1)/2+1`$ time steps are always sufficient we can argue as in the proof of Theorem 1 with Caratheodory’s theorem. The dimension of the convex span of the matrices $`X_jKX_j`$ is at most $`n(n1)/2`$. $`\mathrm{}`$

In the following section we will consider the task to cancel some interactions and keep others. Then $`A:=\stackrel{~}{W}/W`$ has only $`1`$ and $`0`$ as entries. In graph-theoretical language, it is the adjacency matrix of the desired interaction graph.

It is surprising that it is relevant for our lower bound whether the smallest eigenvalue of the adjacency matrix is irrational (see Case 3.). It is not clear whether this is only a feature of our proof or whether there is a true connection to the irrationality of graph spectra.

In the following section we will use graph spectra for deriving lower bounds based on the above theorems. We show that some of them are almost tight by sketching simulation schemes based on well-known results on selective decoupling.

## 3 Applications

In this section we consider the simulation of special interaction graphs that appeared in the literature in various applications. Some interesting models in quantum information theory refer to quite idealized types of Hamiltonians like nearest neighbor interactions. If the natural Hamiltonian contains long range interactions between all nodes one may try to simulate the idealized interaction. Then the problem is to cancel the unwanted terms without destroying the desired interactions. The examples below show that it may cause a large number of time steps to cancel unwanted long-range interactions no matter how fast they are decreasing with the distance. As long as they are not neglected, the control sequences that cancel them may be rather long. Here we restrict our attention to $`\sigma _z\sigma _z`$ interactions between $`n`$ qubits.

### 3.1 Cyclic and open spin chains

We examine first a quantum system of $`n`$ spins equally spaced on a one-dimensional lattice. The interaction between the spins may decrease with the distance between the spins. We consider the problem to simulate a chain with only nearest neighbor interactions, i.e. we have to cancel the interactions between all non-adjacent pairs. The desired interaction graphs can be seen in Fig. 1.

The adjacency matrix of the cyclic spin chain is $`S+S^T`$, where $`S`$ denotes the cyclic shift in $`n`$ dimensions. The eigenvalues of $`S`$ and $`S^T`$ are the $`n`$th roots of unity. The eigenvalues of $`S+S^T`$ are given by

$$2\mathrm{}\left(\mathrm{exp}(i2\pi l/n)\right)\text{ with }\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}l<n.$$

For odd $`n`$ the lower bound on the number of time steps is $`n`$ since the least eigenvalue is irrational for all $`n>3`$. If $`n`$ is even it is rational and has multiplicity $`1`$. Hence the bound is $`n1`$ in this case. There are numbers $`n`$ where this bound is almost tight. This is shown by the following example. Let $`n`$ be an even number with the property that a Hadamard matrix of dimension $`n/2`$ exists. This is for instance the case for each power of $`2`$ (see ). We construct a simulation scheme that consists of $`2`$ subroutines.

The first subroutine simulates the interaction between the pairs $`\{1,2\},\{3,4\},\mathrm{},`$ $`\{n1,n\}`$ and the second simulates $`\{n,1,\},\{2,3\},\mathrm{},\{n2,n1\}`$. Note that all the pairs in the same subroutine are disjoint. The problem to simulate the interactions between disjoint pairs is a special instance of well-known “cluster decoupling”  where the interaction between independent cliques are cancelled and the interactions within the same clique remains. It can be achieved using Hadamard matrices having at least the number of cliques as dimension (compare ). The entries $`\pm 1`$ in column $`j`$ determine which spins are conjugated in the step $`j`$. The dimension is the number of time steps of the decoupling subroutine.

Using this method, we need $`n/2`$ steps in each subroutine. Therefore we have given a simulation with $`n`$ steps. whereas our lower bound is $`n1`$. In general the number of steps for simulating the circle grows only linearly in $`n`$. This shows that the lower bound is quite good even for general $`n`$.

The adjacency matrix corresponding to the open spin chain has the eigenvalues

$$2\mathrm{cos}\left(\frac{\pi }{n+1}i\right),i=1,\mathrm{},n.$$
(7)
The smallest eigenvalue is irrational for all $`n>2`$. By Theorem 3 (Case 3) we conclude that the number of time steps is at least $`n`$.

### 3.2 Square lattice

We consider a quantum system of $`n=l^2`$ spins located on a two-dimensional square lattice. For simplicity assume $`l`$ to be even. We want to simulate a lattice with only nearest neighbor interactions.

The desired interaction graph is shown on the left of Fig. 2. This kind of interaction can for instance be used for preparing the initial state in the ‘One-Way Quantum Computer’ proposed in . The eigenvalues of the corresponding adjacency matrix $`A`$ are known in graph theory :

$$2\mathrm{cos}\left(\frac{\pi }{l+1}i\right)+2\mathrm{cos}\left(\frac{\pi }{l+1}j\right),i,j=1,\mathrm{},l.$$
(8)
We first consider the time overhead. An upper bound is given by $`4`$ since this is the chromatic index of the graph . It is easy to see that the minimal eigenvalue of $`A`$ is given by

$$\lambda _{\mathrm{min}}=2\mathrm{cos}\left(\frac{\pi }{l+1}l\right)+2\mathrm{cos}\left(\frac{\pi }{l+1}l\right).$$
(9)
By Theorem 3 (Case 3) the lower bound on the number of time steps is $`n`$ since the smallest eigenvalue is irrational. Note that this example shows that the complexity measures time overhead and number of time steps may differ significantly.

An upper bound on the number of time steps can be obtained as follows. The graph has $`2(l1)l`$ edges. We can partition the edges into $`4`$ sets of edges such that each set contains only disjoint interacting pairs. These $`4`$ partitions are shown in Fig. 2. The simulation consists of $`4`$ subroutines simulating one the interactions in one of the $`4`$ classes. For each subroutine we choose Hadamard matrices with a dimension that is at least the number of cliques. The numbers of cliques are $`l^2/2`$ or $`l^2/2+l`$ in each subroutine. Since there exist Hadamard matrices for every power of $`2`$ the square lattice graph can always be simulated in $`O(l^2)=O(n)`$ time steps.

### 3.3 Graph codes

The computational power of different $`n`$ spin interactions is not well understood yet. It would be interesting to know which $`n`$ qubit transformations can easily be implemented when a certain Hamiltonian is given. However, one of the few examples where the power of specific Hamiltonians is directly used (without using them to implement 2-qubit gates) is the preparation of states of graph codes proposed in . Here the codes states are obtained by the free time evolution according to a Hamiltonian with $`\sigma _z\sigma _z`$ interactions. The graph representing a code is the interaction graph that can be used for preparing the states. We assume the natural interaction to be equal $`zz`$ interaction between all $`6`$ spins. and would like to simulate the interaction graph in Fig 3.

The eigenvalues of the ‘wheel’ in Fig. 3 can easily be computed by any computer algebra system. They are

$$1+\sqrt{6},\frac{1}{2}\sqrt{5}\frac{1}{2},\frac{1}{2}\sqrt{5}\frac{1}{2},\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}\sqrt{6},\frac{1}{2}\sqrt{5}\frac{1}{2},\frac{1}{2}\sqrt{5}\frac{1}{2}.$$

The minimal eigenvalue is $`1/2\sqrt{5}/2`$. By Theorem 3 (Case 3) the minimal number of time steps is $`6`$. An implementation with $`12`$ time steps is given as follows. The scheme consists of $`3`$ subroutines each consisting of $`4`$ time steps. As above each subroutine simulates the interaction between disjoint cliques and cancels the interaction between different cliques. Subroutine $`1`$ has the cliques $`\{1,2,6\},\{4,5\},\{3\}`$. The clique partitions in subroutine $`2`$ and $`3`$ are $`\{3,4,6\},\{1,5\},\{2\}`$ and $`\{1\},\{4\},\{2,3\},\{5,6\}`$, respectively. In each subroutine decoupling the different cliques can be achieved by Hadamard matrices of dimension $`4`$ since no subroutine has more than $`4`$ cliques. Hence we have $`4`$ time steps in each subroutine.

## 4 Conclusions

We have derived lower bounds on the number of time steps for simulating arbitrary pair-interactions between $`n`$ qudits. Like the lower bounds on the time overhead, they make use of the spectrum of the coupling matrices. However, there is no direct connection between both complexity measures since the time overhead refers to spectral majorization while the bounds on the number of time steps refer to the number of eigenvalues not contained in a certain interval. We have shown an example where the number of time steps is of the order $`n`$ but the time overhead is independent of $`n`$.

## Acknowledgements

Thanks to Debbie Leung for helpful corrections. This work has been founded by the BMBF project “Informatische Prinzipien und Methoden bei der Steuerung komplexer Quantensysteme”.


format: text

text:
# Mirage Cosmology in M-theory

## I Introduction

Recently there has been renewed interest in the cosmological model based on the brane universe since this idea can be applied to string theory. The idea of brane universe is that our observed universe is a three-brane embedded in a higher dimensional space RTA ; RS . Many cosmological models regarding this have been studied. These models can largely be classified into two categories. One is that the brane is a static solution of the underlying theory and the cosmological evolution is due to the time evolution of the energy density on the brane BDL . The other is that the cosmological evolution of the brane universe is due to the motion of the brane in the background of the bulk as well as the matter density on the brane CRP ; Keki ; vsl ; bwilliams .

One interesting model among the second category is the so-called mirage cosmology presented by Kehagias and Kiritsis Keki . The idea is that the motion of the brane through the bulk, ignoring its back reaction to the ambient geometry, induces cosmological evolution on the brane even when there is no matter field on the brane. The crucial mechanism underlying the construction of this formalism is the coupling of the probe brane to the background gauge field. They derived Friedman-like equations for various bulk background field solutions within type II string theory.

This model was studied extensively by others. The mirage cosmology with non-trivial dilaton field was studied by the author kimdil . Since the dilaton as well as the induced metric affects the effective matter density, the cosmological evolution with nontrivial dilaton profile is different from the one without dilaton. The motion of a three-brane in the background of type 0B string theory was examined in Ref. type0 . Brane inflation for tachyonic and non-tachyonic type 0B string theories was studied and it is known that the presence of tachyon slows down the inflation in mirage cosmology. Brane cosmology in the background of D-brane with NS $`B`$ field was studied by Youm youmnsb . The corrections to the Friedman equations due to nonzero NS $`B`$ field were obtained and analyzed for various limits. The mirage cosmology for non-planar probe universe was studied in Ref. youmclo . There the author considered the spherical probe brane wrapped around the sphere part of various background spacetimes and commented its relevance to the giant graviton giantg . It is known that the mirage cosmology approach matches with the familiar junction condition approach when there is just one extra dimension steerparry .

Since type 0B string theory is defined on the world sheet of type IIB theory by performing a nonchiral Gliozzi-Sherk-Olive (GSO) projection gso , so far the study on mirage cosmology is mainly based on the type IIB string theory. In this paper we will extend the idea of mirage cosmology to M-theory. As a concrete example, we consider the M-theory background generated by a stack of non-threshold (M2,M5) bound states. We study the cosmological evolution by mirage matter in this background.

The organization of the paper is as follows. In Sec. II we briefly review the (M2,M5) background. In Sec. III we construct the action of a probe M5-brane under the background ignoring the back reaction. In Sec. IV, we consider the cosmological evolution of the brane. We estimate the range of $`r`$ where the formalism is valid. Effective matter density on the probe brane is expressed as a function of the scale factor. We also discuss a limiting case of the result to compare with the known result from the type IIB string background. Finally we conclude and discuss our results in Sec. V.

## II The brane background

The supergravity background we will consider is the one generated by a stack of parallel non-threshold (M2,M5) bound states malru . The metric for this eleven dimensional supergravity solution can be written as camram

$`ds^2`$ $`=`$ $`f^{1/3}h^{1/3}\left[(dx^0)^2+(dx^1)^2+(dx^2)^2+h\left\{(dx^3)^2+(dx^4)^2+(dx^5)^2\right\}\right]`$ (1)
$`+`$ $`f^{2/3}h^{1/3}\left[dr^2+r^2d\mathrm{\Omega }_4^2\right],`$
where $`d\mathrm{\Omega }_4^2`$ is the metric of a unit 4-sphere and $`f`$ and $`h`$ are given by

$$f=1+\frac{R^3}{r^3},h^1=\mathrm{sin}^2\phi f^1+\mathrm{cos}^2\phi .$$
(2)
The above solution appeared in Ref. rutse  and was interpreted as a two-brane lying within a five-brane. The M5-brane component extends along the directions $`x^0,\mathrm{},x^5`$, while the M2-brane lies along $`x^0,x^1,x^2`$. The angle $`\phi `$ in Eq. (2) carries the mixing of the M2- and M5-branes in the bound state. The radial parameter $`R`$ is defined as $`R^3\mathrm{cos}\phi \pi Nl_p^3`$, where $`l_p`$ is the eleven dimensional Planck length and $`N`$ is the number of the bound states of the stack. We also have a non-vanishing value of the four-form field strength $`F^{(4)}`$ given by

$`F^{(4)}`$ $`=`$ $`\mathrm{sin}\phi _r(f^1)dx^0dx^1dx^2dr3R^3\mathrm{cos}\phi ϵ_{(4)}`$ (3)
$``$ $`\mathrm{tan}\phi _r(hf^1)dx^3dx^4dx^5dr,`$
where $`ϵ_{(4)}`$ represents the volume form of the unit four-sphere $`\mathrm{S}^4`$. We parameterize the metric of $`\mathrm{S}^4`$ as

$$d\mathrm{\Omega }_4^2=\frac{1}{1\rho ^2}d\rho ^2+(1\rho ^2)d\varphi ^2+\rho ^2d\mathrm{\Omega }_2^2,$$
(4)
where $`d\mathrm{\Omega }_2^2`$ is the metric of a unit two-sphere $`\mathrm{S}^2`$ (which we will label $`\theta ^1`$ and $`\theta ^2`$). The ranges of $`\rho `$ and $`\varphi `$ are $`0\rho 1`$ and $`0\varphi 2\pi `$ respectively. Then, the three-form and six-form potential relevant for our calculation can be written as

$`C^{(3)}`$ $`=`$ $`\mathrm{sin}\phi f^1dx^0dx^1dx^2R^3\mathrm{cos}\phi \rho ^3d\varphi ϵ_{(2)}+\mathrm{tan}\phi hf^1dx^3dx^4dx^5,`$ (5)
$`C^{(6)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}\phi \mathrm{cos}\phi f^1R^3\rho ^3dx^0dx^1dx^2d\varphi ϵ_{(2)}`$ (6)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1h\mathrm{cos}^2\phi }{\mathrm{cos}\phi }}f^1dx^0dx^1dx^2dx^3dx^4dx^5`$
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}\phi R^3\rho ^3hf^1dx^3dx^4dx^5d\varphi ϵ_{(2)},`$
where $`ϵ_{(2)}`$ is the volume form of $`\mathrm{S}^2`$. For later use we can write the metric components of Eq. (1) as

$`g_{00}`$ $`=`$ $`g_{11}=g_{22}=f^{1/3}h^{1/3},g_{33}=g_{44}=g_{55}=f^{1/3}h^{2/3}g(r),`$
$`g_{rr}`$ $`=`$ $`f^{2/3}h^{1/3},g_{\rho \rho }=f^{2/3}h^{1/3}{\displaystyle \frac{r^2}{1\rho ^2}},g_{\varphi \varphi }=f^{2/3}h^{1/3}r^2(1\rho ^2),`$ (7)
$`g_{\theta ^1\theta ^1}`$ $`=`$ $`f^{2/3}h^{1/3}r^2\rho ^2,g_{\theta ^2\theta ^2}=f^{2/3}h^{1/3}r^2\rho ^2\mathrm{sin}^2\theta ^1.`$
## III The probe M5-brane action

We consider a probe M5-brane moving in the background of (M2, M5) bound states which shares $`(x^3,x^4,x^5)`$ directions with the background and wraps $`\mathrm{S}^2`$ ($`\theta ^1,\theta ^2`$). The dynamics of M5-brane, ignoring all fermions, is given by the so-called PST action pst . In PST formalism the world volume fields are a three-form field strength $`F`$ and a scalar $`a`$ (the PST scalar). The action consists of three terms

$$S=T_{M5}d^6\xi [_{DBI}+_{H\stackrel{~}{H}}+_{WZ}],$$
(8)
where $`T_{M5}`$ is the tension of the M5-brane $`T_{M5}=1/(2\pi )^5l_p^6`$. The explicit forms of $`_{DBI}`$, $`_{H\stackrel{~}{H}}`$ and $`_{WZ}`$ are given by

$`_{DBI}`$ $`=`$ $`\sqrt{\mathrm{det}(\gamma _{ij}+\stackrel{~}{H}_{ij})},`$ (9)
$`_{H\stackrel{~}{H}}`$ $`=`$ $`{\displaystyle \frac{1}{24(a)^2}}ϵ^{ijklmn}H_{lmn}H_{jkp}\gamma ^{pq}_ia_qa,`$ (10)
$`_{WZ}`$ $`=`$ $`{\displaystyle \frac{1}{6!}}ϵ^{ijklmn}\left\{P[C^{(6)}]_{ijklmn}+10H_{ijk}P[C^{(3)}]_{lmn}\right\},`$ (11)
where $`\gamma `$ is the induced metric on the M5-brane worldvolume

$$\gamma _{ij}(\xi )=g_{\mu \nu }(x)\frac{x^\mu }{\xi ^i}\frac{x^\nu }{\xi ^j},$$
(12)
and $`P[C^{(3)}]`$ and $`P[C^{(6)}]`$ are the pullbacks of the corresponding background potentials. The field $`H`$ and $`\stackrel{~}{H}`$ are defined as

$`H_{ijk}`$ $`=`$ $`F_{ijk}P[C^{(3)}]_{ijk},`$ (13)
$`\stackrel{~}{H}^{ij}`$ $`=`$ $`{\displaystyle \frac{1}{3!\sqrt{\mathrm{det}\gamma }}}{\displaystyle \frac{1}{\sqrt{(a)^2}}}ϵ^{ijklmn}_kaH_{lmn}.`$ (14)
To write down the action explicitly, we take the worldvolume coordinates $`\xi ^i(i=0,1,\mathrm{},5)`$ in the static gauge as

$$\xi ^i=(x^0,x^3,x^4,x^5,\theta ^1,\theta ^2).$$
(15)
In this system of coordinates the variables $`x^1,x^2,r,\rho ,\varphi `$ are functions of $`\xi ^i`$ in general. We assume that these variables depend only on time and there is a translational symmetry along the $`x^1`$ and $`x^2`$ directions. Then the configuration we are interested in is described by

$$r=r(t),\rho =\rho (t),\varphi =\varphi (t),$$
(16)
where $`t=x^0`$. The induced metric $`\gamma _{ij}`$ is calculated, in terms of eleven dimensional spacetime metric, as

$`\gamma _{00}`$ $`=`$ $`|g_{00}|+g_{rr}\dot{r}^2+g_{\rho \rho }\dot{\rho }^2+g_{\varphi \varphi }\dot{\varphi }^2,`$
$`\gamma _{33}`$ $`=`$ $`\gamma _{44}=\gamma _{55}=g(r),`$
$`\gamma _{\theta ^1\theta ^1}`$ $`=`$ $`g_{\theta ^1\theta ^1},\gamma _{\theta ^2\theta ^2}=g_{\theta ^2\theta ^2},`$
where the dot($``$) denotes the derivative with respect to $`t`$. We also assume that the only non-vanishing components of $`H`$ are those of $`P[C^{(3)}]`$, i.e. $`H_{x^3x^4x^5}H_{345}`$ and $`H_{x^0\theta ^1\theta ^2}H_0`$. By fixing the gauge, the auxiliary field $`a`$ can be eliminated from the action at the expense of losing the manifest covariance. Choosing the gauge $`a=x^0=t`$, the only nonzero component of $`\stackrel{~}{H}`$ is

$$\stackrel{~}{H}_{\theta _1\theta _2}=\sqrt{\frac{g_{\theta ^1\theta ^1}g_{\theta ^2\theta ^2}}{g^3}}H_{345}=f^{7/6}h^{4/3}r^2\rho ^2\sqrt{\widehat{g}^{(2)}}H_{345},$$
(17)
with $`\widehat{g}^{(2)}`$ being the determinant of the metric of the unit two-sphere. Using (17) one can calculate $`_{DBI}`$ as

$`_{DBI}`$ $`=`$ $`\sqrt{(|g_{00}|g_{rr}\dot{r}^2g_{\rho \rho }\dot{\rho }^2g_{\varphi \varphi }\dot{\varphi }^2)g_{\theta ^1\theta ^1}g_{\theta ^2\theta ^2}(g^3+H_{345}^2)}`$ (18)
$`=`$ $`fr^3\rho ^2\sqrt{\widehat{g}^{(2)}}\lambda _1\left[r^2f^1r^2\dot{r}^2{\displaystyle \frac{\dot{\rho }^2}{1\rho ^2}}(1\rho ^2)\dot{\varphi }^2\right]^{1/2},`$
where $`\lambda _1`$ is defined as

$$\lambda _1\sqrt{hf^1+H_{345}^2h^1}.$$
(19)
The remaining terms of the action are calculated as

$$_{H\stackrel{~}{H}}+_{WZ}=\frac{1}{2}F_{345}F_0F_{345}P[C^{(3)}]_0+P[C^{(6)}]_{0345}+\frac{1}{2}P[C^{(3)}]_{345}P[C^{(3)}]_0,$$
(20)
where the index $`0`$ means $`x^0\theta ^1\theta ^2`$. The pullbacks of $`C^{(3)}`$ and $`C^{(6)}`$ are

$`P[C^{(3)}]_0`$ $`=`$ $`R^3\rho ^3\mathrm{cos}\phi \sqrt{\widehat{g}^{(2)}}\dot{\varphi },`$
$`P[C^{(3)}]_{345}`$ $`=`$ $`\mathrm{tan}\phi hf^1,`$ (21)
$`P[C^{(6)}]_{0345}`$ $`=`$ $`{\displaystyle \frac{1}{2}}R^3\rho ^3\mathrm{sin}\phi hf^1\sqrt{\widehat{g}^{(2)}}\dot{\varphi }.`$
Substituting Eq. (21) in Eq. (20) the last two terms cancel each other, and we have

$$_{H\stackrel{~}{H}}+_{WZ}=R^3\rho ^3F_{345}\mathrm{cos}\phi \sqrt{\widehat{g}^{(2)}}\dot{\varphi }+\frac{1}{2}F_{345}F_0.$$
(22)
We assume that $`F_0=\sqrt{\widehat{g}^{(2)}}f_0`$ with $`f_0`$ being independent of the angle of the $`\mathrm{S}^2`$. With this ansatz for the electric component of $`F`$, we can integrate out $`\theta _1`$ and $`\theta _2`$ using

$$\sqrt{\widehat{g}^{(2)}}𝑑\theta _1𝑑\theta _2=4\pi \mathrm{\Omega }_2.$$
(23)
Then the action can be reduced to the following four-dimensional (three-brane) effective action

$$S=𝑑t𝑑x^3𝑑x^4𝑑x^5,$$
(24)
with

$`=\mathrm{\Omega }_2T_{M5}\{`$ $``$ $`\sqrt{|g_{00}|g^3g_\theta ^2}\left[1{\displaystyle \frac{g_{rr}}{|g_{00}|}}\dot{r}^2{\displaystyle \frac{g_{\rho \rho }}{|g_{00}|}}\dot{\rho }^2{\displaystyle \frac{g_{\varphi \varphi }}{|g_{00}|}}\dot{\varphi }^2\right]^{\frac{1}{2}}\left[1+{\displaystyle \frac{H_{345}^2}{g^3}}\right]^{\frac{1}{2}}`$ (25)
$`+`$ $`R^3\rho ^3F_{345}\mathrm{cos}\phi \dot{\varphi }+{\displaystyle \frac{1}{2}}F_{345}f_0\},`$
where $`g_\theta =g_{\theta ^1\theta ^1}=f^{2/3}h^{1/3}r^2\rho ^2`$.

## IV Brane cosmology

Since we are interested in the cosmological evolution in terms of $`r`$, we consider the case when $`\rho =\mathrm{constant}`$, i.e. $`\dot{\rho }=0`$. This corresponds to the case when the probe universe is planar. In this case we can rewrite $``$ as

$$=\mathrm{\Omega }_2T_{M5}\left\{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}+G\dot{\varphi }+\frac{1}{2}F_{345}f_0\right\},$$
(26)
where

$`A`$ $`=`$ $`|g_{00}|g_\theta ^2(g^3+H_{345}^2)=fr^4\rho ^4\lambda _1^2,`$
$`B`$ $`=`$ $`g_{rr}g_\theta ^2(g^3+H_{345}^2)=f^2r^4\rho ^4\lambda _1^2,`$
$`D`$ $`=`$ $`g_{\varphi \varphi }g_\theta ^2(g^3+H_{345}^2)=f^2r^6\rho ^4(1\rho ^2)\lambda _1^2,`$
$`G`$ $`=`$ $`R^3\rho ^3F_{345}\mathrm{cos}\phi .`$ (27)
The momenta and hamiltonian, divided by the overall factor $`\mathrm{\Omega }_2T_{M5}`$, are calculated as

$`p_r`$ $`=`$ $`{\displaystyle \frac{}{\dot{r}}}={\displaystyle \frac{B(r)\dot{r}}{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}}},`$
$`p_\varphi `$ $`=`$ $`{\displaystyle \frac{}{\dot{\varphi }}}={\displaystyle \frac{D(r)\dot{\varphi }}{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}}}+G,`$ (28)
$``$ $`=`$ $`\dot{r}p_r+\dot{\rho }p_\rho +F_0{\displaystyle \frac{}{F_0}}={\displaystyle \frac{A(r)}{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}}}.`$
We require the conservation of energy as well as the angular momentum

$``$ $`={\displaystyle \frac{A(r)}{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}}}=E=\mathrm{const},`$ (29)
$`p_\varphi `$ $`={\displaystyle \frac{D(r)\dot{\varphi }}{\sqrt{A(r)B(r)\dot{r}^2D(r)\dot{\varphi }^2}}}+G=\mathrm{}=\mathrm{const}.`$ (30)
If we solve Eqs. (29) and (30) for $`\dot{\varphi }`$ and $`\dot{r}`$, we have

$`\dot{\varphi }^2`$ $`=\left({\displaystyle \frac{A}{D}}\right)^2\left({\displaystyle \frac{\mathrm{}G}{E}}\right)^2,`$ (31)
$`\dot{r}^2`$ $`={\displaystyle \frac{A}{B}}\left\{1{\displaystyle \frac{A}{E^2}}{\displaystyle \frac{D+(\mathrm{}G)^2}{D}}\right\}.`$ (32)
Since $`\dot{r}^20`$ , we have the constraint for the allowed values of $`r`$

$$\frac{A}{B}\left\{1\frac{A}{E^2}\frac{D+(\mathrm{}G)^2}{D}\right\}0.$$
(33)
Using the expressions in Eq. (27), we can estimate the range of $`r`$ where our formalism is valid

$$\frac{r}{R}\frac{E^2R^2(1\rho ^2)}{(\mathrm{}R^3\rho ^3\mathrm{cos}\phi F_{345})^2+\rho ^4(1\rho ^2)R^6\mathrm{cos}^2\phi F_{345}^2}\frac{r_c}{R}.$$
(34)
The induced metric on the three-brane universe (= 5-brane/$`\mathrm{S}^2`$) can be written as

$$ds_{4d}^2=(|g_{00}|+g_{rr}\dot{r}^2+g_{\varphi \varphi }\dot{\varphi }^2)dt^2+g(r)[(dx^3)^2+(dx^4)^2+(dx^5)^2].$$
(35)
Using Eqs. (31) and (32), this reduces to

$$ds_{4d}^2=|g_{00}|\frac{A}{E^2}dt^2+g(r)[(dx^3)^2+(dx^4)^2+(dx^5)^2]d\eta ^2+g(r(\eta ))(d\stackrel{}{x})^2,$$
(36)
where we defined the cosmic time $`\eta `$ as

$$d\eta =\frac{|g_{00}|^{1/2}A^{1/2}}{E}dt=\frac{|g_{00}|g^{3/2}g_\theta (1+\frac{H_{345}^2}{g^3})^{1/2}}{E}dt.$$
(37)
If we define the scale factor as $`a^2g`$, we can calculate, from the analogue of the four-dimensional Friedman equation, the Hubble constant $`H=\dot{a}/a`$

$$\left(\frac{\dot{a}}{a}\right)^2=\frac{1}{4|g_{00}|}\left(\frac{E^2}{B}\frac{A}{B}\frac{A}{D}\frac{(\mathrm{}G)^2}{B}\right)\left(\frac{g^{}}{g}\right)^2,$$
(38)
where the dot denotes the derivative with respect to cosmic time $`\eta `$ and the prime denotes the derivative with respect to $`r`$. The right hand side of Eq. (38) can be interpreted as the effective matter density on the probe 3-brane. Upon substituting the specific forms of $`B,C,D`$ and $`G`$ of Eq. (27) we have

$$\frac{8\pi }{3}\rho _{\mathrm{eff}}=\frac{1}{4|g_{00}|}\left(\frac{E^2}{\rho ^4r^4f^2\lambda _1^2}\frac{1}{f}\frac{1}{(1\rho ^2)r^2f}\frac{(\mathrm{}G)^2}{\rho ^4r^4f^2\lambda _1^2}\right)\left(\frac{g^{}}{g}\right)^2.$$
(39)
Defining the dimensionless variable $`x`$ as $`xr/R`$, we can write the effective matter density explicitly as

$`{\displaystyle \frac{8\pi }{3}}\rho _{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{4R^2(\mathrm{cos}^2\phi )^{1/3}}}{\displaystyle \frac{\left\{1+(1\mathrm{tan}^2\phi )x^3\right\}^2}{(1+x^3)^{7/3}(1+\mathrm{sec}^2\phi x^3)^{7/3}}}`$ (40)
$`\times `$ $`\left\{{\displaystyle \frac{E^2\mathrm{cos}^2\phi }{\rho ^4R^4x^4}}k(x,\phi ,F)1{\displaystyle \frac{(\mathrm{}R^3\rho ^3F_{345}\mathrm{cos}\phi )^2\mathrm{cos}^2\phi }{\rho ^4(1\rho ^2)R^6x^3(1+x^3)}}k(x,\phi ,F)\right\},`$
where

$$k(x,\phi ,F)=[12\mathrm{sin}\phi \mathrm{cos}\phi F_{345}+\mathrm{cos}^2\phi (1+\mathrm{cos}^2\phi x^3)F_{345}^2]^1.$$
(41)
To obtain more transparent expression for the cosmological evolution, we express the effective matter density in terms of scale factor $`a`$ as

$`{\displaystyle \frac{8\pi }{3}}\rho _{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{4R^2}}{\displaystyle \frac{a(f1)^{8/3}}{f^{5/2}}}(12a^3f^{1/3}\mathrm{cos}^2\phi )^2`$ (42)
$`\times `$ $`[{\displaystyle \frac{E^2}{\rho ^4R^4}}{\displaystyle \frac{(f1)^{4/3}}{a^3f^{1/2}}}\{1+({\displaystyle \frac{F_{345}}{a^3}}{\displaystyle \frac{\mathrm{tan}\phi }{f^{1/2}}})^2\}^11`$
$``$ $`{\displaystyle \frac{(\mathrm{}R^3\rho ^3F_{345}\mathrm{cos}\phi )^2}{(1\rho ^2)\rho ^4R^6}}{\displaystyle \frac{(f1)^2}{a^3f^{3/2}}}\{1+({\displaystyle \frac{F_{345}}{a^3}}{\displaystyle \frac{\mathrm{tan}\phi }{f^{1/2}}})^2\}^1],`$
where

$$f=\frac{12a^6\mathrm{cos}^2\phi \mathrm{sin}^2\phi +\sqrt{14a^6\mathrm{cos}^2\phi \mathrm{sin}^2\phi }}{2a^6\mathrm{cos}^4\phi }.$$
(43)
Let us consider a limiting case to compare our expression with the result from type IIB background. We consider the case when there is no gauge field on the worldvolume, i.e., $`F_{345}=0`$. In this case the effective matter density is given by, taking the leading powers of the scale factor,

$$\frac{8\pi }{3}\rho _{\mathrm{eff}}\frac{1}{4R^2(\mathrm{cos}^2\phi )^{1/3}}\left[\frac{E^2}{\rho ^4R^4(\mathrm{cos}^2\phi )^{5/3}}\frac{1}{a^8}1\frac{\mathrm{}^2}{\rho ^4(1\rho ^2)R^6\mathrm{cos}^2\phi }\frac{1}{a^6}\right].$$
(44)
Near the horizon, the effective matter density is proportional to $`\rho _{\mathrm{eff}}a^8`$, which shows the same power behavior as the result from the type IIB background without any gauge field on the worldvolume Keki . Also the $`\mathrm{}^2`$ term has the same sign and power behavior.

## V Discussion

We searched the possibility of constructing the mirage cosmology in M-theory background. We considered the motion of a five-brane in the background formed by a stack of non-threshold (M2,M5) bound states. From the brane action in the PST formalism, we derived a Friedman-like equation. We took a limiting case and compared the result with the one from type IIB background. We conclude that the cosmological evolution by mirage energy density is a possible scenario in M-theory background.

As estimated in Eq. (34), our formalism on mirage cosmology holds for $`rr_c`$. But this does not mean we can extend our result to the initial singularity where the effect of the back reaction is important. In mirage cosmology the initial singularity appears not because the solution is singular but because the effective field theory is not valid in this region. It is just an artifact of the low energy description Keki . The cosmological evolution from our result can be summarized as follows. When the probe brane is near the (M2,M5) bound states, the probe brane expands mainly due to mirage energy density. In this region the universe expands very rapidly ($`\rho _{\mathrm{eff}}a^8`$). As the brane universe moves away from the background bound states ($`rr_c`$), the effect of the background branes to the probe brane will not be strong enough to drive the inflation. Then our formalism on mirage cosmology is not valid any more. In this region the matter density of the probe universe itself will drive the cosmological expansion and the rate of expansion will be slower than the one by mirage energy density.

Although it is an important open problem how to study the back reaction of the probe brane, we did not consider the back reaction of the probe brane to the background geometry. When $`\mathrm{}^2`$ term dominates the effective density is negative and we have contraction rather than inflation. We hope this fact might be improved if we consider the back reaction. In our presentation, we considered the motion of probe brane with constant $`\rho `$ ($`\dot{\rho }=0`$). It would be an interesting topic if one studies the mirage cosmology with constant $`r`$ ($`\dot{r}=0`$). In this case, we expect that one could construct the mirage cosmology of closed universe similar to the case in Ref. youmclo . It is known that a Friedman-type evolution in brane cosmology is equivalent to the formalism of varying speed of light vsl . One can also study this model in this context.

###### Acknowledgements.

We would like to thank S.P. Kim for useful discussion. This work was supported by the Korea Science and Engineering Foundation, No. 2001-1-11200-003-2.


text:
# Using the High Resolution X-ray Spectrum of PSR B0656+14 to Constrain the Chemical Composition of the Neutron Star Atmosphere

## 1 Introduction

The X-ray spectra of many “middle-aged” and isolated neutron stars are dominated by a thermal component, so, besides temperatures and radii, they can in principle be used to measure the atmospheric composition. See Becker & Trümper (1997) for a summary of ROSAT observations of several neutron stars whose X-ray emission could be fitted to blackbody spectra. The composition can, in turn, affect the estimated temperature and size of the emission region (Rajagopal & Romani, 1996; Rutledge et al., 1999) but ROSAT spectra were insufficient to determine if there were spectral features that might be expected if cyclotron lines or heavy elements dominate the atmosphere. For example, Rajagopal et al. (1997) showed that an iron-dominated atmosphere in a magnetic field as strong as $`10^{12}`$ G would show many narrow absorption features due to the distortion of the energy levels.

PSR B0656+14 is an isolated radio pulsar with a period of about 384 ms, found in the Einstein ultrasoft X-ray survey by Cordova et al. (1989). ROSAT observations showed weak pulsations with an amplitude of $``$ 14% (Finley et al., 1992) and that the X-ray spectrum appeared to be well fit with a simple black body model with $`T=9\times 10^5`$ K. From longer ROSAT observations, Possenti et al. (1996) found that the spectrum appeared slightly more complex and that the pulse shape changes with energy. Possenti et al. (1996) fitted a two component model to the pulse light curve and X-ray spectrum, suggesting that the hotter component is emission from the poles and the cooler part comes from the equator but the hotter component could instead be modelled by a steep power law component with $`\alpha =3.5\pm 0.4`$ (where $`f_\nu \nu ^\alpha `$). Edelstein et al. (2000) found a somewhat lower temperature by combining EUVE, ROSAT, and optical data. Greiveldinger et al. (1996) required two black body components as well as a hard power law ($`\alpha =1.5\pm 1.1`$) to describe the ASCA spectra extending to 5 keV. The distance estimated from the pulsar’s dispersion measure is 760 pc (Taylor, Manchester, & Lyne, 1993).

We observed PSR B0656+14 using the Chandra Low Energy Grating Spectrometer (LETGS) in order to detect narrow X-ray lines in absorption that might be observed if the atmosphere is dominated by iron. The energies of these absorption features should depend on the average magnetic field strength, so we used the High Resolution Camera Spectroscopy (HRC-S) detector which provide event timing to better than 10 $`\mu `$s in order to obtain pulse phased spectra. We find no strong spectral lines and show that the LETGS data are consistent with the previous results; the continuum is well fit by a model involving two black body components. Pulse phased spectra are used to search for spectral features that might be phase dependent.

## 2 Observations and Data Reduction

### 2.1 Imaging

PSR B0656+14 was observed with the LETGS on 28 November 1999 (JD - 2451000. = 510.99 - 511.43). The level 0 were processed with CIAO version 2.0b.111CIAO is the Chandra Interactive Analysis of Observations, a software system developed by the Chandra X-ray Center. The exposure time was 38167 s. The count rate in zeroth order was consistent with 0.19 count/s for the entire observation. The bandpass for the zeroth order image is hard to define due to the low energy resolution of the HRC-S, so we estimated the energy range over which one may integrate the spectrum-weighted effective area to obtain 90% of the observed count rate: 0.20 - 0.75 keV. To obtain 99% of the total count rate, the bandpass should be extended to 0.16-1.0 keV. The resultant image was azimuthally symmetric, as expected for a point source.

We tested the radial distribution of events against that of a comparable LETG/HRC-S observation of a point source: Capella (observation ID 1248). Due to the dispersion by the fine and coarse support structure, the point response function of the LETGS zeroth order image depends slightly on the spectrum of the source. The zeroth order image yielded 180440 events within 5′ of the centroid in 85260 s. The spectrum of Capella is soft but is harder than that of the pulsar, so we divided the data into two equal parts based on the pulse height distribution. The average pulse height for the low half is quite similar to that of PSR B0656+14, but the radial profile differs by no more than 1.5% from the average profile using all events so we used the average profile for comparison to the pulsar. The encircled power was modelled as the sum of Gaussian and exponential components, representing the Capella data to better than 1.5% over the 0-5″ range:

$$\mathrm{\Phi }(\theta )=(1f)[1e^{\theta ^2/\sigma _{PSF}^2/2}]+f[1e^{\theta /\theta _0}],$$
(1)
where $`\theta `$ is the angle from the centroid, $`f=0.31`$ is the fraction of the power in the exponential component, $`\sigma _{PSF}`$ is the Gaussian width parameter of the point spread function (PSF), and $`\theta _0`$ is the scaling length of the exponential component. Capella is very bright, so the background has a negligible contribution ($`<`$ 1%) to the radial profile out to 5″ from the source centroid. For the pulsar, background was estimated from an annulus 8-10″ from the source centroid.

The cumulative event radial profiles are given in Fig. 1. The profiles are normalized to unity at 5″. The pulsar’s radial distribution differs from Capella’s by $``$ 8.5% at 0.5″, indicating that the X-ray source may consist of a point source and some distributed emission that is extended on a scale of order 0.5″. A Smirnov test indicates that a 5% difference between the Capella and PSR B0656+14 profiles would be significant at the $`>8\sigma `$ level and that differences of 1.6% are significant at the 2$`\sigma `$ level. We modelled the profile by combining a Gaussian model of the extended emission with a PSF model. The model of the extended emission was constructed using the PSF model but increasing the Gaussian width parameter by combining a second term in quadrature, $`\sigma _{ext}`$. The model for the pulsar’s radial distribution is matched to about 2% when $`\sigma _{ext}=0.32`$″ and the extended component comprises 50% of the total power. The uncertainties are difficult to characterize and are probably dominated by systematic uncertainties in the PSF. Given the possible systematic uncertainties in these early HRC-S imaging data, we estimate that an unresolved source accounts for $``$ 50% of the profile and that the FWHM of the extended emission is in the range of 0.2-0.5″.

### 2.2 Timing

The zeroth order image was used to determine the phase and period of the X-ray pulsations. After applying timing corrections to the solar system barycenter and correcting for the satellite’s position in orbit, we determined that the pulse period was 384.8990 ms with a 1 $`\sigma `$ uncertainty of 0.0010 ms using a $`\chi ^2`$ test on folded light curves. This period is consistent with the predicted value, 384.89970 ms from radio pulse monitoring (Andrew Lyne, 2001, priv. communication). The pulse light curve shown in Fig. 2 was folded using the radio ephemeris. The pulse appears to be asymmetric, as previously found by Possenti et al. (1996). The modulation amplitude is 10 $`\pm `$ 2%, defined as $`(R_{max}R_{min})/(R_{max}+R_{min})`$, where $`R_{max}`$ and $`R_{min}`$ are the maximum and minimum of the pulse light curve, respectively. This value is consistent with the pulse fractions obtained by Finley et al. (1992) (14 $`\pm `$ 2%) and Possenti et al. (1996) (9 $`\pm `$ 1%).

The radio phase is defined so that the peak of the radio pulse occurs at a phase of zero. Fig. 2 shows that the X-ray peak is at about phase 0.85, which is about 55° out of phase with the radio phase. The absolute timing HRC-S has been verified to an accuracy of $`<`$ 0.001 s using the Crab pulsar (Allyn Tennant, 2001, priv. comm.).

### 2.3 Spectroscopic Data Reduction

The spectral data were reduced from standard event lists using IDL using custom processing scripts; the method is quite similar to standard processing using CIAO. The procedure was to: 1) select data based on the criteria suggested by the Chandra X-ray Center (CXC) LETGS calibration team222Chandra calibration information is available at the CXC web site:http://cxc.harvard.edu/cal/. The HRC-S selection recommendation can be found at: http://cxc.harvard.edu/cal/Links/Letg/User/Hrc\_bg/. in order to reduce background, 2) determine the location of zeroth order using one dimensional profiles fitted to Gaussians, 3) rotate events from sky coordinates to compensate for the telescope roll and correct for a 0.54″ offset between HRC-S plates 0 and 1, 4) compute the dispersed grating coordinates ($`m\lambda `$ and $`\theta `$) using the grating dispersion angles and the dispersion relation, 5) select “source” events spatially within oppositely curving parabolas of the dispersion line as suggested by the CXC LETG calibration team 333See http://cxc.harvard.edu/cal/Links/Letg/User/Hrc\_QE/EA/Wads. , 6) eliminate data affected by detector gaps, and 7) bin events at $`\mathrm{\Delta }\lambda =`$0.025Å.

The effective area (EA) of the LETGS has undergone a few revisions since launch both due to a recalculation of the LETG efficiencies and due to a long series of in-flight calibration observations designed to probe the HRC-S quantum efficiency. We started with the updated LETGS EA that was released by the CXC of 31 October 2000444 This effective area is available at http://cxc.harvard.edu/cal/Links/Letg/User/Hrc\_QE/EA/correct\_ea/letgs\_NOGAP\_EA\_001031.mod.. These LETG efficiencies were used, along with the transmission models of the UV ion shield, to determine the effective areas for orders 2 through 5. In a final step, an adjustment to the pre-flight calibration was applied which was derived from an analysis of the spectrum of Mrk 478 (Marshall et al., 2002) and PKS 2155-304 (Marshall et al. 2002, in preparation). The spectra of both targets are well fit by a simple power law with Galactic absorption. The adjustments affect primarily the spectrum at low energies ($`E<0.2`$ keV), where high order contributions are beginning to be important.

Directly totalling the observed fluxes over the 0.2-3.0 keV band gives an observed flux of 1.01 $`\pm 0.03\times 10^{11}`$ erg cm-2 s-1 for an absorbed luminosity of 7.0 $`\pm 0.2\times 10^{32}`$ erg s-1. The data were rebinned adaptively to provide a signal/noise ratio of 5 in each bin over the 0.10 to 2.0 keV range. The spectrum, shown in Fig. 3, was estimated using the first order EA only. The contributions to the observed counts due to high orders are estimated by folding a model for first order through the high order EA and dividing by the first order EA. This procedure has the disadvantage that all the spectral features in the high order EAs are apparent but the advantage is that the result will match the intrinsic spectrum well if first order dominates. For this source, high orders are negligible above 0.20 keV.

## 3 Modelling the Spectrum

Following previous analyses, we modelled the continuum with two blackbody components. Our objective is primarily to define a smooth continuum model to use as a baseline for line searches. We did not include the hard power law component found by Greiveldinger et al. (1996) because its contribution in the 0.2-1.0 keV band is negligible. We exclude the data below 0.15 keV from the fit where uncertainties in the high order grating efficiencies can be important. Gaussian statistics could be used for continuum modelling by rebinning to obtain a signal/noise ratio of at least 5 in each bin, giving reduced $`\chi ^2=1.13`$, acceptable at the 90% level. The distribution of the residuals is not consistent with a Gaussian, however; the Kolmogorov test rejects at the 99.99% confidence level due to systematic skewing of the residuals to positive values. This may result from systematic errors at the 5-10% level in the 0.45-0.55 keV portion of the spectrum.

The best fit temperatures of the two blackbody components were $`8.0\pm 0.3\times 10^5`$ and $`1.6\pm 0.3\times 10^6`$ K and the radii were 22.5 $`\pm `$ 2.1 and 1.7 $`\pm `$ 1.0 km, respectively. The best-fit interstellar medium (ISM) column density, $`N_H`$, was $`1.73\pm 0.18\times 10^{20}`$ cm-2. These parameters are close to those found by Greiveldinger et al. (1996). The best fit model is plotted against the fitted data in Fig. 3. The count spectrum (Fig. 4) was binned at 0.125 Å resolution in order to search for narrow spectral features against the continuum model. No significant absorption features were found. A fit to a nonmagnetized H atmosphere (Zavlin, Pavlov, & Shibanov, 1996) gave a much lower temperature, 230000 K, a very large radius, 413 $`\pm `$ 53 km and a much larger $`N_H`$, $`2.6\pm 0.2\times 10^{20}`$ cm-2.

The data were divided into pulse maximum and minimum by phase: 0.05-0.55 and 0.55-1.05 (see Fig. 2). The difference spectrum (Fig. 5) shows a dip in the 45-50 Å region. The pulse amplitude spectrum was defined as the ratio of the difference spectrum to the model of the pulse-averaged spectrum and is also shown in Fig. 5. The result is rather flat, consistent with being constant over this wavelength band.

## 4 Discussion and Summary

Models of a neutron star atmosphere dominated by Fe in a strong magnetic field, such as computed by Rajagopal et al. (1997), should show deep narrow features which we do not observe. For a hydrogen atmosphere, however, there are very few spectral features but the radius of the fitted nonmagnetized H atmosphere is much too large. In order to reduce the apparent radius to 13 km (appropriate for a neutron star with a mass of 1.4 $`M_{\mathrm{}}`$ and a radius of 10 km), the distance would have to be absurdly small – 24 pc – given the large $`N_H`$ and the large distance from the dispersion measure. For the blackbody model, however, the distance would have to be about 440 pc, which is consistent with the estimate from the pulsar’s dispersion measure. For a magnetized H atmosphere, results are similar to the case with $`B=0`$ but are not quite so extreme. Meyer et al. (1994) found that the ROSAT spectrum of Geminga fit a blackbody with $`T=7.6\times 10^5`$ K while a magnetized H atmosphere gave $`T=5.0\times 10^5`$ K for $`B=4.7\times 10^{12}`$ G. The magnetic field estimated for PSR B0656+14 is the same (Taylor, Manchester, & Lyne, 1993), so one might expect a similar reduction when fitting the PSR B0656+14 spectrum. In order to generate the same observed flux, the radius of the emission region would have to increase by a factor of 2.3 to 52 $`\pm `$ 5 km so the distance would have to be reduced to 190 pc. A distance to less than 200 pc can be definitely ruled out using the dispersion measure and a more refined model of the local ISM (J. Cordes, 2002, private communication).

We find no significant absorption features in the pulse-averaged or pulse-phased spectra over the 0.15-0.80 keV band. Thus, we can rule out electron and proton cyclotron resonance lines, which should be rather deep and possibly broadened, according to models by Zane et al. (2001) and Ho & Lai (2001). The proton cyclotron line would be expected at an observed energy $`E_o=E_p/(1+z)=\mathrm{}eB/(m_pc)/(1+z)=0.0063/(1+z)B_{12}`$ keV, where $`B_{12}=B/(10^{12}`$ G) and $`z=(12GM/Rc^2)^{1/2}1`$ is the surface redshift. Substituting the electron mass for the proton mass, $`m_p`$, gives the electron cyclotron line energy $`E_e/(1+z)=\mathrm{}eB/(m_ec)/(1+z)=1.16/(1+z)B_{12}`$ keV. Thus, we eliminate the ranges $`0.017<B_{12}<0.090`$ and $`30<B_{12}<165`$ if $`z=0.30`$ (for a neutron star of 1.4 M and a radius of 10 km). These ranges do not exclude the magnetic field estimate from the dipole braking model: $`B_{12}=4.7`$ (Taylor, Manchester, & Lyne, 1993). The lower limits can be reduced if the bulk of the emission comes from the equatorial zone where the magnetic field is $`50`$% of the polar value or if the absorbing plasma is far off of the neutron surface.

For this dipole field estimate, a light element atmosphere should show H or He ionization edges in the 0.15-0.80 keV bandpass. Lai & Salpeter (1997) give an estimate for the ionization edge energy and its dependence on the atomic number $`Z`$:

$$E_Z=4.4Z^2[\mathrm{ln}(\frac{B}{Z^2B_0})]^2\mathrm{eV}$$
(2)
where $`B_0=2.35\times 10^9`$ G. A He i edge would be expected in the 0.5-0.7 keV range. We do not detect this edge, which may not be surprising, as helium may have settled in the atmosphere to high optical depths. For surface temperatures of 9-22 $`\times 10^5`$ K (allowing for variation between the pole and the equator), He is almost completely neutral while H is 7-22% ionized. So, we expect to find the H i edge between 0.20 and 0.25 keV (depending on the surface gravity) but it is not observed. The spectrum is dominated by the cooler thermal component, however, which is presumably emission from the equatorial zone in which $`B`$ is $`2\times `$ weaker than at the pole. Thus, the edge may be found in the 0.16-0.21 keV region, where the signal/noise ratio is poor and high order contributions to the LETGS spectrum become important.

Most of a pulsar’s spin-down luminosity is thought to be carried off as a relativistic electron/positron wind. When this wind is confined by external pressure, a shock forms and the relativistic particles radiate synchrotron emission in the form of a pulsar wind nebula (PWN). The Crab Nebula (Hester et al. 1995, Weisskopf et al. 2000) is the best known example. PWNe have been observed around at least 6 pulsars so far (Becker & Trümper, 1997; Lu et al., 1997). Claims from ASCA observations of vastly extended PWNe ($`>`$ 10-20′) in this and other pulsars with luminosities of the order of $`10^{32}`$ erg s-1 have not been confirmed by ROSAT and BeppoSAX (Becker et al., 1999). We also do not detect such a large extended nebula around PSR B0656+14. The image appears slightly extended, consistent with a point source comprising $`50`$% of the total 0.1-3 keV luminosity. The remainder, with an unabsorbed luminosity of $`3\times 10^{32}`$ erg s-1 can be modelled with a Gaussian with a FWHM of 0.75″ corresponding to a size of 8.5 $`\times 10^{15}`$ cm. This extended emission is consistent with the result of Becker et al. (1999), who found that a PWN around PSR B0656+14 must have an extent of less than 10″. Based on its spin-down age, $`10^5`$ yr, PSR B0656+14 is an order of magnitude older than the Vela pulsar, the oldest pulsar around which a PWN has been detected (Becker & Trümper, 1997). If we assume a similar emission efficiency for PSR B0656+14 as for the Vela pular ($`0.04\%`$), we can set an upper limit to the expected luminosity of the PWN of 1.5$`\times 10^{33}`$ erg s-1, consistent with the observed value.

We are extremely grateful to Prof. Claude R. Canizares, the Principal Investigator of the Chandra High Energy Transmission Grating Spectrometer, for allocating part of his guaranteed time to this observation. We thank the referee for comments that have resulted in significant improvements to the paper. This work has been supported in part under NASA contract SAO SV1-61010.


text:
# Constraints on Neutrino Degeneracy from the Cosmic Microwave Background and Primordial Nucleosynthesis

## I Introduction

The present relic neutrino number asymmetry is not directly observable. Hence, there is no firm experimental basis for postulating that the lepton number for each species is zero. Charge neutrality, however, demands that any universal net lepton number beyond the net baryon number must reside entirely within the neutrino sector. It has been suggested that the total lepton number could be large in the context of the $`SU(5)`$ and $`SO(10)`$ grand unified theories , or supersymmetric baryogenesis  based upon the Affleck-Dine scenario . Such mechanisms might generate lepton number asymmetry up to ten orders of magnitude larger than the baryon number asymmetry. Furthermore, even if one demands that $`BL0`$, it is possible for the lepton numbers $`L_l`$ of individual neutrino species to be large compared to the baryon number of the universe, $`B`$, as long as the net total lepton number is small.

Moreover, there presently exists at least some marginal cosmological evidence for neutrino asymmetry. For example, neutrinos with large lepton asymmetry and masses $`0.07`$ eV might be required to explain the existence of cosmic rays with energies in excess of the Greisen-Zatsepin-Kuzmin cutoff . Also, degenerate, massive (2.4 eV) neutrinos might be required  to provide a good fit to the power spectrum of large scale structure in mixed dark matter models. It is thus important to carefully scrutinize the limits which cosmology places on the allowed values of possible universal lepton asymmetry. Indeed, a number of recent papers  have addressed this issue with varying degrees of complexity. The present work differs from those in several details as summarized below. It represents an independent examination of this issue.

## II Present Approach

In a recent paper  we considered new constraints imposed on neutrino degeneracy from primordial nucleosynthesis. Particular attention was paid to the neutrino decoupling temperatures before the nucleosynthesis epoch. Of relevance to the present work is that we have shown that neutrinos can decouple at a higher temperature than estimated in earlier studies . This means that more particle degrees of freedom could be present at neutrino decoupling. This causes the relic neutrino temperature to be lower by simple entropy considerations. A smaller relic neutrino energy density implies that larger neutrino degeneracies may be allowed. For example, we have shown that interesting regions of the model parameters for big-bang nucleosynthesis (BBN) are allowed such that substantial lepton asymmetry and baryon density (even $`\mathrm{\Omega }_bh^21`$ where $`h`$ is the present value of the Hubble constant in units of 100 km sec-1 Mpc-1) are possible while still satisfying the adopted abundance constraints from primordial nucleosynthesis.

A stronger constraint, however, on lepton asymmetry comes from the power spectrum of fluctuations in the cosmic microwave background (CMB) which we now address in the present paper. We apply a likelihood analysis of neutrino-degenerate models to the combined latest BOOMERANG , DASI  and MAXIMA-1  results. We note, however, that a recent analysis  of the implications of neutrino oscillations derived from a combination of the atmospheric, and solar neutrino constraints implies much tighter limits on degeneracy for all neutrino flavors. If neutrino oscillation parameters are in the range of the large mixing angle solution then an upper limit of $`|\xi _{\nu _i}|_{}^<0.07`$ applies to all neutrino flavors. The limits derived here do not assume any particular model for neutrino mixing, and should be taken as independent of, and complementary to, those constraints.

The implications of the CMB data for neutrino-degenerate cosmologies have been noted in a number of recent papers . The constraints on the effective number of relativistic particles can also arise in other contexts, such as cosmic quintessence . The present work, however, differs from those in several respects. For one, we consider the most recent combined data sets, not just the first year BOOMERANG data as in . and generate a marginalized likelihood function for the neutrino degeneracy and other cosmological parameters. Many of the existing studies have marginalized over a more limited set of cosmological parameters. For example, in  only $`\Omega _\mathrm{b}h^2`$ and and neutrino degeneracy were marginalized to set limits while other cosmological parameters were set to various fixed values. In , for example, no likelihood analysis was made. In  a likelihood analysis was made but without window functions.

For the present work we use the updated Radpack package  described below which includes all relevant window functions. Another difference between our analysis and other works is that our marginalization utilizes a global minimum search algorithm  rather than a discrete grid of cosmological parameters. Marginalizing parameters for each fixed value of one parameter requires at least 1000 model calculations to get $`10^4`$ accuracy for the $`\chi ^2`$ minima even by using this algorithm. Nevertheless, in this way we are sure to identify the true marginalized likelihood functions.

The most similar recent likelihood analysis to that described here is in the work of . Our analysis differs from  in several respects. In the present work we make use of our deduced new family of baryon densities and lepton asymmetries allowed by BBN to assign weak top-hat priors on the derived likelihood functions. This differs from that in  in which separate Gaussian likelihood functions were evaluated for the nucleosynthesis constraints and the CMB. A total likelihood function was then defined by marginalizing over the product of these two functions. We prefer our method because the uncertainties in the BBN constraints are dominated by systematic errors. Systematic errors are not equivalent to random Gaussian errors. We thus prefer weak top-hat priors as a more realistic representation of the systematic errors in the BBN constraints.

One other important difference is that we adopt a strong Gaussian prior of $`h=0.72\pm 0.08`$ based upon the Hubble Key Project results . In  weak a top-hat prior of $`h=0.65\pm 0.20`$ was adopted. As noted above, another difference between the present work and all previous results is that we consider carefully the change in background neutrino temperatures as a function of degeneracy. Although this is a small effect for low degeneracies, it can slightly affect the upper limits.

## III Neutrino-Degenerate BBN

Neutrino degeneracy affects BBN in two ways. The inclusion of a small amount of electron $`\nu _e`$ degeneracy resets the equilibrium neutron to proton ratio at weak-reaction freezeout to $`n/p=\mathrm{exp}\{\mathrm{\Delta }m/T\xi _{\nu _e}\}`$. This can cause a reduction in the primordial helium abundance. Indeed, it has been argued  that the apparent conflict between the low helium abundance inferred from HII regions of metal poor galaxies and the low Lyman-$`\alpha `$ deuterium abundance may even require $`\nu _e`$ degeneracy for its resolution. The present deuterium-absorption limits on $`\mathrm{\Omega }_bh^20.020\pm 0.002(2\sigma )`$ requires a large primordial helium abundance of $`Y_p{}_{}{}^{>}0.25`$ and substantial destruction of primordial 7Li in stars. These conditions tax even the most generous adopted limits from observed light-element abundances . Thus, a modification of BBN which allows for large values of $`\mathrm{\Omega }_bh^2`$ while still satisfying the constraints from light-element abundances is worth investigating. Such conditions are easily satisfied by neutrino-degenerate models.

The inclusion of either $`\nu _\mu `$ or $`\nu _\tau `$ degeneracy on the other hand only enhances the background energy density and therefore the universal expansion rate. During the radiation dominated epoch, relativistic neutrinos contribute a large fraction of the mass energy. Thus, even a small modification of the neutrino energy density can significantly affect the expansion.

The energy density $`\rho _\nu `$ due to degenerate neutrinos (or any other fermions) are described by the usual Fermi-Dirac distribution functions $`f_\nu =[\mathrm{exp}\left(E/T_\nu \xi _\nu \right)+1]^1`$, where the neutrino degeneracy parameter is defined by $`\xi _\nu \mu _\nu /T_\nu `$, and $`\mu _\nu `$ is the neutrino chemical potential. Thus, we have

$$\rho _\nu +\rho _{\overline{\nu }}=\frac{1}{2\pi ^2}_0^{\mathrm{}}𝑑pp^2E_\nu (f_\nu (p)+f_{\overline{\nu }}(p)).$$
(1)
where, $`p`$ denotes the magnitude of the 3-momentum, and $`E_\nu =\sqrt{p^2+m_\nu ^2}`$, with $`m_\nu `$ the neutrino mass. Here and throughout the paper we use natural units( $`\mathrm{}=c=k_B=1`$).

For the present discussion it is sufficient to only consider massless neutrinos. (Possible limits on neutrino-degenerate models with massive neutrinos are considered in ). The energy density in massless neutrinos becomes

$$\rho _\nu +\rho _{\overline{\nu }}=\frac{7}{8}\frac{\pi ^2}{15}\underset{i}{}T_{\nu _i}^4\left[1+\frac{15}{7}\left(\frac{\xi _{\nu _i}}{\pi }\right)^4+\frac{30}{7}\left(\frac{\xi _{\nu _i}}{\pi }\right)^2\right],$$
(2)
from which it is clear that degeneracy in any neutrino species tends to increase the energy density. The associated increased expansion rate tends to increase the neutrino decoupling temperature. This causes an increase in the primordial helium and other light-element abundances.

Since $`\xi _{\nu _\mu }`$ and $`\xi _{\nu _\tau }`$ primarily affect the expansion rate, they are roughly interchangeable as far as their effects on nucleosynthesis or the CMB are concerned. Furthermore, it now seems likely  that the mixing parameters for $`\nu _\mu `$ and $`\nu _\tau `$ involve a large mixing angle and small $`\delta m^2`$. In this case it is plausible that the muon and tau neutrinos were interconverted in the early universe and would therefore obtain nearly an identical degeneracy parameter if an asymmetry exists. Thus, we adopt a conservative model in which the $`\mu `$ and $`\tau `$ neutrinos are equally degenerate, $`|\xi _{\nu _\mu }|=|\xi _{\nu _\tau }|\xi _{\nu _{\mu ,\tau }}`$.

As shown in , for each value of $`\xi _{\nu _{\mu ,\tau }}`$ there is a unique range of $`\xi _{\nu _e}`$ and $`\Omega _\mathrm{b}h^2`$ which satisfies the combined deuterium and primordial helium constraints. The allowed family of neutrino-degenerate models employed in this work is summarized in Figure 1.

This figure differs slightly from the family of allowed solutions given in  in that we have adopted the newer D/H constraint from  (i.e. D/H $`=3.0(\pm 0.4)\times 10^5`$) and slightly different limits on the primordial helium abundance ($`0.228Y_p0.248`$). In the limit of the standard nondegenerate big bang ($`\xi _{\nu _{\mu ,\tau }}=\xi _{\nu _e}=0`$) our limits on $`\Omega _\mathrm{b}h^2`$ reduce to those of , i.e. $`\Omega _\mathrm{b}h^2=0.021\pm 0.002`$. The allowed shaded regions in Figure 1 will be adopted as weak top-hat priors in the CMB likelihood analysis described below. These regions include both the uncertainties from the abundance constraints described above and the uncertainties in the BBN model predictions .

## IV CMB Power Spectrum

Having defined the family of allowed priors from BBN we can now do a likelihood search for optimum cosmological parameters which fit the CMB data. Several recent works  have explained how neutrino degeneracy can dramatically alter the power spectrum of the CMB. For massless neutrinos it can be shown  that the only effect of neutrino degeneracy is to increase the background pressure and energy density of relativistic particles. The essence of this constraint is that degenerate neutrinos increase the energy density in radiation at the time of photon decoupling and delay the time of matter-radiation energy-density equality. This mainly causes an increase in the amplitude of the first acoustic peak in the CMB power spectrum at $`l200`$. For example, based upon a $`\chi ^2`$ analysis  of 19 experimental points and window functions, it was concluded in  that $`\xi _\nu 6`$ for a single degenerate neutrino species with an $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ cosmology.

However, in the existing CMB constraint calculations  only small degeneracy parameters with the standard relic neutrino temperatures were studied in the derived constraint. Hence, the possible effect of a diminished relic neutrino temperature at high degeneracy needs to be considered. To investigate this we have done calculations of the CMB power spectrum, $`\mathrm{\Delta }T^2=l(l+1)C_l/2\pi `$ based upon the CMBFAST code of Seljak & Zaldarriago . We have explicitly modified this code to account for the contribution of massless degenerate neutrinos with varying relic neutrino temperatures $`T_{\nu _i}`$ for each species .

The experimental uncertainties are non-Gaussian, but can be well represented by an offset log-normal distribution . As in  we have evaluated the $`\chi ^2`$ goodness of fit for a range of theoretical power spectra $`C_l`$ as follows: We define the goodness of fit by

$$\chi ^2=\underset{i,j}{}\left(Z_i^t\overline{Z}_i^d\right)M_{ij}^Z\left(Z_j^tZ_j^d\right)+\chi _{cal}^2,$$
(3)
where separate summation over the different data sets is implied. For each set of binned power data $`Z_i^d`$ we utilize the published off-set log-normal data from the three data sets.

$$Z_i^d\mathrm{ln}(D_i+x_i)$$
(4)
Where $`D_i`$ is the measured band power. The corresponding binned theoretical power spectra are

$$Z_i^t\mathrm{ln}\left[\underset{l}{}ϵ_i(W_{il}/l)(𝒞_l+x_i)\right],$$
(5)
where the $`ϵ_i`$ are the published calibration uncertainties taken to be 8%, for MAXIMA-1 and DASI, and 20% for BOOMERANG. Window functions $`W_{il}`$ for the three data sets are available on the world wide web. The error matrix is simply

$$M_{ij}^Z=M_{ij}(D_i+x_i)(D_j+x_j)$$
(6)
where $`M_{ij}`$ is the weight matrix for the band powers $`D_i`$. The effect of the calibration uncertainty on the goodness of fit is obtained from

$$\chi _{cal}^2=\underset{i}{}\frac{(ϵ_i1)^2}{\sigma _i^2},$$
(7)
where $`\sigma _i`$ is the experimental uncertainty. The total $`\chi ^2`$ evaluated in this way can be converted  into a likelihood function for each parameter $`x`$ marginalized over the remaining parameter set $`\stackrel{}{y}`$

$$(x)=P_{prior}(x,\stackrel{}{y})exp(\chi ^2/2)𝑑\stackrel{}{y}.$$
(8)
In neutrino-degenerate models which satisfy the constraints from primordial nucleosynthesis , increasing the neutrino-degeneracy must be accompanied by a commensurate increase in baryon density. Fits to the CMB power spectrum for large degeneracy , therefore show a suppression of the the second acoustic peak due to baryon drag .

Indeed, such suppression of the second acoustic peak seemed to be present in the first reported power spectra based upon the balloon-based CMB sky maps from the BOOMERANG  and MAXIMA-1  collaborations. This remains true for the likelihood analysis based upon a MAXIMA-1 data  which indicates $`\mathrm{\Omega }_bh^2=0.030_{0.010}^{+0.018}(2\sigma )`$. However, in the most recent data sets from BOOMERANG  and DASI  the second peak has become much better defined. Both the BOOMERANG and DASI data sets now imply $`\mathrm{\Omega }_bh^2=0.022_{0.003}^{+0.004}(1\sigma )`$ ($`\eta _{10}=6.00_{0.81}^{+1.10}`$). This value is close to the value implied by the cosmic deuterium abundance in high-redshift Lyman-$`\alpha `$ clouds observed along the line of sight to background quasars  $`\mathrm{\Omega }_bh^2=0.020\pm 0.001(1\sigma )`$ ($`\eta _{10}=5.46\pm 0.27`$). Hence, the newer data imply at most a marginal requirement for a larger baryon density or neutrino degeneracy. Indeed, these new data tighten constraints on the possibility of degenerate cosmological neutrinos. In the present paper we explore the new limits on possible neutrino degeneracy implied by the combined data sets and our BBN constraints.

## V Results

We limit our consideration to flat $`\mathrm{\Omega }_{tot}=\mathrm{\Omega }_M+\mathrm{\Omega }_\mathrm{\Lambda }=1`$ cosmological models with ionization parameter $`\tau =0`$. This is sufficient for our purposes since the likelihood functions so deduced are not expected to be much different if $`\mathrm{\Omega }_{tot}`$ or $`\tau `$ are varied (cf. ). This is because $`\tau `$ and the spectrum tilt $`n`$ are nearly degenerate parameters, i.e. changing one is equivalent to changing the other. Moreover, $`\mathrm{\Omega }_{tot}`$ is generally tightly constrained to be near unity anyway.

There are then nine parameters over which we marginalize. These are $`(\mathrm{\Omega }_bh^2,\xi _{\nu _{\mu ,\tau }},\xi _{\nu _e}\mathrm{\Omega }_\mathrm{\Lambda },h,n,ϵ_i)`$. We utilize a strong Gaussian prior for $`h=0.72\pm 0.08`$ and for the calibration uncertainties $`ϵ_i`$ as listed above. Also, as noted above, we adopt weak top-hat priors when marginalizing over $`\mathrm{\Omega }_bh^2`$ and $`\xi _{\nu _e}`$ designated by the shaded regions of Figure 1 for each value of $`\xi _{\nu _{\mu ,\tau }}`$. In  it has been argued that without some priors on $`\mathrm{\Omega }_M`$ (through flatness, $`h`$, etc.) it is difficult to place bounds on the amount of relativistic matter. Hence, the model constraints adopted here are probably required to break the parameter degeneracy between relativistic and nonrelativistic matter. Ultimately, however, high resolution sky maps such as the Planck mission will be able to determine separately the amounts of relativistic and nonrelativistic matter.

Figure 2 illustrates one of the main results of this study. Shown are contours of constant $`\mathrm{\Delta }\chi ^2`$ in the $`\xi _{\nu _{\mu ,\tau }}`$ vs. $`n`$ plane for three values of the cosmological constant (i.e. $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65,0.75,`$ and $`0.8`$) and for fixed $`h=0.75`$ as noted. For $`\mathrm{\Omega }_\mathrm{\Lambda }0.75`$, a minimum in $`\chi ^2`$ develops for values of $`\xi _{\nu _{\mu ,\tau }}11.5`$. Indeed, for a simple 2 parameter search with fixed values, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.75`$ and $`h=0.75`$, neutrino degeneracy is preferred at the level of more than $`3\sigma `$ over a nondegenerate model. For smaller values of $`\mathrm{\Omega }_\mathrm{\Lambda }`$, this minimum for neutrino-degenerate models becomes even more pronounced.

A second minimum also develops for higher degeneracy ($`\xi _{\nu _{\mu ,\tau }}11.4`$) as noted in . This is due to the large change in particle degrees of freedom for neutrinos which decouple just above the QCD transition. However, the goodness of fit is so poor ($`\mathrm{\Delta }\chi ^2500`$) that it would not be apparent in the contours drawn on Figure 2. Hence, this large-degeneracy solution is definitely ruled out by the current CMB power spectrum.

Figure 3 shows the marginalized likelihood distributions for three of the cosmological parameters ($`\xi _{\nu _{\mu ,\tau }},\mathrm{\Omega }_\mathrm{\Lambda }`$$`n`$) considered here. For the present study the likelihood functions for $`\xi _{\nu _e}`$, $`\mathrm{\Omega }_bh^2`$ and $`\mathrm{\Omega }_M`$ are related to these since $`\xi _{\nu _e}`$ and $`\mathrm{\Omega }_bh^2`$ are functions of $`\xi _{\nu _{\mu ,\tau }}`$ and $`\mathrm{\Omega }_M=1\mathrm{\Omega }_\mathrm{\Lambda }`$. From Figure 3 we deduce optimum values of $`\mathrm{\Omega }_\mathrm{\Lambda }=0.74_{0.11}^{+0.08}`$ and $`n=0.93\pm 0.02`$. A slight preference for finite neutrino degeneracy is evident $`\xi _{\nu _{\mu ,\tau }}=1.0_{1.0(0.5\sigma )}^{+0.8(1\sigma )}`$. This preference, however, is not particularly significant. For now, the data mainly imply $`(2\sigma )`$ upper limits on neutrino degeneracy of $`\xi _{\nu _{\mu ,\tau }}2.1`$. This value implies upper limits of $`\xi _{\nu _e}0.30`$ and $`\Omega _\mathrm{b}h^20.030`$ from Figure 1. For a single large-degeneracy neutrino species, these limits become $`|\xi _{\nu _\mu }|`$ or $`|\xi _{\nu _\tau }|=1.4_{1.4(0.5\sigma )}^{+1.1(1\sigma )}`$ with a $`2\sigma `$ upper limit of $`\xi _{\nu _{\mu ,\tau }}2.5`$.

Our results are slightly more stringent than the results from  who found an equivalent single species upper limit based upon the CMB data alone of $`\xi _{\nu _\mu }`$ or $`\xi _{\nu _\tau }<2.9`$. This is at first surprising given that we have adopted weak top-hat (instead of Gaussian) priors for the BBN constraint. We have traced the main reason for the more stringent upper limits derived here to our adoption of a strong Gaussian prior on $`h`$. A larger neutrino degeneracy is possible if larger values of $`h`$ are permitted. Figure 4 shows contours of constant $`\mathrm{\Delta }\chi ^2`$ in the $`H_0`$ vs. $`n`$ plane for $`\mathrm{\Omega }_\mathrm{\Lambda }=0.75`$ models with $`\xi _{\nu _{\mu ,\tau }}=0`$, $`1.0`$ and $`1.5`$ as labeled. This illustrates the sensitivity of the degenerate solution to the assumed prior for $`h`$. If weaker priors on $`h`$ are adopted, or if new larger values of $`h`$ in the upper range of the present Key-Project uncertainty are ever determined, the neutrino-degenerate models could become strongly preferred over the non-degenerate models. The $`\xi _{\nu _{\mu ,\tau }}=0`$ non-degenerate solution is only the preferred minimum, for all values of $`\mathrm{\Omega }_\mathrm{\Lambda }`$, when $`h0.70`$. This is consistent with the results of .

Figure 5 shows some optimum model power spectra compared with the combined data sets. The solid line shows our optimum degenerate model for which ($`\Omega _\mathrm{b}h^2`$, $`\xi _{\nu _{\mu ,\tau }}`$, $`\xi _{\nu _e}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`h`$, $`n`$) $`=`$ ($`0.021`$, $`1.0`$, $`0.09`$, $`0.74`$, $`0.74`$, $`0.93`$). For this parameter set we obtain a total $`\chi ^2=29.8`$ for 29 degrees of freedom implying a nearly perfect fit. For comparison, the dotted line shows the best non-degenerate ($`\xi _{\nu _{\mu ,\tau }}=\xi _{\nu _e}=0`$) model \[($`\mathrm{\Omega }_bh^2`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`h`$, $`n`$) $`=`$ ($`0.021`$, $`0.62`$, $`0.62`$,$`1.0`$) (dotted line)\] from . For illustration we also show the large-degeneracy minimum \[($`\mathrm{\Omega }_bh^2`$, $`\xi _{\nu _{\mu ,\tau }}`$, $`\xi _{\nu _e}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`h`$, $`n`$) $`=`$ ($`0.052`$, $`11.4`$, $`0.74`$, $`0.45`$, $`0.80`$, $`0.72`$) (dot-dashed line)\].

## VI Conclusions

In neutrino-degenerate models the larger baryon density associated with the observed low deuterium abundance can be more easily accommodated than in non-degenerate models. Moreover, neutrino-degenerate models provide a slightly improved goodness of fit for the latest CMB power spectra from BOOMERANG, DASI, and MAXIMA-1.

Using cosmological models consistent with the constraints from light-element abundances as a function of the neutrino degeneracy parameter $`\xi _{\nu _{\mu ,\tau }}`$, we have shown that a slight maximum in the likelihood function forms for neutrino-degenerate models with $`\xi _{\nu _{\mu ,\tau }}1`$. However, the improvement over the nondegenerate models is only at the level of about 0.5$`\sigma `$. Although this minimum is not particularly statistically significant for the present data set and assumed priors, it could become much more pronounced should larger values of $`h`$ and/or smaller values of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ ever be established near their current $`1\sigma `$ limits.

The present data place $`2\sigma `$ limits for two identical large-degeneracy neutrino species is $`\xi _{\nu _{\mu ,\tau }}2.1`$, which implies $`\xi _{\nu _e}0.30`$. For only one species with large degeneracy, the limit becomes $`|\xi _{\nu _\mu }|`$ or $`|\xi _{\nu _\tau }|2.5`$. This is slightly more restrictive than the limits deduced in other studies.

Finally, we remark that, since neutrino degeneracy is now limited to such small values, the present work has established that the effects of the changing neutrino decoupling temperature with increased degeneracy has little effect. Hence, previous studies which neglected this effect are justified.

###### Acknowledgements.

One of the authors (GJM) wishes to acknowledge the hospitality of the National Astronomical Observatory of Japan where some of this work was done. This work has been supported in part by the Grant-in-Aid for Scientific Research (10640236, 10044103, 11127220, 12047233, 13002001,13640313) of the Ministry of Education, Science, Sports, and Culture of Japan. Work at Univ. Notre Dame supported by DoE Nuclear Theory Grant (DE-FG02-95-ER40934. Work at Univ. Oklahoma supported by NSF CAREER grant AST-0094335.