format: text

text:
# Some Applications for an Euler Angle Parameterization of 𝑆𝑈(𝑁) and 𝑈(𝑁)

## I Introduction

We know that, in general, $`N`$-dimensional density matrices can be written in the form

$$\rho =\frac{1}{N}\left(1\mathrm{l}_N+\sqrt{\frac{N(N1)}{2}}(𝐧𝝀)\right).$$
(1)
This representation is a convenient one since for pure states it yields,

$$𝐧𝐧=1,\text{ and }𝐧𝐧=𝐧,$$
(2)
where the “star” product is defined according to

$$(𝐚𝐛)_k=\sqrt{\frac{N(N1)}{2}}\left(\frac{1}{N2}\right)d_{ijk}a_ib_j.$$
(3)
The way to see this is to just use the following relationship from Gell-Mann (1964); Greiner and Müller (1989)

$$\lambda _i\lambda _j=\frac{2}{N}\delta _{ij}1\mathrm{l}_N+(d_{ijk}+ic_{ijk})\lambda _k,$$
(4)
where the $`c_{ijk}`$ and the $`d_{ijk}`$ are the structure coefficients for the $`SU(N)`$ Lie algebra in question. For $`N`$-dimensional density matrices there exists a simple procedure to calculate the components of $`𝐧`$: beginning with equation (1) we use the trace condition on the elements of the algebra

$$Tr[\lambda _i\lambda _j]=2\delta _{ij}$$
(5)
to generate

$`n_i`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{2N(N1)}}}Tr[(N\rho 1\mathrm{l}_N)\lambda _i]`$ (6)
$`=`$ $`\sqrt{{\displaystyle \frac{N}{2(N1)}}}Tr[\rho \lambda _i]`$
We can use this to represent $`𝐧`$ in terms of other parameterizations such as when $`\rho =U\rho _dU^{}`$222Substitution of $`\rho =U\rho _dU^{}`$ into equation (6) yields
$$n_i=\sqrt{\frac{N}{2(N1)}}Tr[(U\rho _dU^{})\lambda _i]$$
which upon calculation gives $`𝐧`$ as a function of the $`N(N1)`$ group parameters (denoted by $`\alpha _i`$) and the $`N1`$ “rotations” (denoted by $`\theta _i`$) parameterizing $`\rho _d`$. or to evaluate (1) when $`\rho `$ is explicitly given as in the following example.

### I.1 Example Calculation: Bell’s States

We begin by defining the four Bell states Bell (1964):

$$\psi _1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ 0\\ 0\\ 1\end{array}\right),\psi _2=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ 0\\ 0\\ 1\end{array}\right),\psi _3=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ 1\\ 1\\ 0\end{array}\right),\psi _4=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ 1\\ 1\\ 0\end{array}\right),$$
(7)
which can also be represented as the following density matrices

$$\begin{array}{cccc}\psi _1\psi _1^{}\rho _{BS1}=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 1\end{array}\right),& \hfill \psi _2\psi _2^{}\rho _{BS2}=\frac{1}{2}\left(\begin{array}{cccc}1& \hfill 0& 0& \hfill 1\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 1& \hfill 0& 0& \hfill 1\end{array}\right),& & \\ \psi _3\psi _3^{}\rho _{BS3}=\frac{1}{2}\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 1& 1& \hfill 0\\ 0& \hfill 1& 1& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \hfill \psi _4\psi _4^{}\rho _{BS4}=\frac{1}{2}\left(\begin{array}{cccc}0& 0\hfill & 0& \hfill 0\\ 0& 1\hfill & 1& \hfill 0\\ 0& 1\hfill & 1& \hfill 0\\ 0& 0\hfill & 0& \hfill 0\end{array}\right).& & \end{array}$$
(8)
These four density matrices represent the four possible EPR pairs for two qubit systems in $`SU(4)`$ and represent an orthonormal basis for the entire two qubit, pure state, state space (i. e. for the vector space but not for $`\rho `$). By their definition, they are maximally entangled states, i.e. nonfactorizable superpositions of product states, and thus impart non-local correlations between the behavior of the two qubits that make up these four states. For, if these states were factorizable

$$\rho =\rho _A\rho _B,$$
(9)
then the probability of joint detection would also factorize

$$|\rho |^2=|\rho _A|^2|\rho _B|^2,$$
(10)
and thus the measurements would be independent of each other.

Now although the Bell states cannot be decomposed into a set of product states, we can decompose their density matrix representation given in equation (8) into the one given by equation (1) by using equation (6) with $`N=4`$,

$$\rho =\frac{1}{4}(1\mathrm{l}_4+\sqrt{6}𝐧𝝀),$$
(11)
thus yielding

$`\rho _{BS1}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(1\mathrm{l}_4+\sqrt{6}({\displaystyle \frac{1}{\sqrt{6}}}\lambda _3+{\displaystyle \frac{1}{3\sqrt{2}}}\lambda _8+\sqrt{{\displaystyle \frac{2}{3}}}\lambda _9{\displaystyle \frac{1}{3}}\lambda _{15})),`$
$`\rho _{BS2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(1\mathrm{l}_4+\sqrt{6}({\displaystyle \frac{1}{\sqrt{6}}}\lambda _3+{\displaystyle \frac{1}{3\sqrt{2}}}\lambda _8\sqrt{{\displaystyle \frac{2}{3}}}\lambda _9{\displaystyle \frac{1}{3}}\lambda _{15})),`$
$`\rho _{BS3}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(1\mathrm{l}_4+\sqrt{6}({\displaystyle \frac{1}{\sqrt{6}}}\lambda _3+\sqrt{{\displaystyle \frac{2}{3}}}\lambda _6{\displaystyle \frac{1}{3\sqrt{2}}}\lambda _8+{\displaystyle \frac{1}{3}}\lambda _{15})),`$
$`\rho _{BS4}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(1\mathrm{l}_4+\sqrt{6}({\displaystyle \frac{1}{\sqrt{6}}}\lambda _3\sqrt{{\displaystyle \frac{2}{3}}}\lambda _6{\displaystyle \frac{1}{3\sqrt{2}}}\lambda _8+{\displaystyle \frac{1}{3}}\lambda _{15})).`$ (12)
We can see that for the first two Bell states the only non-zero components of $`𝐧`$ are $`n_3`$, $`n_8`$, $`n_9`$ and $`n_{15}`$ whereas for the last two Bell states one switches $`n_9`$ for $`n_6`$ (with its corresponding Lie algebra component). The rest of the components of $`𝐧`$ are zero.

## II General Entangling Operations on Two Qubits

From (Wooters (2001); Zýczkowski et al. (1998); Zýczkowski (1999) and references within) we know that entangled (and entangleable) two qubit density matrices must satisfy333 Note that this is just a necessary condition for a state to either be already entangled, or to possibly be entangled under some group operation. Obviously, pure states satisfy this criterion automatically, but it is the mixed state situation to which this criterion is more often applied. Separable two qubit density matrices have $`Tr[\rho ^2]\frac{1}{3}`$.

$$Tr[\rho ^2]>\frac{1}{3}.$$
(13)
If we start with equation (1), assigning $`N=4`$

$$\rho =\frac{1}{4}\left(1\mathrm{l}_4+\sqrt{6}(𝐧𝝀)\right),$$
(14)
then equation (13) yields

$$Tr[\rho ^2]=Tr\left[\frac{1}{16}\left(1\mathrm{l}_4+\sqrt{6}(𝐧𝝀)\right)^2\right]=\frac{1}{4}(1+3𝐧^2)>\frac{1}{3}$$
(15)
implying that $`𝐧>1/3`$ for entangled states. On the other hand if we begin with $`\rho =U\rho _dU^{}`$ where $`USU(4)`$, we can see that

$$Tr[\rho ^2]=Tr[(U\rho _dU^{})(U\rho _dU^{})]=Tr[U\rho _d^2U^{}]=Tr[\rho _d^2]>\frac{1}{3}$$
(16)
where $`\rho _d`$ is

$$\rho _d=\left(\begin{array}{cccc}\mathrm{sin}^2(\theta _1)\mathrm{sin}^2(\theta _2)\mathrm{sin}^2(\theta _3)& 0& 0& 0\\ 0& \mathrm{cos}^2(\theta _1)\mathrm{sin}^2(\theta _2)\mathrm{sin}^2(\theta _3)& 0& 0\\ 0& 0& \mathrm{cos}^2(\theta _2)\mathrm{sin}^2(\theta _3)& 0\\ 0& 0& 0& \mathrm{cos}^2(\theta _3)\end{array}\right)$$
(17)
from Tilma et al. (2002) and we have exploited the knowledge that $`UU^{}=1\mathrm{l}_4`$. Evaluation of this trace yields the following demands on the ranges of $`\theta _1`$, $`\theta _2`$, and $`\theta _3`$ (here expressed as functions of $`\mathrm{sin}(\theta _i)`$ and $`\mathrm{cos}(\theta _i)`$):

$$Tr[\rho _d^2]=\mathrm{cos}(\theta _3)^4+(\mathrm{cos}(\theta _2)^4+\frac{1}{4}(3+\mathrm{cos}(4\theta _1))\mathrm{sin}(\theta _2)^4)\mathrm{sin}(\theta _3)^4>\frac{1}{3}.$$
(18)
Depending on the ranges of $`\theta _i`$ one could either have a pure or mixed state that would satisfy equation (13).

### II.1 Pure State Entanglement

One can see that for $`\theta _i=\pi /2`$, ($`i=1,2,3`$), equation (18) would be satisfied and thus we would have a $`\rho _d`$ which could be entangled. Therefore using equation (17) with $`\theta _i=\pi /2`$, ($`i=1,2,3`$), we can define a pure state as444Remember that any 4 by 4 matrix with one element along the diagonal equal to unity and the rest equal to zero will be invariant under a $`U(3)`$ subgroup of $`SU(4)`$ and would therefore represent a pure state. The representation given in equation (19) however is ideally suited for the following calculations based on the group parameters remaining after evaluating the coset $`SU(4)/U(3)`$ (see Tilma and Sudarshan (2002b) for more details) since it is those remaining 6 parameters (equally split between three $`\lambda _3`$ phases and three $`\lambda _{(i1)^2+1}`$ ($`i=2,3,4`$) rotations) that are explicitly contained in the pure state volume measure given in Tilma and Sudarshan (2002b) which parameterizes the space $`\text{P}^3`$ via the corresponding Fubini-Study metric (see Byrd (1999, 1998) for the $`SU(3)`$ case) which directly acts upon the $`\rho _{\{1,1\}}`$ element. This pure state representation is also consistent with the generalized $`\rho _d`$, used in Tilma et al. (2002); Tilma and Sudarshan (2002a) given previously.

$$\rho _d=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$
(19)
Using equation (6) we can then calculate the components of $`𝐧`$ in terms of the twelve $`\alpha `$ and three $`\theta `$ parameters when $`\rho `$ is given by $`\rho =U\rho _dU^{}`$ and $`USU(4)`$. From this it would be possible to determine the actual parameter values that would generate the Bell states given in equation (8) by solving the fifteen simultaneous equations implied by the representation of $`𝐧`$ in terms of the Euler parameters. Another, more instructive way is to apply successive unitary operations, $`USU(4)`$, to the two qubit pure state given in equation (19) until one achieves the requisite Bell state.555Because of the $`U(3)`$ invariance of the two qubit pure state given in equation (19), we only have to look at those operations in the coset $`SU(4)/U(3)`$ which we know from Tilma and Sudarshan (2002b) to be represented by
$$e^{i\lambda _3\alpha _1}e^{i\lambda _2\alpha _2}e^{i\lambda _3\alpha _3}e^{i\lambda _5\alpha _4}e^{i\lambda _3\alpha _5}e^{i\lambda _{10}\alpha _6}$$
but it is instructive to see, through the use of the full group $`SU(4)`$, that only those operations generated from the set $`\{\lambda _2,\lambda _5,\lambda _{10}\}`$ yield pure state entanglement.

#### II.1.1 Bell States One and Two

To begin, we first act upon our pure state with the group operations,

$$U=e^{i\lambda _{10}\alpha }\text{ and }U^{}=e^{i\lambda _{10}\alpha },$$
(20)
yielding

$$\rho =U\rho _dU^{}=e^{i\lambda _{10}\alpha }\rho _de^{i\lambda _{10}\alpha },$$
(21)
which in matrix notation is

$`\rho =`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\alpha )& 0& 0& \mathrm{sin}(\alpha )\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ \mathrm{sin}(\alpha )& 0& 0& \mathrm{cos}(\alpha )\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}\mathrm{cos}(\alpha )& 0& 0& \mathrm{sin}(\alpha )\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ \mathrm{sin}(\alpha )& 0& 0& \mathrm{cos}(\alpha )\end{array}\right)`$
$`=`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\alpha )^2& 0& 0& \mathrm{cos}(\alpha )\mathrm{sin}(\alpha )\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \mathrm{cos}(\alpha )\mathrm{sin}(\alpha )& 0& 0& \mathrm{sin}(\alpha )^2\end{array}\right).`$ (22)
Taking the partial transpose of the above density matrix yields

$$\rho ^{pt}=\left(\begin{array}{cccc}\mathrm{cos}(\alpha )^2& 0& 0& 0\\ 0& 0& \mathrm{cos}(\alpha )\mathrm{sin}(\alpha )& 0\\ 0& \mathrm{cos}(\alpha )\mathrm{sin}(\alpha )& 0& 0\\ 0& 0& 0& \mathrm{sin}(\alpha )^2\end{array}\right),$$
(23)
which has an eigenvalue decomposition equal to

$$\{\chi _1,\chi _2,\chi _3,\chi _4\}=\{\mathrm{cos}(\alpha )^2,\mathrm{cos}(\alpha )\mathrm{sin}(\alpha ),\mathrm{cos}(\alpha )\mathrm{sin}(\alpha ),\mathrm{sin}(\alpha )^2\},$$
(24)
and where the constant term in the characteristic polynomial is666The constant term is just the zeroth order coefficient of $`\chi `$ in the characteristic equation.

$$\text{Det}(\rho ^{pt}1\mathrm{l}_4\chi )\mathrm{cos}(\alpha )^4\mathrm{sin}(\alpha )^4.$$
(25)
Recalling that $`\alpha `$ varies from $`0`$ to $`\pi /2`$ we can see that for $`0<\alpha <\frac{\pi }{2}`$ we have an entangled density matrix $`\rho `$. In particular we can see that if $`\alpha =\pi /4`$ we generate the second Bell state

$`\rho (\alpha ={\displaystyle \frac{\pi }{4}})=`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\frac{\pi }{4})^2& 0& 0& \mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& 0& 0& \mathrm{sin}(\frac{\pi }{4})^2\end{array}\right)`$
$`=`$ $`\left(\begin{array}{cccc}\frac{1}{2}& 0& 0& \frac{1}{2}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \frac{1}{2}& 0& 0& \frac{1}{2}\end{array}\right).`$ (26)
Now, in general, if we use the Euler angle parameterization of $`SU(4)`$ given in Tilma et al. (2002) and take the most general $`USU(4)`$ to be given as

$`U=`$ $`U(\alpha _1,0,\alpha _3,0,\alpha _5,\alpha _6,\alpha _7,0,\alpha _9,0,\alpha _{11},0,\alpha _{13},\alpha _{14},\alpha _{15})`$
$`=`$ $`e^{i\lambda _3\alpha _1}e^{i\lambda _3\alpha _3}e^{i\lambda _3\alpha _5}e^{i\lambda _{10}\alpha _6}e^{i\lambda _3\alpha _7}e^{i\lambda _3\alpha _9}e^{i\lambda _3\alpha _{11}}e^{i\lambda _3\alpha _{13}}e^{i\lambda _8\alpha _{14}}e^{i\lambda _{15}\alpha _{15}},`$ (27)
we would then have a density matrix equal to

$`\rho =`$ $`U\rho _dU^{}`$
$`=`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\alpha _6)^2& 0& 0& e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _6)\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _6)& 0& 0& \mathrm{sin}(\alpha _6)^2\end{array}\right),`$ (28)
whose partial transpose is

$$\rho ^{pt}=\left(\begin{array}{cccc}\mathrm{cos}(\alpha _6)^2& 0& 0& 0\\ 0& 0& e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _6)& 0\\ 0& e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _6)& 0& 0\\ 0& 0& 0& \mathrm{sin}(\alpha _6)^2\end{array}\right),$$
(29)
which yields an eigenvalue decomposition and a constant term in the characteristic polynomial equivalent to equations (24) and (25) but that does not generate the second Bell state when $`\alpha _6=\pi /4`$

$`\rho (\alpha _6={\displaystyle \frac{\pi }{4}})=`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\frac{\pi }{4})^2& 0& 0& e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ e^{i(\alpha _1+\alpha _3+\alpha _5)}\mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& 0& 0& \mathrm{sin}(\frac{\pi }{4})^2\end{array}\right)`$
$`=`$ $`\left(\begin{array}{cccc}\frac{1}{2}& 0& 0& \frac{1}{2}e^{i(\alpha _1+\alpha _3+\alpha _5)}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \frac{1}{2}e^{i(\alpha _1+\alpha _3+\alpha _5)}& 0& 0& \frac{1}{2}\end{array}\right).`$ (30)
What we would have, on the other hand, is the following. Define $`\alpha _1+\alpha _3+\alpha _5=\beta `$, where $`0\beta 5\pi `$ since $`\alpha _1`$ runs from $`0`$ to $`\pi `$ and $`\alpha _3`$ and $`\alpha _5`$ each run from $`0`$ to $`2\pi `$.777We are now using the covering ranges for $`SU(4)`$, defined in Tilma et al. (2002). Then

$`e^{\pm i(\alpha _1+\alpha _3+\alpha _5)}=`$ $`e^{\pm i\beta }`$
$`=`$ $`\mathrm{cos}(\beta )\pm i\mathrm{sin}(\beta )`$
$`=`$ $`1\text{ when }\beta =\pi ,3\pi ,5\pi `$
$`=`$ $`1\text{ when }\beta =0,2\pi ,4\pi .`$ (31)
Thus when $`\alpha _1+\alpha _3+\alpha _5=\pi `$, $`3\pi `$, or $`5\pi `$ we get the first Bell state and when $`\alpha _1+\alpha _3+\alpha _5=0`$, $`2\pi `$, or $`4\pi `$ we get the second Bell state. Intermediate values of $`\alpha _1+\alpha _3+\alpha _5`$ can be equated to “intermediate” Bell states; states which have an equivalent density matrix representation as the first and second Bell state, but which are not equal to any type of convex sum of said Bell states.

Now if, instead of equation (20) we were to choose

$$U=e^{i\lambda _2\alpha }\text{ or }U=e^{i\lambda _5\alpha }$$
(32)
and apply it to $`\rho _d`$ (as in equation (21)) we would not generate the other two Bell states, or for that matter any entangled density matrix $`\rho `$, even with the most general U (as in equation (27)). The question now is, what combination of the exponentiation of $`\lambda _2`$ with $`\lambda _5`$, and/or $`\lambda _{10}`$ will entangle the pure state density matrix $`\rho _d`$ yielding the other two Bell states. It is to this question we now proceed.

#### II.1.2 Bell States Three and Four

To begin, we first act upon our pure state with the group operations,

$$U=e^{i\lambda _5\mu }e^{i\lambda _2\nu }\text{ and }U^{}=e^{i\lambda _2\nu }e^{i\lambda _5\mu },$$
(33)
yielding

$$\rho =U\rho _dU^{}=e^{i\lambda _5\mu }e^{i\lambda _2\nu }\rho _de^{i\lambda _2\nu }e^{i\lambda _5\mu },$$
(34)
which in matrix notation is888Recall that since $`\lambda _2`$ and $`\lambda _5`$ do not commute, the other possible group operation $`e^{i\lambda _2\nu }e^{i\lambda _5\mu }\rho _de^{i\lambda _5\mu }e^{i\lambda _2\nu }`$ will not generate this matrix. But, because $`[\lambda _2,\lambda _5]=[\lambda _5,\lambda _2]`$, the subsequent work after this step will be similar for either $`e^{i\lambda _2\nu }e^{i\lambda _5\mu }\rho _de^{i\lambda _5\mu }e^{i\lambda _2\nu }`$ or $`e^{i\lambda _5\mu }e^{i\lambda _2\nu }\rho _de^{i\lambda _2\nu }e^{i\lambda _5\mu }`$.

$`\rho =`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\mu )\mathrm{cos}(\nu )& \mathrm{cos}(\mu )\mathrm{sin}(\nu )& \mathrm{sin}(\mu )& 0\\ \mathrm{sin}(\nu )& \mathrm{cos}(\nu )& 0& 0\\ \mathrm{cos}(\nu )\mathrm{sin}(\mu )& \mathrm{sin}(\mu )\mathrm{sin}(\nu )& \mathrm{cos}(\mu )\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}\mathrm{cos}(\mu )\mathrm{cos}(\nu )& \mathrm{sin}(\nu )& \mathrm{cos}(\nu )\mathrm{sin}(\mu )& 0\\ \mathrm{cos}(\mu )\mathrm{sin}(\nu )& \mathrm{cos}(\nu )& \mathrm{sin}(\mu )\mathrm{sin}(\nu )& 0\\ \mathrm{sin}(\mu )& 0& \mathrm{cos}(\mu )\\ 0& 0& 0& 1\end{array}\right)`$
$`=`$ $`\left(\begin{array}{cccc}\mathrm{cos}(\mu )^2\mathrm{cos}(\nu )^2& \mathrm{cos}(\mu )\mathrm{cos}(\nu )\mathrm{sin}(\nu )& \mathrm{cos}(\mu )\mathrm{cos}(\nu )^2\mathrm{sin}(\mu )& 0\\ \mathrm{cos}(\mu )\mathrm{cos}(\nu )\mathrm{sin}(\nu )& \mathrm{sin}(\nu )^2& \mathrm{cos}(\nu )\mathrm{sin}(\mu )\mathrm{sin}(\nu )& 0\\ \mathrm{cos}(\mu )\mathrm{cos}(\nu )^2\mathrm{sin}(\mu )& \mathrm{cos}(\nu )\mathrm{sin}(\mu )\mathrm{sin}(\nu )& \mathrm{cos}(\nu )^2\mathrm{sin}(\mu )^2& 0\\ 0& 0& 0& 0\end{array}\right).`$ (35)
Now the above density matrix does not look like either of the remaining two Bell states; unless we demand that $`\mu =\pi /2`$. Then we get

$$\rho =\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \mathrm{sin}(\nu )^2& \mathrm{cos}(\nu )\mathrm{sin}(\nu )& 0\\ 0& \mathrm{cos}(\nu )\mathrm{sin}(\nu )& \mathrm{cos}(\nu )^2& 0\\ 0& 0& 0& 0\end{array}\right),$$
(36)
which has the same form as the remaining two Bell states. As before, taking the partial transpose of the above density matrix yields

$$\rho ^{pt}=\left(\begin{array}{cccc}0& 0& 0& \mathrm{cos}(\nu )\mathrm{sin}(\nu )\\ 0& \mathrm{sin}(\nu )^2& 0& 0\\ 0& 0& \mathrm{cos}(\nu )^2& 0\\ \mathrm{cos}(\nu )\mathrm{sin}(\nu )& 0& 0& 0\end{array}\right),$$
(37)
which has an eigenvalue decomposition equal to

$$\{\chi _1,\chi _2,\chi _3,\chi _4\}=\{\mathrm{cos}(\nu )^2,\mathrm{cos}(\nu )\mathrm{sin}(\nu ),\mathrm{cos}(\nu )\mathrm{sin}(\nu ),\mathrm{sin}(\nu )^2\},$$
(38)
and where the constant term in the characteristic polynomial is

$$\text{Det}(\rho ^{pt}1\mathrm{l}_4\chi )\mathrm{cos}(\nu )^4\mathrm{sin}(\nu )^4.$$
(39)
Recalling that $`\nu `$ varies from $`0`$ to $`\pi /2`$ we can see that for $`0<\nu <\frac{\pi }{2}`$ we have an entangled density matrix $`\rho `$. In particular we can see that if $`\nu =\pi /4`$ we generate the third Bell state

$`\rho (\nu ={\displaystyle \frac{\pi }{4}})=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \mathrm{sin}(\frac{\pi }{4})^2& \mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& 0\\ 0& \mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& \mathrm{cos}(\frac{\pi }{4})^2& 0\\ 0& 0& 0& 0\end{array}\right),`$
$`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \frac{1}{2}& \frac{1}{2}& 0\\ 0& \frac{1}{2}& \frac{1}{2}& 0\\ 0& 0& 0& 0\end{array}\right).`$ (40)
Notice, if we had instead used

$$U\rho U^{}=e^{i\lambda _2\psi }e^{i\lambda _5\varphi }\rho _de^{i\lambda _5\varphi }e^{i\lambda _2\psi }$$
(41)
as our initial starting point, we would have instead produced the following density matrix

$$\rho =\left(\begin{array}{cccc}\mathrm{cos}(\varphi )^2\mathrm{cos}(\psi )^2& \mathrm{cos}(\varphi )^2\mathrm{cos}(\psi )\mathrm{sin}(\psi )& \mathrm{cos}(\varphi )\mathrm{cos}(\psi )\mathrm{sin}(\varphi )& 0\\ \mathrm{cos}(\varphi )^2\mathrm{cos}(\psi )\mathrm{sin}(\psi )& \mathrm{cos}(\varphi )^2\mathrm{sin}(\psi )^2& \mathrm{cos}(\varphi )\mathrm{sin}(\varphi )\mathrm{sin}(\psi )& 0\\ \mathrm{cos}(\varphi )\mathrm{cos}(\psi )\mathrm{sin}(\varphi )& \mathrm{cos}(\varphi )\mathrm{sin}(\varphi )\mathrm{sin}(\psi )& \mathrm{sin}(\varphi )^2& 0\\ 0& 0& 0& 0\end{array}\right).$$
(42)
One can see that in order for this matrix to be equivalent to (36), one must demand that $`\psi =\pi /2`$, thus yielding

$$\rho =\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \mathrm{cos}(\varphi )^2& \mathrm{cos}(\varphi )\mathrm{sin}(\varphi )& 0\\ 0& \mathrm{cos}(\varphi )\mathrm{sin}(\varphi )& \mathrm{sin}(\varphi )^2& 0\\ 0& 0& 0& 0\end{array}\right).$$
(43)
For simplicity, we shall generalize this group operation and not the original $`e^{i\lambda _5\mu }e^{i\lambda _2\nu }`$ calculation.

Therefore, in general, if we use the Euler angle parameterization of $`SU(4)`$ given in Tilma et al. (2002) and take $`USU(4)`$ to be given as

$`U=`$ $`U(\alpha _1,{\displaystyle \frac{\pi }{2}},\alpha _3,\alpha _4,\alpha _5,0,\alpha _7,0,\alpha _9,0,\alpha _{11},0,\alpha _{13},\alpha _{14},\alpha _{15})`$
$`=`$ $`e^{i\lambda _3\alpha _1}e^{i\lambda _2\frac{\pi }{2}}e^{i\lambda _3\alpha _3}e^{i\lambda _5\alpha _4}e^{i\lambda _3\alpha _5}e^{i\lambda _3\alpha _7}e^{i\lambda _3\alpha _9}e^{i\lambda _3\alpha _{11}}e^{i\lambda _3\alpha _{13}}e^{i\lambda _8\alpha _{14}}e^{i\lambda _{15}\alpha _{15}},`$ (44)
we would generate the following density matrix

$`\rho =`$ $`U\rho _dU^{}`$
$`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \mathrm{cos}(\alpha _4)^2& e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\alpha _4)\mathrm{sin}(\alpha _4)& 0\\ 0& e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\alpha _4)\mathrm{sin}(\alpha _4)& \mathrm{sin}(\alpha _4)^2& 0\\ 0& 0& 0& 0\end{array}\right),`$ (45)
whose partial transpose is

$$\rho ^{pt}=\left(\begin{array}{cccc}0& 0& 0& e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\alpha _4)\mathrm{sin}(\alpha _4)\\ 0& \mathrm{cos}(\alpha _4)^2& 0& 0\\ 0& 0& \mathrm{sin}(\alpha _4)^2& 0\\ e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\alpha _4)\mathrm{sin}(\alpha _4)& 0& 0& 0\end{array}\right),$$
(46)
which yields an eigenvalue decomposition and a constant term in the characteristic polynomial equivalent to equations (38) and (39) but that does not generate the third Bell state when $`\alpha _4=\pi /4`$

$`\rho (\alpha _4={\displaystyle \frac{\pi }{4}})=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \mathrm{cos}(\frac{\pi }{4})^2& e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& 0\\ 0& e^{i(\alpha _1\alpha _3)}\mathrm{cos}(\frac{\pi }{4})\mathrm{sin}(\frac{\pi }{4})& \mathrm{sin}(\frac{\pi }{4})^2& 0\\ 0& 0& 0& 0\end{array}\right)`$
$`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \frac{1}{2}& \frac{1}{2}e^{i(\alpha _1\alpha _3)}& 0\\ 0& \frac{1}{2}e^{i(\alpha _1\alpha _3)}& \frac{1}{2}& 0\\ 0& 0& 0& 0\end{array}\right).`$ (47)
What we would have, as before with the second Bell state, is the following. Define $`\alpha _1\alpha _3=\gamma `$, where $`2\pi \gamma 0`$ since $`\alpha _1`$ and $`\alpha _3`$ again run from $`0`$ to $`\pi `$ and $`0`$ to $`2\pi `$ respectively. Then

$`e^{\pm i(\alpha _1\alpha _3)}=`$ $`e^{\pm i\gamma }`$
$`=`$ $`\mathrm{cos}(|\gamma |)\pm i\mathrm{sin}(|\gamma |)`$
$`=`$ $`\mathrm{cos}(|\gamma |)i\mathrm{sin}(|\gamma |)`$
$`=`$ $`1\text{ when }|\gamma |=\pi `$
$`=`$ $`1\text{ when }|\gamma |=0,2\pi .`$ (48)
Thus when $`|\alpha _1\alpha _3|=\pi `$ we get the fourth Bell state and when $`|\alpha _1\alpha _3|=0`$ or $`2\pi `$ we get the third Bell state. Intermediate values of $`|\alpha _1\alpha _3|`$ can be thought of as “intermediate” Bell states; states which have an equivalent density matrix representation as the third and fourth Bell state, but which are not equal to any type of convex sum of said Bell states.

#### II.1.3 General Two Qubit Pure State Entanglement

The natural extension of the previous work is to look at the case when we use the Euler angle parameterization of $`SU(4)`$ given in Tilma et al. (2002) and take $`USU(4)`$ to be given as999We only need to look at the first 6 group operations of $`U`$ because they are the ones that “parameterize” the coset $`SU(4)/U(3)=\text{P}^3`$ Tilma and Sudarshan (2002b).

$`U=`$ $`U(\alpha _1,\alpha _2,\alpha _3,\alpha _4,\alpha _5,\alpha _6,\alpha _7,0,\alpha _9,0,\alpha _{11},0,\alpha _{13},\alpha _{14},\alpha _{15})`$
$`=`$ $`e^{i\lambda _3\alpha _1}e^{i\lambda _2\alpha _2}e^{i\lambda _3\alpha _3}e^{i\lambda _5\alpha _4}e^{i\lambda _3\alpha _5}e^{i\lambda _{10}\alpha _6}e^{i\lambda _3\alpha _7}e^{i\lambda _3\alpha _9}e^{i\lambda _3\alpha _{11}}e^{i\lambda _3\alpha _{13}}e^{i\lambda _8\alpha _{14}}e^{i\lambda _{15}\alpha _{15}}.`$ (49)
We would then generate the following density matrix (see $`\rho `$ on next page).

One can see immediately that in order to obtain the general form for Bell states three and four, $`\alpha _2`$ and $`\alpha _6`$ must be set to $`\pi /2`$ and zero respectively. Similarly, in order to obtain the general form for Bell states one and two, $`\alpha _2`$ and $`\alpha _4`$ must be set to zero.

Generally, though, the two sets of eigenvalues of the partial transpose of equation (50)

$`\mathrm{\Psi }_\pm `$ $`=\pm e^{\frac{i}{2}\left(2\alpha _1+\alpha _5\right)}\mathrm{cos}(\alpha _4)\mathrm{cos}(\alpha _6)\times \mathrm{\Delta },`$
$`\mathrm{\Phi }_\pm `$ $`={\displaystyle \frac{1}{2}}\pm {\displaystyle \frac{1}{2}}e^{\frac{i}{2}(2\alpha _1+\alpha _5)}\sqrt{(e^{i(2\alpha _1+\alpha _5)}4\mathrm{cos}(\alpha _4)^2\mathrm{cos}(\alpha _6))}\times \mathrm{\Delta },`$ (51)
where $`\mathrm{\Delta }`$ is

$`\mathrm{\Delta }=`$ $`\left(e^{i\alpha _5}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)+e^{2i\alpha _1}\mathrm{cos}(\alpha _2)\mathrm{sin}(\alpha _6)\right)^{\frac{1}{2}}`$
$`\times \left(e^{2i\alpha _1}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)+e^{i\alpha _5}\mathrm{cos}(\alpha _2)\mathrm{sin}(\alpha _6)\right)^{\frac{1}{2}},`$ (52)
indicate that the phase parameter $`\alpha _3`$ does not contribute. We claim that one only needs three rotations, and one overall phase in order to carry out a general entangling operation on $`\rho _d`$. This can be seen if one expands the constant term from the characteristic polynomial for this situation

$`\text{Det}(\rho ^{pt}1\mathrm{l}_4ϵ)`$ $`e^{2i(2\alpha _1+\alpha _5)}\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^4`$
$`\times \left(e^{i\alpha _5}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)+e^{2i\alpha _1}\mathrm{cos}(\alpha _2)\mathrm{sin}(\alpha _6)\right)^2`$
$`\times \left(e^{2i\alpha _1}\mathrm{cos}(\alpha _6)\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)+e^{i\alpha _5}\mathrm{cos}(\alpha _2)\mathrm{sin}(\alpha _6)\right)^2`$ (53)
and defines $`\eta =2\alpha _1\alpha _5`$ as the cumulative phase

$`ϵ^0`$ $`\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^8\mathrm{sin}(\alpha _2)^4\mathrm{sin}(\alpha _4)^4`$
$`2e^{i\eta }\mathrm{cos}(\alpha _2)\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^7\mathrm{sin}(\alpha _2)^3\mathrm{sin}(\alpha _4)^3\mathrm{sin}(\alpha _6)`$
$`2e^{i\eta }\mathrm{cos}(\alpha _2)\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^7\mathrm{sin}(\alpha _2)^3\mathrm{sin}(\alpha _4)^3\mathrm{sin}(\alpha _6)`$
$`4\mathrm{cos}(\alpha _2)^2\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^6\mathrm{sin}(\alpha _2)^2\mathrm{sin}(\alpha _4)^2\mathrm{sin}(\alpha _6)^2`$
$`e^{2i\eta }\mathrm{cos}(\alpha _2)^2\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^6\mathrm{sin}(\alpha _2)^2\mathrm{sin}(\alpha _4)^2\mathrm{sin}(\alpha _6)^2`$
$`e^{2i\eta }\mathrm{cos}(\alpha _2)^2\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^6\mathrm{sin}(\alpha _2)^2\mathrm{sin}(\alpha _4)^2\mathrm{sin}(\alpha _6)^2`$
$`2e^{i\eta }\mathrm{cos}(\alpha _2)^3\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^5\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)\mathrm{sin}(\alpha _6)^3`$
$`2e^{i\eta }\mathrm{cos}(\alpha _2)^3\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^5\mathrm{sin}(\alpha _2)\mathrm{sin}(\alpha _4)\mathrm{sin}(\alpha _6)^3`$
$`\mathrm{cos}(\alpha _2)^4\mathrm{cos}(\alpha _4)^4\mathrm{cos}(\alpha _6)^4\mathrm{sin}(\alpha _6)^4.`$ (54)
From the previous work it is obvious that only the parameters $`\alpha _{2i}`$ ($`i=1,2,3`$) and $`\eta `$, the overall phase running from $`0`$ to $`2\pi `$, are needed to parameterize all entangling operations that can be done on an initial pure state (given in equation (19)). These three rotations and one phase thus can be used to parameterize the following manifold

$$_{\text{E}_{SU(4)}}=\frac{SU(4)}{U(3)\times U(1)_{SU(2)_\eta }}=\frac{\text{P}^3}{U(1)_{SU(2)_{2\alpha _1\alpha _5}}},$$
(55)
the volume of which can be found using the material from Tilma and Sudarshan (2002b)

$`V_{_{\text{E}_{SU(4)}}}=V\left({\displaystyle \frac{SU(4)}{U(3)\times U(1)_{SU(2)_\eta }}}\right)=V\left({\displaystyle \frac{\text{P}^3}{U(1)_{2\lambda _3}}}\right)`$ $`={\displaystyle \frac{V_{\text{P}^3}}{2V_{U(1)_{\lambda _3}}}}`$
$`={\displaystyle \frac{\pi ^3/6}{2\pi \times 2\pi }}`$
$`={\displaystyle \frac{\pi }{24}}.`$ (56)
What we have been able to derive in equation (56) is the volume of the manifold of all operations on two qubit pure states which produce entanglement. Since these operations act upon a pure state (again, see equation (19)), which is just a point in $`\text{P}^3`$ and thus of measure zero, one may conclude that, up to the volume of a measure zero set, the volume of the manifold of all operations on two qubit pure states which produce entanglement is equivalent to the volume of the set of all entangled two qubit pure states

$`\{\text{Set of All Entangled Two Qubit Pure States}\}`$ $`_{\text{E}_{SU(4)}}`$
$`V_{\text{Set of All Entangled Two Qubit Pure States}}`$ $`V_{_{\text{E}_{SU(4)}}}={\displaystyle \frac{\pi }{24}}.`$ (57)
This volume is less than the numerically estimated value calculated by Zyczkowski et. al. in Zýczkowski et al. (1998)(eq. 30) and referenced in Zýczkowski (1999) for the lower bound of the volume of entangled two qubit states (1-0.863) by approximately 6 one-thousandths.

Unfortunately, whereas our calculation was done with only pure states in mind, the Zyczkowski et. al. calculation was numerically done using a mixed state product measure defined on $`U(N)`$ and randomly chosen density matrices Zýczkowski et al. (1998); Zýczkowski (1999) which satisfied the Peres-Horodecki criterion for separability Peres (1996); Horodecki et al. (1996). Therefore although a tantalizing conclusion, we must concede that until a more general calculation is performed using our mixed state product measure (defined in Tilma and Sudarshan (2002b)) we cannot definitively state that we have calculated the exact volume of the set of all entangled two qubit pure states.

### II.2 Mixed State Entanglement

Using common mathematical software, one can see that for $`\theta _i\pi /2`$, ($`i=1,2,3`$) equation (18) could still be satisfied if

$$\mathrm{sin}(\theta _1)=\sqrt{\frac{1}{2}+\frac{\sqrt{32789757}}{12482}}\text{ and }$$

$$1>\mathrm{sin}(\theta _2)>\frac{1}{\sqrt{2}}\sqrt{\frac{1+\sqrt{\mathrm{sin}(\theta _1)^2\mathrm{sin}(\theta _1)^4}}{1\mathrm{sin}(\theta _1)^2+\mathrm{sin}(\theta _1)^4}}\text{ and }1>\mathrm{sin}(\theta _3)>\frac{\sqrt{3}}{2}$$
(58)
or

$$\mathrm{sin}(\theta _1)=\sqrt{\frac{1}{2}+\frac{\sqrt{\frac{785323439}{3}}}{37446}}\text{ and }$$

$$1>\mathrm{sin}(\theta _2)>\frac{1}{3\sqrt{2}}\sqrt{\frac{9+\sqrt{3}\sqrt{1+28\mathrm{sin}(\theta _1)^228\mathrm{sin}(\theta _1)^4}}{1\mathrm{sin}(\theta _1)^2+\mathrm{sin}(\theta _1)^4}}\text{ and }1>\mathrm{sin}(\theta _3)>\frac{\sqrt{3}}{2}$$
(59)
or

$$\frac{1}{\sqrt{2}}<\mathrm{sin}(\theta _1)<\sqrt{\frac{1}{2}+\frac{\sqrt{\frac{785323439}{3}}}{37446}}\text{ and }$$

$$\frac{79}{100}<\mathrm{sin}(\theta _2)\frac{1}{3\sqrt{2}}\sqrt{\frac{9+\sqrt{3}\sqrt{1+28\mathrm{sin}(\theta _1)^228\mathrm{sin}(\theta _1)^4}}{1\mathrm{sin}(\theta _1)^2+\mathrm{sin}(\theta _1)^4}}\text{ and }$$

$$1>\mathrm{sin}(\theta _3)>\frac{1}{\sqrt{6}}\sqrt{\frac{3+\sqrt{3}\sqrt{1+4\mathrm{sin}(\theta _2)^24\mathrm{sin}(\theta _2)^4+4\mathrm{sin}(\theta _1)^2\mathrm{sin}(\theta _2)^44\mathrm{sin}(\theta _1)^4\mathrm{sin}(\theta _2)^4}}{1\mathrm{sin}(\theta _2)^2+\mathrm{sin}(\theta _2)^4\mathrm{sin}(\theta _1)^2\mathrm{sin}(\theta _2)^4+\mathrm{sin}(\theta _1)^4\mathrm{sin}(\theta _2)^4}}$$
(60)
or

$$1>\mathrm{sin}(\theta _1)>\sqrt{\frac{1}{2}+\frac{\sqrt{32789757}}{12482}}\text{ and }1>\mathrm{sin}(\theta _2)>\frac{79}{100}\text{ and }1>\mathrm{sin}(\theta _3)>\frac{\sqrt{3}}{2}$$
(61)
or

$$\sqrt{\frac{1}{2}+\frac{\sqrt{32789757}}{12482}}>\mathrm{sin}(\theta _1)>\sqrt{\frac{1}{2}+\frac{\sqrt{\frac{785323439}{3}}}{37446}}\text{ and }1>\mathrm{sin}(\theta _2)>\frac{79}{100}\text{ and }1>\mathrm{sin}(\theta _3)>\frac{\sqrt{3}}{2}$$
(62)
or

$$\frac{1}{\sqrt{2}}<\mathrm{sin}(\theta _1)<\sqrt{\frac{1}{2}+\frac{\sqrt{\frac{785323439}{3}}}{37446}}\text{ and }$$

$$1>\mathrm{sin}(\theta _2)>\frac{1}{3\sqrt{2}}\sqrt{\frac{9+\sqrt{3}\sqrt{1+28\mathrm{sin}(\theta _1)^228\mathrm{sin}(\theta _1)^4}}{1\mathrm{sin}(\theta _1)^2+\mathrm{sin}(\theta _1)^4}}\text{ and }1>\mathrm{sin}(\theta _3)>\frac{\sqrt{3}}{2}.$$
(63)
Therefore we could generate a $`\rho _d`$ which could be entangled but would no longer be a pure state. Unfortunately, the fact that $`\rho _d`$ is no longer a pure state also means that we would have to look at the most general $`USU(4)`$ acting on $`\rho _d`$ in order to determine which successive unitary operations $`U\rho _dU^{}`$ would produce entanglement. This is a rather lengthy and complicated calculation and is beyond the scope of this section.

What we can do though is make an educated guess as to the volume of entangled two qubit mixed states by noticing that it is a product of the volume of the 3-dimensional symplex of eigenvalues of $`\rho _d`$ (with appropriate ranges) and the volume of the flag manifold $`SU(4)/U(1)_{SU(2)}\times U(1)_{SU(3)}\times U(1)_{SU(4)}`$.101010Explained in detail in Tilma and Sudarshan (2002b). These volumes can be calculated by using the mixed state product measure given in Tilma and Sudarshan (2002b) with the necessary ranges for the eigenvalues (given above) and the covering ranges for $`SU(4)`$ (given in Tilma et al. (2002)). Using this material we arrive at

$$dV_{E_{ms}}=\alpha _s\mathrm{\Lambda }_1^{s1}\mathrm{\Lambda }_2^{s1}\mathrm{\Lambda }_3^{s1}(1\underset{i=1}{\overset{3}{}}\mathrm{\Lambda }_i)^{s1}d\mathrm{\Lambda }_1\mathrm{}d\mathrm{\Lambda }_3\times d\left(\frac{SU(4)}{U(1)_{SU(2)}\times U(1)_{SU(3)}\times U(1)_{SU(4)}}\right)d\alpha _{12}\mathrm{}d\alpha _1.$$
(64)
A general (and rather naive) evaluation of this measure for our situation yields

$$V_{E_{ms}}=\alpha _s\frac{\left(a_{L}^{}{}_{}{}^{s}a_{U}^{}{}_{}{}^{s}\right)\left(b_{L}^{}{}_{}{}^{s}b_{U}^{}{}_{}{}^{s}\right)\left(c_{L}^{}{}_{}{}^{s}c_{U}^{}{}_{}{}^{s}\right)\left(d_{L}^{}{}_{}{}^{s}d_{U}^{}{}_{}{}^{s}\right)}{s^4}\times \frac{\pi ^6}{12}=\omega \frac{\pi ^6}{12},$$
(65)
where $`\{a_U,a_L\}`$,$`\{b_U,b_L\}`$,$`\{c_U,c_L\}`$ and $`\{d_U,d_L\}`$ are the squared values of the above maximal and minimal ranges ($`\{d_U,d_L\}`$ comes from $`\mathrm{\Lambda }_41_{i=1}^3\mathrm{\Lambda }_i`$). Since the symplex measure is assumed to be non-zero, and using the work contained in Tilma and Sudarshan (2002b) we can hypothesize that $`\omega `$ has the following bounds (dependent on the value of $`s`$ and recalling that $`\alpha _s>0`$ and $`s>0`$)

$$0<\left(a_{L}^{}{}_{}{}^{s}a_{U}^{}{}_{}{}^{s}\right)\left(b_{L}^{}{}_{}{}^{s}b_{U}^{}{}_{}{}^{s}\right)\left(c_{L}^{}{}_{}{}^{s}c_{U}^{}{}_{}{}^{s}\right)\left(d_{L}^{}{}_{}{}^{s}d_{U}^{}{}_{}{}^{s}\right)<4^{4s}(1+4^s).$$
(66)
Notice that the right side of the above inequality approaches 0 when both $`s\mathrm{}`$ and $`s0`$, thus we can conclude that the numerator of (65) will be $`<1`$ and therefore, the value of $`\omega `$ (be it either $`>1`$ or $`<1`$) will be completely dependent on the explicit choice of the value of $`s`$.

The important point to recognize here is not the symplex calculations but rather the flag manifold volume. Through the Euler parameterization of $`SU(N)`$ and $`U(N)`$ given in Tilma and Sudarshan (2002a, b) we have been able to generate the appropriate representation of the “truncated” Haar measure which is crucial to any mixed state volume calculation. It is this factor which is not “user dependent”; i. e. dependent on the initial distribution chosen for the ($`N1`$)-dimensional symplex, and therefore not completely subject to disagreements between researchers studying entanglement.111111Although disagreements in numerical values are found, they are mostly due to variations in the ranges of the $`N(N1)`$ parameters which define the measure (see for example Slater (2002) and references within).

## III Qubit/Qutrit Entanglement

### III.1 Pure State Entanglement

By following the same procedure as was done in the two qubit case, we can derive the manifold of operations that produce entanglement of an initial qubit/qutrit pure state

$$\rho _d=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$$
(67)
via a $`USU(6)`$ given in Tilma and Sudarshan (2002a)

$`U=`$ $`U(\alpha _1,\alpha _2,\alpha _3,\alpha _4,\alpha _5,\alpha _6,\alpha _7,\alpha _8,\alpha _9,\alpha _{10},\alpha _{11},0,\alpha _{13},0,\alpha _{15},0,\alpha _{17},`$
$`\times 0,\alpha _{19},0,\alpha _{21},0,\alpha _{23},0,\alpha _{25},0,\alpha _{27},0,\alpha _{29},0,\alpha _{31},\alpha _{32},\alpha _{33},\alpha _{34},\alpha _{35})`$
$`U=`$ $`e^{i\lambda _3\alpha _1}e^{i\lambda _2\alpha _2}e^{i\lambda _3\alpha _3}e^{i\lambda _5\alpha _4}e^{i\lambda _3\alpha _5}e^{i\lambda _{10}\alpha _6}e^{i\lambda _3\alpha _7}e^{i\lambda _{17}\alpha _8}e^{i\lambda _3\alpha _9}e^{i\lambda _{26}\alpha _{10}}`$
$`\times e^{i\lambda _3\alpha _{11}}e^{i\lambda _3\alpha _{13}}e^{i\lambda _3\alpha _{15}}e^{i\lambda _3\alpha _{17}}e^{i\lambda _3\alpha _{19}}e^{i\lambda _3\alpha _{21}}e^{i\lambda _3\alpha _{23}}`$
$`\times e^{i\lambda _3\alpha _{25}}e^{i\lambda _3\alpha _{27}}e^{i\lambda _3\alpha _{29}}e^{i\lambda _3\alpha _{31}}e^{i\lambda _{15}\alpha _{33}}e^{i\lambda _{24}\alpha _{34}}e^{i\lambda _{35}\alpha _{35}},`$ (68)
by taking the partial transpose of $`U\rho _dU^{}`$ and evaluating the corresponding eigenvalues through the Peres-Horodecki criterion.121212We are forced to generate the six eigenvalue equations rather than evaluate the constant term from the characteristic polynomial because in this case we have two eigenvalues equal to zero thus negating the constant term’s effectiveness. Doing this work, which did not yield well to simplification, generates the following hypothesized manifold

$$_{\text{E}_{SU(6)}}=\frac{SU(6)}{U(5)\times U(1)_{SU(2)_\kappa }}=\frac{\text{P}^5}{U(1)_{SU(2)_\kappa }},$$
(69)
where $`\kappa 2(2\alpha _1+\alpha _3+\alpha _7+\alpha _9)`$ is the cumulative phase.131313After taking into account the degeneracy in the eigenvalues, full simplification of the resulting 4th order characteristic polynomial was not possible without making certain numerical assumptions. Thus, the actual representation of $`\kappa `$ as a function of $`\alpha _1,\alpha _3,\alpha _7`$ and $`\alpha _9`$ was not possible; thus the equivalence. The volume of the manifold is

$`V_{_{\text{E}_{SU(6)}}}=V\left({\displaystyle \frac{SU(6)}{U(5)\times U(1)_{SU(2)_\kappa }}}\right)`$ $`={\displaystyle \frac{V_{\text{P}^5}}{10V_{U(1)_{\lambda _3}}}}`$
$`={\displaystyle \frac{\pi ^5}{5!20\pi ^4}}`$
$`={\displaystyle \frac{\pi }{2400}}.`$ (70)
If we make the same argument for the volume of this manifold to be equivalent to the volume of entangled qubit/qutrit pure states, as we did in the two qubit case, then, in this case, our values are within the ranges specified in Zýczkowski et al. (1998); Zýczkowski (1999) but since we do not include the possibility of bound entangled states (states which are entangled, but which have positive partial transposes) we must concede that our value given in equation (70) is probably too small.

### III.2 Mixed State Entanglement

Again, as in the two qubit case, one can see that for a general $`\rho _d`$ for a qubit/qutrit system Tilma and Sudarshan (2002a)

$`\rho _d=\text{diag}\{`$ $`\mathrm{sin}(\theta _1)^2\mathrm{sin}(\theta _2)^2\mathrm{sin}(\theta _3)^2\mathrm{sin}(\theta _4)^2\mathrm{sin}(\theta _5)^2,\mathrm{cos}(\theta _1)^2\mathrm{sin}(\theta _2)^2\mathrm{sin}(\theta _3)^2\mathrm{sin}(\theta _4)^2\mathrm{sin}(\theta _5)^2,`$
$`\mathrm{cos}(\theta _2)^2\mathrm{sin}(\theta _3)^2\mathrm{sin}(\theta _4)^2\mathrm{sin}(\theta _5)^2,\mathrm{cos}(\theta _3)^2\mathrm{sin}(\theta _4)^2\mathrm{sin}(\theta _5)^2,\mathrm{cos}(\theta _4)^2\mathrm{sin}(\theta _5)^2,\mathrm{cos}(\theta _5)^2\}`$ (71)
when $`\theta _i\pi /2`$, ($`i=1,\mathrm{},5`$) the generalization of equation (13) for the qubit/qutrit case Zýczkowski et al. (1998); Zýczkowski (1999)

$$Tr[\rho ^2]>\frac{1}{5}$$
(72)
where

$`Tr[\rho ^2]=`$ $`Tr[(U\rho _dU^{})^2]=Tr[\rho _d^2]`$
$`=`$ $`\mathrm{\hspace{0.17em}1}2\mathrm{sin}(\theta _5)^2(1+\mathrm{sin}(\theta _5)^2\mathrm{sin}(\theta _4)^2\mathrm{sin}(\theta _5)^2+\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2`$
$`\mathrm{sin}(\theta _3)^2\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2+\mathrm{sin}(\theta _3)^4\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2\mathrm{sin}(\theta _2)^2\mathrm{sin}(\theta _3)^4\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2`$
$`+\mathrm{sin}(\theta _2)^4\mathrm{sin}(\theta _3)^4\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2\mathrm{sin}(\theta _1)^2\mathrm{sin}(\theta _2)^4\mathrm{sin}(\theta _3)^4\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2`$
$`+\mathrm{sin}(\theta _1)^4\mathrm{sin}(\theta _2)^4\mathrm{sin}(\theta _3)^4\mathrm{sin}(\theta _4)^4\mathrm{sin}(\theta _5)^2)`$ (73)
could still be satisfied. In this case then we would have for the entangled mixed state product measure (under appropriate ranges for $`\mathrm{\Lambda }_i`$ and $`\alpha _i`$)

$`dV_{E_{ms}}=`$ $`\alpha _s\mathrm{\Lambda }_1^{s1}\mathrm{}\mathrm{\Lambda }_5^{s1}(1{\displaystyle \underset{i=1}{\overset{5}{}}}\mathrm{\Lambda }_i)^{s1}d\mathrm{\Lambda }_1\mathrm{}d\mathrm{\Lambda }_5`$
$`\times d\left({\displaystyle \frac{SU(6)}{U(1)_{SU(2)}\times U(1)_{SU(3)}\times U(1)_{SU(4)}\times U(1)_{SU(5)}\times U(1)_{SU(6)}}}\right)d\alpha _{30}\mathrm{}d\alpha _1.`$ (74)
Another general (and again, rather naive) evaluation of this measure for our situation yields

$`V_{E_{ms}}`$ $`=\alpha _s{\displaystyle \frac{\left(a_{L}^{}{}_{}{}^{s}a_{U}^{}{}_{}{}^{s}\right)\left(b_{L}^{}{}_{}{}^{s}b_{U}^{}{}_{}{}^{s}\right)\left(c_{L}^{}{}_{}{}^{s}c_{U}^{}{}_{}{}^{s}\right)\left(d_{L}^{}{}_{}{}^{s}d_{U}^{}{}_{}{}^{s}\right)\left(e_{L}^{}{}_{}{}^{s}e_{U}^{}{}_{}{}^{s}\right)\left(f_{L}^{}{}_{}{}^{s}f_{U}^{}{}_{}{}^{s}\right)}{s^6}}\times {\displaystyle \frac{\pi ^{15}}{34560}}`$
$`=\omega {\displaystyle \frac{\pi ^{15}}{34560}},`$ (75)
where $`\{a_U,a_L\}`$ through $`\{f_U,f_L\}`$ are the squared values of the maximal and minimal ranges of $`\mathrm{\Lambda }_i`$ that satisfy (72). As before, the symplex measure must be assumed to be non-zero, therefore using the work contained in Tilma and Sudarshan (2002b) we can hypothesize that $`\omega `$ has the following bounds (dependent on the value of $`s`$ and recalling that $`\alpha _s>0`$ and $`s>0`$)

$$0<\left(a_{L}^{}{}_{}{}^{s}a_{U}^{}{}_{}{}^{s}\right)\left(b_{L}^{}{}_{}{}^{s}b_{U}^{}{}_{}{}^{s}\right)\left(c_{L}^{}{}_{}{}^{s}c_{U}^{}{}_{}{}^{s}\right)\left(d_{L}^{}{}_{}{}^{s}d_{U}^{}{}_{}{}^{s}\right)\left(e_{L}^{}{}_{}{}^{s}e_{U}^{}{}_{}{}^{s}\right)\left(f_{L}^{}{}_{}{}^{s}f_{U}^{}{}_{}{}^{s}\right)<6^{6s}(1+6^s).$$
(76)
Again we notice that the right side of the above inequality approaches 0 when both $`s\mathrm{}`$ and $`s0`$, we can again conclude that the numerator of (III.2) will be $`<1`$ and therefore, the value of $`\omega `$ (be it either $`>1`$ or $`<1`$) will be completely dependent on the explicit choice of the value of $`s`$. Also, as before in the two qubit case, the important point to recognize here is not the symplex calculations but rather the flag manifold volume.

## IV Conclusions

In this paper we have applied our $`SU(N)`$ and $`U(N)`$ parameterizations to the two qubit and qubit/qutrit system in order to explicit calculate the manifold of operations which entangle two qubit and qubit/qutrit pure states. We have also been able to give the volume of this manifold, as well as the hypothesized volume for the set of all entangled two qubit and qubit/qutrit pure and mixed states. In the pure state case, the values were within the ranges given by Zýczkowski et al. (1998); Zýczkowski (1999) but in the qubit/qutrit case, because we did not take into account the possibility of bound entangled states (which do not appear in the two qubit case) our volume is most likely smaller than the actual volume for the set of all entangled qubit/qutrit pure states.

Work is continuing on the mixed state situation; explicitly in calculating the volume of the mixed state manifold without having to know the exact probability distribution on the ($`N1`$)-dimensional symplex. Extensions of the pure state work to two qutrit systems is also ongoing.

## Acknowledgments

We would like to thank Dr. M. Byrd for his editorial help on the various representational conventions for density matrices as well as Anil Shaji for invaluable assistance in calculating the various manifold volumes given here.

## Appendix A $`SU(4)`$ Lie Algebra

From Tilma et al. (2002) we know that the Gell-Mann type basis for the Lie algebra of $`SU(4)`$ is given by the following set of matrices Greiner and Müller (1989):

$$\begin{array}{cccc}\lambda _1=\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),& \hfill \lambda _2=\left(\begin{array}{cccc}0& \hfill i& 0& \hfill 0\\ i& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \lambda _3=\left(\begin{array}{cccc}1& \hfill 0& 0& \hfill 0\\ 0& \hfill 1& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \\ \lambda _4=\left(\begin{array}{cccc}0& 0\hfill & 1& \hfill 0\\ 0& 0\hfill & 0& \hfill 0\\ 1& 0\hfill & 0& \hfill 0\\ 0& 0\hfill & 0& \hfill 0\end{array}\right),& \hfill \lambda _5=\left(\begin{array}{cccc}0& \hfill 0& i& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ i& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \lambda _6=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 1& \hfill 0\\ 0& \hfill 1& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \\ \lambda _7=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& i& \hfill 0\\ 0& \hfill i& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \hfill \lambda _8=\frac{1}{\sqrt{3}}\left(\begin{array}{cccc}1& \hfill 0& 0& \hfill 0\\ 0& \hfill 1& 0& \hfill 0\\ 0& \hfill 0& 2& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\end{array}\right),& \lambda _9=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 1\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 1& \hfill 0& 0& \hfill 0\end{array}\right),& \\ \lambda _{10}=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill i\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ i& \hfill 0& 0& \hfill 0\end{array}\right),& \hfill \lambda _{11}=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 1\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 1& 0& \hfill 0\end{array}\right),& \lambda _{12}=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill i\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill i& 0& \hfill 0\end{array}\right),& \\ \lambda _{13}=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 1\\ 0& \hfill 0& 1& \hfill 0\end{array}\right),& \hfill \lambda _{14}=\left(\begin{array}{cccc}0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill 0\\ 0& \hfill 0& 0& \hfill i\\ 0& \hfill 0& i& \hfill 0\end{array}\right),& \lambda _{15}=\frac{1}{\sqrt{6}}\left(\begin{array}{cccc}1& \hfill 0& 0& \hfill 0\\ 0& \hfill 1& 0& \hfill 0\\ 0& \hfill 0& 1& \hfill 0\\ 0& \hfill 0& 0& \hfill 3\end{array}\right).& \end{array}$$
(77)
Using these matrices one can then generate the various group operations given in section II. Similarly, in Tilma and Sudarshan (2002a) one can see how to construct the $`N^21`$ elements of the $`SU(N)`$ Lie algebra necessary for general $`SU(N)`$ group operations.


text:
# A Wave-Mechanical Approach to Cosmic Structure Formation

## 1 Introduction

The local Universe displays a rich hierarchical pattern of galaxy clustering that encompasses a large range of length scales, culminating in rich clusters and superclusters; the early Universe, however, was almost smooth, with only slight ripples seen in the cosmic microwave background radiation (e.g. Smoot et al. 1992). Models of the evolution of structure link these observations by appealing to the effect of gravitational instability, e.g. Lahav et al. (2002). Low-amplitude primordial perturbations of wave-like character within a largely homogeneous universe become amplified by their own self-gravity, first linearly which then, as their amplitudes build up, results in non-linear collapse. The linear theory of perturbation growth for cosmological density fluctuations is well-established, and is founded on tried-and-tested representation of matter on cosmological scales as a cold fluid (e.g. Peebles 1980). The non-linear regime is much more complicated and generally not amenable to analytic solution. Instead, numerical $`N`$–body simulations using particles rather than waves to trace cosmic structure, have led the way towards an understanding of strongly developed clustering. Analytic approaches based on particles, chiefly the Zel’dovich (1970) approximation, have also proved useful. Although these agree to an extent with fully numerical calculations (Coles, Melott & Shandarin 1993), they break down when particle trajectories cross to form singularities known as caustics.

In this paper we adopt an approach to cosmic structure that has advantages over both fluid and particle descriptions. Following Widrow & Kaiser (1993) and Coles (2002), we construct a formalism in which the dynamics of gravitational is couched in the language of quantum mechanics and governed by a Schrödinger wave equation. The advantages of this approach will be made apparent as we explore its behaviour in simple examples which can be solved exactly and, within the Schödinger approach, to various levels of approximation. Our aim is to construct a useful tool for evolving density fluctuations into the mildly non-linear regime, beyond the breakdown of linear perturbation theory.

The importance of the search for analytic understanding of structure formation is often overlooked in the light of the tremendous advances that have been made recently in computational astrophysics. One motivation for a wider range of theoretical techniques is practical. High-resolution $`N`$-body experiments are CPU-intensive and it is difficult to find sufficient resources to use them to explore a large range of choices of initial fluctuations, different cosmological parameters, and so on. A neat analytic approximation can do this job more efficiently, at least to the point of narrowing down the field of contenders for fuller study. The other motivation is methodological: the complex “cosmic web” of large-scale structure requires explanation rather than reproduction, and explanation rests on using analytic approaches to identify the important physics as much as possible.

The layout of this paper is as follows. In Section 2 we run through the standard basics of gravitational instability in an expanding universe, including perturbation theory and the Zel’dovich approximation. We then, in Section 3, present the details of the Schrödinger approach and derive a spherically-symmetric solution. In Section 4 we apply the method to examples of one-dimensional collapse showing how various approximations within the wave-mechanical framework behave when compared to the real data. We sum up in Section 5.

## 2 Cosmological Structure Formation

### 2.1 The Background Cosmology

Space-times consistent with the Cosmological Principle can be described by the Robertson–Walker metric

$$\mathrm{d}s^2=c^2\mathrm{d}t^2a^2(t)\left(\frac{\mathrm{d}r^2}{1\kappa r^2}+r^2\mathrm{d}\theta ^2+r^2\mathrm{sin}^2\theta \mathrm{d}\varphi ^2\right),$$
(1)
where $`\kappa `$ is the spatial curvature, scaled so as to take the values $`0`$ or $`\pm 1`$. The case $`\kappa =0`$ represents flat space sections, and the other two cases are space sections of constant positive or negative curvature, respectively. The time coordinate $`t`$ is cosmological proper time and $`a(t)`$ is the cosmic scale factor. The dynamics of a Friedmann-Robertson-Walker universe are determined by the Einstein gravitational field equations which can be written

$`3\left({\displaystyle \frac{\dot{a}}{a}}\right)^2`$ $`=`$ $`8\pi G\rho {\displaystyle \frac{3\kappa c^2}{a^2}}+\mathrm{\Lambda }c^2,`$ (2)
$`{\displaystyle \frac{\ddot{a}}{a}}`$ $`=`$ $`{\displaystyle \frac{4\pi G}{3}}\left(\rho +3{\displaystyle \frac{p}{c^2}}\right)+{\displaystyle \frac{\mathrm{\Lambda }c^2}{3}},`$ (3)
$`\dot{\rho }`$ $`=`$ $`3{\displaystyle \frac{\dot{a}}{a}}\left(\rho +{\displaystyle \frac{p}{c^2}}\right).`$ (4)
These equations determine the time evolution of the cosmic scale factor $`a(t)`$ (the dots denote derivatives with respect to cosmological proper time $`t`$) and therefore describe the global expansion or contraction of the universe. The behaviour of these models can further be parametrized in terms of the Hubble parameter $`H=(\dot{a}/a)`$ and the density parameter $`\mathrm{\Omega }=8\pi G\rho /3H^2`$, a suffix $`0`$ representing the value of these quantities at the present epoch when $`t=t_0`$.

### 2.2 Fluid Treatment

In the standard treatment of the Jeans Instability one begins with the dynamical equations governing the behaviour of a self-gravitating fluid. These are the Euler equation

$$\frac{(𝐯)}{t}+(𝐯)𝐯+\frac{1}{\rho }p+\varphi =0;$$
(5)
the continuity equation

$$\frac{\rho }{t}+(\rho 𝐯)=0,$$
(6)
expressing the conservation of matter; and the Poisson equation

$$^2\varphi =4\pi G\rho ,$$
(7)
describing Newtonian gravity. If the length scale of the perturbations is smaller than the effective cosmological horizon $`d_H=c/H`$, a Newtonian treatment of cosmic structure formation is still expected to be valid in expanding world models. In an expanding cosmological background, the Newtonian equations governing the motion of gravitating particles can be written in terms of

$$𝐱𝐫/a(t)$$
(8)
(the comoving spatial coordinate, which is fixed for observers moving with the Hubble expansion),

$$𝐯\dot{𝐫}H𝐫=a\dot{𝐱}$$
(9)
(the peculiar velocity field, representing departures of the matter motion from pure Hubble expansion), $`\rho (𝐱,t)`$ (the matter density), and $`\varphi (𝐱,t)`$ (the peculiar Newtonian gravitational potential, i.e. the fluctuations in potential with respect to the homogeneous background) determined by the Poisson equation in the form

$$_{𝐱}^{}{}_{}{}^{2}\varphi =4\pi Ga^2(\rho \rho _0)=4\pi Ga^2\rho _0\delta .$$
(10)
In this equation and the following the suffix on $`_x`$ indicates derivatives with respect to the new comoving coordinates. Here $`\rho _0`$ is the mean background density, and

$$\delta \frac{\rho \rho _0}{\rho _0}$$
(11)
is the density contrast. Using these variables the Euler equation becomes

$$\frac{(a𝐯)}{t}+(𝐯_𝐱)𝐯=\frac{1}{\rho }_𝐱p_𝐱\varphi .$$
(12)
The first term on the right-hand-side of equation (12) arises from pressure gradients, and is neglected in dust-dominated cosmologies. Pressure effects may nevertheless be important in the the (collisional) baryonic component of the mass distribution when nonlinear structures eventually form. The second term on the right-hand side of equation (12) is the peculiar gravitational force, which can be written in terms of $`𝐠=_𝐱\varphi /a`$, the peculiar gravitational acceleration of the fluid element. If the velocity flow is irrotational, $`𝐯`$ can be rewritten in terms of a velocity potential $`\varphi _v`$:

$$𝐯=_𝐱\varphi _v/a.$$
(13)
This is expected to be the case in the cosmological setting because (a) there are no sources of vorticity in these equations and (b) vortical perturbation modes decay with the expansion. We also have a revised form of the continuity equation:

$$\frac{\rho }{t}+3H\rho +\frac{1}{a}_𝐱(\rho 𝐯)=0.$$
(14)
### 2.3 Linear Perturbation Theory

The procedure for handling perturbations is to linearise the Euler, continuity and Poisson equations by perturbing physical quantities defined as functions of Eulerian coordinates, i.e. relative to an unperturbed coordinate system. Expanding $`\rho `$, $`𝐯`$ and $`\varphi `$ perturbatively and keeping only the first-order terms in equations (12) and (14) gives the linearised continuity equation:

$$\frac{\delta }{t}=\frac{1}{a}_𝐱𝐯,$$
(15)
which can be inverted, with a suitable choice of boundary conditions, to yield

$$\delta =\frac{1}{aHf}\left(_𝐱𝐯\right).$$
(16)
The function $`f\mathrm{\Omega }_0^{0.6}`$; this is simply a fitting formula to the full solution (Peebles 1980). The linearised Euler and Poisson equations are

$$\frac{𝐯}{t}+\frac{\dot{a}}{a}𝐯=\frac{1}{\rho a}_𝐱p\frac{1}{a}_𝐱\varphi ,$$
(17)
$$_{𝐱}^{}{}_{}{}^{2}\varphi =4\pi Ga^2\rho _0\delta ;$$
(18)
$`|𝐯|,|\varphi |,|\delta |1`$ in equations (16), (17) & (18). From these equations, and if one ignores pressure forces, it is easy to obtain an equation for the evolution of $`\delta `$:

$$\ddot{\delta }+2H\dot{\delta }\frac{3}{2}\mathrm{\Omega }H^2\delta =0.$$
(19)
For a spatially flat universe dominated by pressureless matter, $`\rho _0(t)=1/6\pi Gt^2`$ and equation (19) admits two linearly independent power law solutions $`\delta (𝐱,t)=D_\pm (t)\delta (𝐱)`$, where $`D_+(t)a(t)t^{2/3}`$ is the growing mode and $`D_{}(t)t^1`$ is the decaying mode.

The above considerations apply to the evolution of a single Fourier mode of the density field $`\delta (𝐱,t)=D_+(t)\delta (𝐱)`$. What is more likely to be relevant, however, is the case of a superposition of waves, resulting from some kind of stochastic process in which the density field consists of a superposition of such modes with different amplitudes. A statistical description of the initial perturbations is therefore required, and any comparison between theory and observations will also have to be statistical. Many versions of the inflationary scenario for the very early universe (Guth 1981; Guth & Pi 1982; Linde 1982; Albrecht & Steinhardt 1982) predict the initial density fluctuations to take the form of a Gaussian random field in which the initial Fourier modes of the perturbation field have random phases.

The linearised equations of motion provide an excellent description of gravitational instability at very early times when density fluctuations are still small ($`\delta 1`$). The linear regime of gravitational instability breaks down when $`\delta `$ becomes comparable to unity, marking the commencement of the quasi-linear (or weakly non-linear) regime. During this regime the density contrast may remain small ($`\delta <1`$), but the phases of the Fourier components $`\delta _𝐤`$ become substantially different from their initial values resulting in the gradual development of a non-Gaussian distribution function even if the primordial density field is Gaussian owing to mode–mode interactions.

Perturbation theory fails when non-linearity develops, but the fluid treatment is intrinsically approximate anyway. A fuller treatment of the problem requires a solution of the Boltzmann equation for the full phase-space distribution of the system $`f(𝐱,𝐯,t)`$ coupled to the Poisson equation (7) that determines the gravitational potential. In cases where the matter component is collisionless, the Boltzmann equation takes the form of a Vlasov equation:

$$\frac{f}{t}=\underset{i=1}{\overset{3}{}}\left(\frac{\varphi }{x_i}\frac{f}{v_i}v_i\frac{f}{x_i}\right).$$
(20)
The fluid approach outline above can only describe cold material where the velocity dispersion of particles is negligible. But even if the dark matter is cold, there may be hot components of baryonic material whose behaviour needs also to be understood. Moreover, the fluid approach assumes the existence of a single fluid velocity at every spatial position. It therefore fails when orbits cross and multi-streaming generates a range of particle velocities through a given point.

### 2.4 The Zel’dovich Approximation

Fortunately the formation of the main elements of large-scale structure can nevertheless be understood using relatively simple tools, such as the Zel’dovich approximation (Zel’dovich 1970). Let the initial (i.e. Lagrangian) coordinate of a particle in this unperturbed distribution be $`𝐪`$. Now each particle is subjected to a displacement corresponding to a density perturbation. In the Zel’dovich approximation the Eulerian coordinate of the particle at time $`t`$ is

$$𝐫(t,𝐪)=a(t)[𝐪b(t)_𝐪\mathrm{\Phi }_0(𝐪)],$$
(21)
where $`𝐫=a(t)𝐱`$, with $`𝐱`$ a comoving coordinate, and we have made $`a(t)`$ dimensionless by dividing throughout by $`a(t_i)`$, where $`t_i`$ is some reference time which we take to be the initial time. Notice that particles in the Zel’dovich approximation execute a kind of inertial motion on straight line trajectories. The derivative on the right hand side is taken with respect to the Lagrangian coordinates. The dimensionless function $`b(t)`$ describes the evolution of a perturbation in the linear regime, with the condition $`b(t_i)=0`$, and therefore solves the equation

$$\ddot{b}+2\left(\frac{\dot{a}}{a}\right)\dot{b}4\pi G\rho b=0;$$
(22)
cf. equation (19). For a flat matter–dominated universe we have $`bt^{2/3}`$ as before. The quantity $`\mathrm{\Phi }_0(𝐪)`$ is proportional to a velocity potential, of the type introduced above, i.e. a quantity of which the velocity field is the gradient:

$$𝐕=\frac{d𝐫}{dt}H𝐫=a\frac{d𝐱}{dt}=a\dot{b}_q\mathrm{\Phi }_0(𝐪);$$
(23)
this means that the velocity field is irrotational. The quantity $`\mathrm{\Phi }_0(𝐪)`$ is related to the density perturbation in the linear regime by the relation

$$\delta =b_𝐪^2\mathrm{\Phi }_0,$$
(24)
which is a simple consequence of Poisson’s equation.

The Zel’dovich approximation, though simple, has a number of interesting properties. First, it is exact for the case of one dimensional perturbations up to the moment of shell crossing. As we have mentioned above, it also incorporates irrotational motion, which is required to be the case if it is generated only by the action of gravity (due to the Kelvin circulation theorem). For small displacements between $`𝐫`$ and $`a(t)𝐪`$, one recovers the usual (Eulerian) linear regime: in fact, equation (22) defines a unique mapping between the coordinates $`𝐪`$ and $`𝐫`$ (as long as trajectories do not cross); this means that $`\rho (𝐫,t)d^3r=\rho (t_i)d^3q`$ or

$$\rho (𝐫,t)=\frac{\rho (t)}{|J(𝐫,t)|},$$
(25)
where $`|J(𝐫,t)|`$ is the determinant of the Jacobian of the mapping between $`𝐪`$ and $`𝐫`$: $`𝐫/𝐪`$. Since the flow is irrotational the matrix $`J`$ is symmetric and can therefore be locally diagonalised. Hence

$$\rho (𝐫,t)=\rho (t)\underset{i=1}{\overset{3}{}}[1+b(t)\alpha _i(𝐪)]^1:$$
(26)
the quantities $`1+b(t)\alpha _i`$ are the eigenvalues of the matrix $`J`$ (the $`\alpha _i`$ are the eigenvalues of the deformation tensor). For times close to $`t_i`$, when $`|b(t)\alpha _i|1`$, equation (26) yields

$$\delta (\alpha _1+\alpha _2+\alpha _3)b(t),$$
(27)
which is just the law of perturbation growth obtained by solving equation (19).

At some time $`t_{sc}`$, when $`b(t_{sc})=1/\alpha _j`$, a singularity appears and the density becomes formally infinite in a region where at least one of the eigenvalues (in this case $`\alpha _j`$) is negative. This condition corresponds to the situation where two points with different Lagrangian coordinates end up at the same Eulerian coordinate. In other words, particle trajectories have crossed and the mapping between Lagrangian and Eulerian space is no longer unique. A region where the shell–crossing occurs is called a caustic. For a fluid element to be collapsing, at least one of the $`\alpha _j`$ must be negative. If more than one is negative, then collapse will occur first along the axis corresponding to the most negative eigenvalue. If there is no special symmetry, one therefore expects collapse to be generically one–dimensional, from three dimensions to two. We shall make use of the fact that in one-dimensional cases, the Zel’dovich solution is exact until the moment of shell-crossing.

The Zel’dovich approximation matches very well the evolution of density perturbations in full $`N`$–body calculations until the point where shell crossing occurs (Coles, Melott & Shandarin 1993). After this, the approximation breaks down completely. Particles continue to move through the caustic in the same direction as they did before. Particles entering a pancake from either side merely sail through it and pass out the opposite side. The pancake therefore appears only instantaneously and is rapidly smeared out. In reality, the matter in the caustic would feel the strong gravity there and be pulled back towards it before it could escape through the other side. Since the Zel’dovich approximation is only kinematic it does not account for these close–range forces and the behaviour in the strongly non–linear regime is therefore described very poorly. Furthermore, this approximation cannot describe the formation of shocks and phenomena associated with pressure. Attempts to understand properties of non-linear structure using a fluid description therefore resort to further approximations (Sahni & Coles 1995) to extend the approach beyond shell-crossing. One method involves the so–called adhesion model (Gurbatov, Saichev & Shandarin 1989). In this model one assumes that the particles stick to each other when they enter a caustic region because of an artificial viscosity which is intended to simulate the action of strong gravitational effects inside the overdensity forming there. This ‘sticking’ results in a cancellation of the component of the velocity of the particle perpendicular to the caustic. If the caustic is two–dimensional, the particles will move in its plane until they reach a one–dimensional interface between two such planes. This would then form a filament. Motion perpendicular to the filament would be cancelled, and the particles will flow along it until a point where two or more filaments intersect, thus forming a node. The smaller is the viscosity term, the thinner will be the sheets and filaments, and the more point–like will be the nodes. Outside these structures, the Zel’dovich approximation is still valid to high accuracy. See Shandarin & Zel’dovich (1989) for further discussion.

### 2.5 The Strongly Nonlinear Regime

Further into the non-linear regime, bound structures form. The baryonic content of these objects may then become important dynamically: hydrodynamical effects (e.g. shocks), star formation and heating and cooling of gas all come into play. The spatial distribution of galaxies may therefore be very different from the distribution of the (dark) matter, even on large scales. Attempts are only just being made to model some of these processes with cosmological hydrodynamics codes, such as those based on Smoothed Particle Hydrodynamics (SPH; Monaghan 1992), but it is some measure of the difficulty of understanding the formation of galaxies and clusters that most studies have only just begun to attempt to include modelling the detailed physics of galaxy formation. In the front rank of theoretical efforts in this area are the so-called semi-analytical models which encode simple rules for the formation of stars within a framework of merger trees that allows the hierarchical nature of gravitational instability to be explicitly taken into account.

The spatial distribution of particles obtained using the adhesion approximation represents a sort of “skeleton” of the real structure: non–linear evolution generically leads to the formation of a quasi–cellular structure, which is a kind of “tessellation” of irregular polyhedra having pancakes for faces, filaments for edges and nodes at the vertices . This skeleton, however, evolves continuously as structures merge and disrupt each other through tidal forces; gradually, as evolution proceeds, the characteristic scale of the structures increases. In order to interpret the observations we have already described, one can think of the giant “voids” as being the regions internal to the cells, while the cell nodes correspond to giant clusters of galaxies. While analytical methods, such as the adhesion model, are useful for mapping out the skeleton of structure formed during the non–linear phase, they are not adequate for describing the highly non–linear evolution within the densest clusters and superclusters. In particular, the adhesion model cannot be used to treat the process of merging and fragmentation of pancakes and filaments due to their own (local) gravitational instabilities, which must be done using full numerical computations.

## 3 Wave Mechanics and Structure Formation

### 3.1 Introduction

The combination of Eulerian fluid-based perturbation theory, Lagrangian approaches like the Zel’dovich approximation and brute-force $`N`$-body numerics has been successful at establishing a standard framework within which the basic features of large-scale structure can be described and understood. There remain significant obstacles to fuller analytic description. The three most interesting from the perspective followed in this paper are:

1. Perturbation Theory. Standard perturbation methods do not guarantee a density field that is everywhere positive. This basic problem is obvious in the standard case where there is a Gaussian distribution of initial fluctuations. When the variance of the distribution of $`\delta `$ is of order unity (its mean is, by definition, always unity) then a Gaussian distribution assigns a non-zero probability to regions with $`\delta <1`$, i.e. with $`\rho <0`$.
2. Shell-crossing. Methods such as the Zel’dovich approximation and its variations break down at shell-crossing, where the density becomes infinite.
3. Gas Pressure. Analytic techniques for modelling the effects of gas pressure are scarce, even in relatively simple systems such as the Lyman-$`\alpha `$ absorbing systems (Matarrese & Mohayee 2002). Fully numerical techniques, such as SPH, are generally expected to be the way forward but even they are not yet able to incorporate all relevant gas physics.

### 3.2 The Widrow–Kaiser Approach

A novel approach to the study of collisionless matter, with applications to structure formation, was suggested by Widrow & Kaiser (1993). It involves re-writing of the fluid equations given in Section 2 in the form of a non-linear Schrödinger equation. The equivalence between this and the fluid approach has been known for some time; see Spiegel (1980) for historical comments. Originally the interest was to find a fluid interpretation of quantum mechanical effects, but in this context we shall use it to describe an entirely classical system.

To begin with, consider the continuity equation and Euler equation for a curl-free flow in which $`𝐯=\varphi )`$, in response to some general potential $`V`$. In this case the continuity equation can be written

$$\frac{\rho }{t}+(\rho \varphi )=0.$$
(28)
It is convenient to take the first integral of the Euler equation, giving

$$\frac{\varphi }{t}+\frac{1}{2}(\varphi )^2=V,$$
(29)
which is usually known as the Bernoulli equation. The trick then is to make a transformation of the form

$$\psi =\alpha \mathrm{exp}(i\varphi /\nu ),$$
(30)
where $`\rho =\psi \psi ^{}=|\psi |^2=\alpha ^2`$; the wavefunction $`\psi (𝐱,t)`$ evidently complex. Notice that the dimensions of $`\nu `$ are the same as $`\varphi `$, namely $`[L^2T^1]`$. After some algebra it emerges that the two equations above can be written in one equation of the form

$$i\nu \frac{\psi }{t}=\frac{\nu ^2}{2}^2\psi +V\psi +P\psi ,$$
(31)
where

$$P=\frac{\nu ^2}{2}\frac{^2\alpha }{\alpha }.$$
(32)
To accommodate gravity we need to couple equation (31) to the Poisson equation in the form

$$^2\mathrm{\Phi }=4\pi G\psi \psi ^{},$$
(33)
taking $`V`$ to be $`\mathrm{\Phi }`$.

This system looks very similar to a Schrödinger equation, except for the extra “operator” $`P`$. The role of this term becomes clearer if one leaves it out of equation (31) and works backwards. The result is an extra term in the Bernouilli equation that resembles a pressure gradient. This is often called the “quantum pressure” that arises when one tries to understand a quantum system in terms of a classical fluid behaviour. Leaving it out to model a collisionless fluid can be justified only if $`\alpha `$ varies only slowly on the scales of interest. On the other hand one can model situations in which one wishes to model genuine effects of pressure by adjusting (or omitting) this term in the wave equation. Widrow & Kaiser (1993) advocated simply writing

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}^2\psi +m\varphi (𝐱)\psi ,$$
(34)
i.e. simply ignoring the quantum pressure term. In this spirit, the constant $`\mathrm{}`$ is taken to be an adjustable parameter that controls the spatial resolution $`\lambda `$ through a de Broglie relation $`\lambda =\mathrm{}/mv`$. In terms of the parameter $`\nu `$ used above, we have $`\nu =\mathrm{}/m`$, giving the correct dimensions for Planck’s constant. Note that the wavefunction $`\psi `$ encodes the velocity part of phase space in its argument through the ansatz

$$\psi (𝐱)=\sqrt{\rho (𝐱)}\mathrm{exp}[i\theta (𝐱)/\mathrm{}],$$
(35)
where $`\theta (𝐱)=𝐩(𝐱)`$, the local ‘momentum field’. This formalism thus yields an elegant description of both the density and velocity fields in a single function.

It should be clear how this approach bypasses immediately the first two items of the list given in Section 3.1. First, the construction (30) requires that $`\rho =|\psi |^2`$ is positive. Secondly, since particles are not treated as point-like entities with definite trajectories, shell-crossing does not have the catastrophic from in this approach that it does in the Zel’dovich approximation. Note that no singularities occur in the wavefunction at any time. Although the simplest ansatz (35) does assume a unique velocity at every fluid location, it is possible to construct more complex representations that allow for multi-streaming (Widrow & Kaiser 1993) although we shall not use them in this paper.

### 3.3 A Spherically-Symmetric Solution

In order to illustrate how the wave equation approach maps onto a fluid description, it is instructive to examine a familiar problem: a static spherically symmetric self-gravitating cloud in hydrostatic balance. This also allows us to explain how item (3) of the list in Section 3.1 can be addressed. We model the pressure gradients needed to hold an object in equilibrium against its self-gravity by the addition to the potential of a term of the form $`\kappa ^2|\psi |^2`$ (with $`\kappa `$ an appropriately-chosen constant), which converts the original equation (34) into a nonlinear Schrödinger equation:

$$i\nu \frac{\psi }{t}=\frac{\nu ^2}{2}^2\psi +\mathrm{\Phi }(𝐱)\psi +\kappa ^2|\psi |^2\psi +P\psi $$
(36)
(Sulem & Sulem 1999), where the final term is the quantum pressure. We will now look for solutions of the form

$$\psi =R\mathrm{exp}\left(\frac{i\varphi }{\nu }\right),$$
(37)
but we are looking for a static solution so obviously the fluid flow velocity should be zero. We therefore need to set $`\varphi `$ to be constant in the above transformation so that $`\varphi `$ can be zero everywhere. Without loss of generality we can exploit the global U(1) symmetry of quantum mechanics to set the phase to zero. Hence we can set $`\psi =R`$ and consequently $`^2\psi =^2R`$. The quantum pressure term has the form

$$P=\frac{\nu ^2}{2}\frac{^2R}{R},$$
(38)
so that this now cancels the first term on the right-hand-side of equation (36). In this event, and requiring a static solution for $`\psi `$ we get

$$\mathrm{\Phi }\psi +\kappa ^2|\psi |^2\psi =0$$
(39)
to be solved simultaneous with Poisson’s equation. We therefore obtain a simple Helmholtz equation for $`\mathrm{\Phi }`$:

$$^2\mathrm{\Phi }+\left(\frac{4\pi G}{\kappa ^2}\right)\mathrm{\Phi }=0.$$
(40)
The solution for spherically symmetry is straightforward:

$$\mathrm{\Phi }=\mathrm{\Phi }_c\frac{\mathrm{sin}\lambda r}{\lambda r},$$
(41)
where

$$\lambda ^2=\frac{4\pi G}{\kappa ^2}$$
(42)
and $`\mathrm{\Phi }_c`$ is constant. There is a particularly neat relation for $`\rho `$ using equation (40) yielding

$$4\pi G\rho +\frac{4\pi G}{\lambda ^2}\mathrm{\Phi }=0,$$
(43)
which gives

$$\rho (r)=\rho _c\frac{\mathrm{sin}\lambda r}{\lambda r}.$$
(44)
The resulting density profile will be recognized as that of an $`n=1`$ polytrope, i.e. a fluid in which the gas pressure is related to the density via $`p=K\rho ^2`$, with $`K`$ a constant.

The important point emerging from this analysis is that one can only obtain a static solution in the presence of the quantum pressure term. A self-consistent solution requires this term to be included and its effect is to stabilize the system. In order to apply the wave mechanical approach to structure formation, therefore, one has to be more careful with this term than Widrow & Kaiser (1993) in their equation (34). We investigate this issue further in the following Section.

## 4 One-dimensional Collapse

In this section we look at some simple one-dimensional examples of gravitational instability using the Schrödinger equation. We do this by comparing the solutions we obtain with the Zel’dovich approximation, which gives an exact solution up to the point of shell-crossing (Shandarin & Zel’dovich 1989). We also want to consider the effects of varying the parameter $`\nu `$ in the Schrödinger equation and how this changes the spatial and velocity resolution of the solution, as well as elucidating more clearly the effects the “quantum pressure” term.

### 4.1 The “Free-Particle” Solution to the Schrödinger Equation

The equation for the trajectory of a particle using the Zel’dovich approximation is identical in form to that of inertial (or ballistic) motion, with the displace replaced by the comoving co-ordinate $`x`$ and the time replaced by $`b`$ as defined in equation (22). The most obviously analogous approach to take with the Schrödinger equation is obtained be setting the potential $`V`$ equal to zero, replacing $`t`$ with $`b`$ and taking spatial derivatives with respect to the comoving co-ordinate $`x`$. Equation (31) then becomes

$$i\nu \frac{\psi }{b}=\frac{\nu ^2}{2}_𝐱^2\psi ,$$
(45)
or, in one-dimension,

$$i\nu \frac{\psi }{b}=\frac{\nu ^2}{2}\frac{^2\psi }{x^2}.$$
(46)
This is the “free-particle” Schrödinger equation and the solutions are given by

$$\psi (x,t)=\underset{k}{}A_ke^{(kx\omega b)},$$
(47)
with $`\omega =k^2\nu /2`$. The boundary conditions are given by the initial value of the wave function at $`b=0`$:

$$\psi (x,0)=\underset{k}{}A_ke^{ikx}.$$
(48)
The values of $`A_k`$ can be calculated using the Fast Fourier Transform of the initial wavefunction, and the value of the wavefunction at any future time can be calculated directly using equation (47).

In order to show the equivalence between the Zel’dovich approximation and Schrödinger method we consider the two solutions when the initial velocity potential is given by a single sine wave. In the Schrödinger equation we start with an initial wave function given by

$$\psi (x,0)=\sqrt{\rho _0}\mathrm{exp}(i\mathrm{cos}x/\nu ).$$
(49)
This is equivalent to setting the initial density contrast $`\delta =0`$ and the velocity $`v=\mathrm{sin}x`$. In one dimension, the Zel’dovich approximation for the position of a particle at time $`t`$ is given by,

$$x(q,t)=qb\frac{d\mathrm{\Phi }_0(q)}{dq},$$
(50)
cf. equation (21). Conservation of mass is enforced by

$$\rho _0dq=\rho dx$$
(51)
and therefore the density contrast is given by,

$$\delta =\left(\frac{x}{q}\right)^11=\left[1b\frac{d^2\mathrm{\Phi }_0(q)}{dq^2}\right]^11.$$
(52)
If the velocity potential $`\varphi _0(q)=\mathrm{cos}q`$ then the density contrast at time $`t`$ is then given by

$$\delta =(1+b\mathrm{cos}x)^11$$
(53)
and shell-crossing occurs when $`b=1`$ when the density contrast becomes infinite. Figure 1 shows the density contrast obtained using both of these methods at three different times and for three different values of $`\nu `$. When $`b=0.3`$ density growth is still approximately linear and the density contrast given by the Schrödinger and Zel’dovich methods are almost identical. When $`b=0.9`$, just before shell-crossing, the density growth has become highly non-linear. The very high density peaks seen in the Zel’dovich approximation are not seen in the Shroödinger solution.

The value of $`\nu `$ corresponds to a “de Broglie wavelength” which gives the spatial resolution of the solution. The effects of this can be seen most clearly when $`b=2.0`$ after shell-crossing has occurred. Although the Schrödinger and Zel’dovich methods show qualitatively different behaviour after shell-crossing, this is not physical, and the overall behaviour of the results is in any case similar, particularly for small values of $`\nu `$.

### 4.2 The Schrödinger equation with time-independent potential

Setting the potential $`V(x)=0`$ in the Schrödinger equation gives a good agreement with the Zel’dovich approximation. However, in both methods the particles continue on their trajectories after shell-crossing with no reference to the self-gravity of the collapsing system. Here we look at how the inclusion of a potential term in the Schrödinger equation might give a more realistic solution after shell-crossing. This approach is motivated by the realisation that, at least during linear perturbation growth, the gravitational potential corresponding to a density fluctuation barely changes. In a different context this has motivated the so-called “frozen-potential” approximation (Matarrese et al. 1992; Brainerd, Scherrer & Villumsen 1993; Bagla & Padmanabhan 1994), which is similar in spirit to the approach we adopt here.

The simplest way to include the potential is to include its initial value, but not to allow this to vary with time. For this example we again ignore the quantum pressure term. We therefore need to solve the equation

$$i\nu \frac{\psi }{t}=\frac{\nu ^2}{2}\frac{^2\psi }{x^2}+\mathrm{cos}(x)\psi .$$
(54)
It should be noted that in this example no account is taken of the expansion of the universe. This means that a direct comparison with this Zel’dovich approximation is not as straightforward as in the previous example and in particular the time coordinate is not the same in the two methods. However a comparison may be made by considering the growth at early times in the linear regime. In the Zel’dovich approximation the density contrast is given by

$$\delta b\frac{^2\varphi }{x^2}b\mathrm{cos}x,$$
(55)
and the velocity in comoving coordinates is given by

$$v=\frac{\mathrm{d}x}{\mathrm{d}t}\dot{b}\frac{\mathrm{d}\varphi _0}{\mathrm{d}x}\dot{b}\mathrm{sin}x.$$
(56)
In a non-expanding universe with a constant potential given by $`V(x)=V_0\mathrm{cos}x`$ the density contrast at early times is

$$\delta \frac{V_0}{2}t^2\mathrm{cos}x.$$
(57)
We can therefore compare the Zel’dovich approximation at time $`b`$ with the Schrödinger equation at time $`t=\sqrt{2b/V_0}`$. Starting with the initial wavefunction given by

$$\psi (x,0)=\sqrt{\rho _0},$$
(58)
and evolving the wavefunction according to equation (54) using Cayley’s method (Goldberg, Schey & Schwartz 1967) leads to the results shown in the first column of Figure 2. These are for $`\nu =0.1`$ at times corresponding to the values of $`b=0.3`$, $`0.9`$ and $`2.0`$ used before. Comparison with Figure 1 shows that the extensive smearing of density peaks which occurs after shell-crossing in both the “free-particle” and Zel’dovich solution is no longer seen in the case where the initial potential is included.

### 4.3 The Schrödinger equation with time-varying potential

In the second column of Figure 2 we show the solution to the Schrödinger equation with the full potential calculated using the Poisson equation (33). We start be considering the equation without the quantum pressure term. In this case we need to solve the coupled Schrödinger and Poisson equations:

$$i\nu \frac{\psi }{t}=\frac{\nu ^2}{2}\frac{^2\psi }{x^2}+V(x)\psi ;$$
(59)
$$\frac{^2V}{x^2}=4\pi G\psi \psi ^{}.$$
(60)
To set up the initial wavefunction we choose a single sine wave for the density profile,

$$\delta =\delta _0\mathrm{cos}x.$$
(61)
In the linear regime the density contrast at time $`t`$ is given by

$$\delta \delta _0\mathrm{exp}(\sqrt{4\pi G}t)\mathrm{cos}x$$
(62)
and the velocity is

$$v=\sqrt{4\pi G}\delta _0\mathrm{exp}(\sqrt{4\pi G}t)\mathrm{sin}x.$$
(63)
Accordingly we compare the Zel’dovich approximation at time $`b`$ with the Schrödinger equation at time $`t=1/\sqrt{4\pi G}\mathrm{ln}(b/\delta _0)`$. We start with the initial wavefunction,using equation (35), given by,

$$\psi (x,0)=\sqrt{\rho _0(1\delta _0\mathrm{cos}x)}\mathrm{exp}\left(i\sqrt{4\pi G}\delta _0\mathrm{cos}x/\nu \right),$$
(64)
and evolve the wavefunction according to equation (59). In this case the potential $`V(x)`$ is calculated by solving the Poisson equation numerically at each time step. The results are shown in the second column of Figure 2 for $`\nu =0.1`$ at times corresponding to $`b=0.3`$, $`0.9`$ and $`2.0`$ used before.

The third column of Figure 2 shows the results for the same initial wavefunction, but this time including the quantum pressure term, as in equation (31). The inclusion of this term makes a significant difference to the results at later times, particularly after shell-crossing, but is slight earlier on during the evolution. This means that useful approximations for the mildly non-linear regime can still be be obtained even if this term is ignored.

### 4.4 Solutions in an expanding background

Although we do not make further use of it in this paper, it is useful to end this section with a few comments about what arises when one places the system in an expanding background. To see what happens, let us define a scaled density $`\chi =\rho /\rho _0=(1+\delta )`$ and take $`\mathrm{\Omega }=1`$. The continuity equation then becomes

$$\frac{\chi }{a}+(\chi \varphi )=0,$$
(65)
where the velocity potential $`\varphi `$ is now such that $`𝐮=d𝐱/dt=\varphi `$ and $`a`$ is the scale factor. It is convenient to take the first integral of the Euler equation, giving

$$\frac{\varphi }{a}+\frac{1}{2}(\varphi )^2=\frac{3}{2a}(\varphi +\theta ),$$
(66)
where $`\theta =2\mathrm{\Phi }/3a^3H^2`$ and $`\mathrm{\Phi }`$ is the gravitational potential. After some more algebra the system again becomes a Schrödinger-like wave equation, but in $`a`$ rather than in $`t`$ and using $`\psi ^2=\chi `$, such that

$$i\nu \frac{\psi }{a}=\frac{\nu ^2}{2}^2\psi +V\psi +P\psi ,$$
(67)
with $`V=\varphi +\theta `$ and $`P`$ as before. It is relatively straightforward to incorporate the extra terms.

## 5 Discussion

In this paper we have set out an approach to the formation of cosmic structure from evolving cosmological density fluctuations that relies on a transformation of the evolution equations into a Schrödinger-Poisson system. One advantage of this new formalism was pointed out by Coles (2002): it yields a rather more convincing approach to understanding the origin of spatial intermittency in the galaxy distribution than that offered by Jones (1999). It also makes a connection in the underlying physics with other systems that display similar phenomena. In this paper we have argued further that it provides an analytical approach that offers benefits over standard approaches of both Lagrangian and Eulerian types. On the one hand, it always guarantees a positive density and on the other it does not produce singularities at shell-crossing events. It is clearly important to investigate both of these issues further. Szapudi & Kaiser (2003) have already completed a study of perturbation theory and its statistical ramifications within this approximation. Dealing fully with post-shell-crossing phenomena is beyond the scope of this paper as it requires a more sophisticated representation of the wave-function in velocity space; see Widrow & Kaiser (1993) for discussion.

In this paper we showed how this approach works in practice by applying it to a static spherically-symmetric case and to some examples of one-dimensional collapse. In doing this we looked at a “free-particle” solution of the system that resembles the Zel’dovich approximation in some respects. Incorporating a fixed gravitational potential yields a solution that has some similarities to the “frozen-potential” approximation (Matarrese et al. 1992). We think the Schrödinger formalism provides a useful tool for modelling fluctuations in the mildly non-linear regime. In future work we will look at its use in modelling the intergalactic medium in order to understand properties of Lyman-$`\alpha `$ absorption systems.

## Acknowledgements

We are grateful to Roya Mohayaee and Peter Watts for related discussions.