text:
# Entanglement in The Second Quantization Formalism

## I Introduction

Entanglement is a phenomenon that has become of great importance in modern applications of quantum mechanics . We usually think of entanglement as existing between different degrees of freedom of two or more particles. There has been a number of important advances recently in understanding entanglement in systems containing a small number of particles. If we are to use this appropriately, however, we need to extend our analysis to realistic systems containing a large number of particles. When it comes to large systems in quantum mechanics, the concept of a particle actually fades away and is replaced by the notion of “an excitation of a given mode of the field representing the particle”. Individual particles actually become indistinguishable (this is, of course, also true for a small number of identical particles, but indistinguishability may not play any significant role). In addition to that, the concepts of different particle statistics (fermions versus bosons) then also becomes directly relevant and it is important to understand its relation to entanglement. The most convenient and appropriate formalism to deal with all the issues involving a large number of systems is the second quantization. So we need to understand what entanglement means in this setting if we are to be able to harness solid state and condensed matter systems for information processing purposes. Another benefit of the second quantization is that it is also the correct formalism for the (relativistic) quantum field theory. We can therefore exploit some already known results in this field in the hope of reaching a more complete understanding of entanglement.

In this paper we will see that entanglement is highly dependent on the choice of modes  (for an alternative treatment see ). The choice of modes is mainly dictated by the physics of the given situation, and the general rule is that in order for modes to be entangled they need to interact in some way. The crux is that what happens to be an interacting hamiltonian for one set of modes may not be so for a different set of modes. Alternatively stated, the vacuum state of one set of modes, may not be the vacuum state when the modes are changed. It is this simple fact that will be discussed in this paper and shown to be at the root of many different phenomena. More interestingly, however, we will see that the modes need not interact in order to produce entanglement - entanglement can be a consequence of indistinguishability of quantum particles and the existence of quantum statistics.

In the next section we review the formalism of second quantization. In section $`3`$ we present a very simple illustrative example of two harmonic oscillators and show how entanglement arises through their interaction. The conclusions of this section are very general and will then be applied to the Bose condensation as well as the Unruh effect in section $`4`$. In section $`5`$ we show that interaction is, however, not necessary to generate entanglement. We present a Fermi gas of non-interacting particles at zero temperature and show that the spin degree of freedom of individual fermions can be entangled even though the fermions themselves are non-interacting. We show that an analogous collection of non-interacting bosons, on the other hand, cannot generate any such entanglement in the internal degrees of freedom. We finally conclude and discuss some future directions in section $`6`$.

## II Second Quantization

First we briefly review the formalism of the second quantization . We will approach the second quantization formalism from the perspective of many particle systems, rather than as a means of achieving relativistic quantum theory. Suppose we have the Schrödinger equation for a single quantum system in one dimension (for simplicity):

$$\widehat{H}(x)\mathrm{\Psi }(x,t)=i\mathrm{}\frac{}{t}\mathrm{\Psi }(x,t)$$
(1)
Note that although this will lead to the non-relativistic second quantization, much of what we say will be true relativistically as well. A formal solution to this equation is

$$\mathrm{\Psi }(x,t)=\underset{n}{}b_n(t)\psi _n(x).$$
(2)
As usual, by substituting this back into the Schrödinger equation, we obtain its time independent version

$$\widehat{H}(x)\psi _n(x)=E_n\psi _n(x),$$
(3)
where $`E_n`$ are the corresponding energies. In order to convert this into a many-particle equation, we apply the formalism of second quantization, which means that we effectively have to “upgrade” the wave-function $`\mathrm{\Psi }`$ into an operator. Formally, we write

$$\widehat{\mathrm{\Psi }}(x,t)=\underset{n}{}\widehat{b}_n(t)\psi _n(x)$$
(4)
The $`\widehat{b}_n`$ operators are the well-known annihilation operators. The conjugate of the above equation becomes

$$\widehat{\mathrm{\Psi }}^{}(x,t)=\underset{n}{}\widehat{b}_n^{}(t)\psi _n^{}(x)$$
(5)
where $`b_n^{}`$ are the creation operators. In this paper we only consider the time independent creation and annihilation field operators. Nothing much would change conceptually if the field operators were to be time dependent, although mathematically the whole analysis would become much more involved. The Hamiltonian also has to be second quantized, and the new Hamiltonian is given by the average of the old (first quantized) Hamiltonian

$$\stackrel{~}{H}=𝑑x\widehat{\mathrm{\Psi }}^{}(x,t)\widehat{H}(x)\widehat{\mathrm{\Psi }}(x,t)$$
(6)
By invoking the orthogonality rules, $`\psi _n|\psi _m=\delta _{mn}`$, the second quatized Hamiltonian becomes

$$\stackrel{~}{H}=\underset{n}{}E_n\widehat{b}_n^{}\widehat{b}_n$$
(7)
and this is the same as a set of independent harmonic oscillators. Now, the fields here are clearly non-interacting and therefore we may conclude that the $`b`$ modes are disentangled. But, we have to be careful. While there is, clearly, no entanglement in the $`b`$ modes, there may be entanglement between some other modes. To illustrate this let us have a look at the simplest case of two harmonic oscillators. Note that the index $`n`$ can also contain the internal degrees of freedom (such as spin or polarization), however, not until section $`5`$ will we need to take this complication into account.

## III Interacting Harmonic Oscillators

Suppose that we have two non-interacting harmonic oscillators so that the total Hamiltonian is (we will omit hats from now on):

$$H=\mathrm{}(\omega \lambda )b_1^{}b_1+\mathrm{}(\omega +\lambda )b_2^{}b_2$$
(8)
Since oscillators are non-interacting their eigenstates are, therefore, composed of disentangled direct products of the eigenstates of the individual harmonic oscillators. For example, the ground state is just a product of the ground states of the oscillator $`1`$ and the oscillator $`2`$, $`|0_b=|0_{b_1}|0_{b_2}`$, so that it is annihilated by both the annihilation operators of individual oscillators

$`b_1|0_b`$ $`=`$ $`0`$ (9)
$`b_2|0_b`$ $`=`$ $`0.`$ (10)
However, let us make a change of modes and see what happens to the new eigenstates. Let the new creation and annihilation operators be

$`a_1={\displaystyle \frac{b_1+b_2^{}}{\sqrt{2}}}`$ (11)
$`a_2={\displaystyle \frac{b_1^{}b_2}{\sqrt{2}}}`$ (12)
(We could have made a more general change, but the one above is sufficient to illustrate the point here). The original Hamiltonian now becomes

$$H=\mathrm{}\omega (a_1^{}a_1+a_2^{}a_2)+\mathrm{}\lambda (a_1^{}a_2+a_2^{}a_1)$$
(13)
and it is transparent that the $`a`$ modes actually represent two coupled harmonic oscillators. It is, therefore, not surprising that the ground state of this system is actually entangled when written in the basis of the original harmonic oscillators. It is given by (up to a normalization which is not relevant for our discussion)

$$|0_a\underset{n}{}[i\mathrm{tanh}\frac{4\lambda }{\omega }]^n|n_{b_1}|n_{b_2},$$
(14)
where $`|n_{b_i}`$ is the state containing $`n`$ quanta in the oscillator $`b_i`$. This state is clearly entangled with respect to eigenstates of $`b`$ s, and it is written in the Schmidt decomposition form making it easy to detect entanglement. The amount of entanglement is obtained by tracing out one oscillator and calculating the entropy of the resulting mixture , which is known as the entropy of entanglement. The resulting entanglement is

$$E=\frac{2\mathrm{ln}(\mathrm{coth}\frac{4\lambda }{\omega })}{\frac{1}{\mathrm{tanh}^2\frac{4\lambda }{\omega }}1}\mathrm{ln}(1\mathrm{tanh}^2\frac{4\lambda }{\omega }).$$
(15)
To make the result look more familiar to us (and to prepare ourselves for the next section) we introduce the following change of variables

$$[\mathrm{tanh}\frac{4\lambda }{\omega }]^2=e^{\frac{\mathrm{}\omega }{kT}},$$
(16)
where we think of $`T`$ as the corresponding (effective) temperature; the entanglement is now given by

$$E=\frac{\mathrm{}\omega }{kT}\frac{1}{e^{\frac{\mathrm{}\omega }{kT}}1}\mathrm{ln}(1e^{\frac{\mathrm{}\omega }{kT}})$$
(17)
which is just the entropy of a thermal state of a harmonic oscillator, i.e. it is the familiar Planck distribution. We see that if $`T0`$ (i.e. the interaction tends to zero, $`\lambda 0`$), then $`E0`$, and entanglement disappears as expected. The fact that this is a thermal spectrum features strongly in the Unruh effect and the Hawking radiation from a black hole  as will be seen in the next section.

Let us now generalize our consideration to an arbitrary field which can be written as superposition of the creation and annihilations operators for the various modes

$$\psi (x)=\underset{n}{}(\psi _n(x)a_n+\psi _n^{}(x)a_n^{}).$$
(18)
Here we basically have infinitely many harmonic oscillators to sum over. The vacuum state for the $`a`$ operators is defined as $`a_n|0=0`$ for all $`n`$. Now, suppose that we expand the field in terms of a different set of creation and annihilation operators, like so

$$\psi (x)=\underset{n}{}(\stackrel{~}{\psi }_n(x)b_n+\stackrel{~}{\psi }_n^{}(x)b_n^{}).$$
(19)
The new operators can always be expanded as a combination of the old operators as

$`b_m`$ $`=`$ $`{\displaystyle \underset{n}{}}\alpha _{mn}a_n+\beta _{mn}a_n^{}`$ (20)
$`a_m`$ $`=`$ $`{\displaystyle \underset{n}{}}\alpha _{mn}^{}b_n\beta _{mn}^{}b_n^{}`$ (21)
(the new mode functions $`\stackrel{~}{\psi }`$ can also be expressed in terms of the old ones , but this is not important for our discussion). Due to the unitarity requirement, the coefficients of transformation have to posses the following properties:

$`{\displaystyle \underset{k}{}}`$ $`(\alpha _{ik}\alpha _{jk}^{}\beta _{ik}\beta _{jk}^{})=\delta _{ij}`$ (22)
$`{\displaystyle \underset{k}{}}`$ $`(\alpha _{ik}\beta _{jk}\beta _{ik}\alpha _{jk})=0.`$ (23)
Now, the point of our argument here is that the new vacuum, defined as $`b_m|\stackrel{~}{0}=0`$ (for all $`m`$), will not in general be the same as the old vacuum (i.e. it will not be annihilated by the old annihilation operators, $`a_n`$). The new vacuum will, for example, have the following number of particles

$$\stackrel{~}{0}|a_i^{}a_i|\stackrel{~}{0}=\underset{j}{}|\beta _{ij}|^2$$
(24)
with respect to the $`i`$th mode of $`a`$ (and this means that the vacuum of the $`b`$ modes contains some particles with respect to the $`a`$ modes). Therefore, since the vacuum $`|\stackrel{~}{0}`$ is an overall pure state, we can conclude that the amount of entanglement between the $`i`$th mode and all the other modes taken together is given by the von Neumman entropy of entanglement

$$E_i=\underset{j}{}|\beta _{ij}|^2\mathrm{ln}|\beta _{ij}|^2$$
(25)
We see that this entanglement is zero only when one of the $`\beta `$s is a unity and the rest are zero, i.e. the new modes are the same as the old ones. The maximum entanglement, on the other hand, is achieved when all the betas are equal. Our two modes at the beginning of the section are a very special case of this analysis.

In summary, there are two main messages in this section: one is that entanglement is highly dependent on the way we decide to divide the whole system into different modes; two is that interaction in general gives rise to entanglement. There are numerous examples of various kinds of mode transformations in quantum optics (for example, in two-mode squeezing), condensed matter physics as well as relativistic quantum field theory which lead to entanglement generation as described above. We discuss a few of them next.

## IV Bose Condensates and The Unruh Effect

We first briefly review the second quantization treatment of interacting systems and then apply it to several cases. The Hamiltonian in the first quantized form is now given by

$$H(x,x^{})=H_0(x)+V(x,x^{})$$
(26)
where $`H_0(x)`$ is the basic (free) Hamiltonian without interactions and the interaction itself is encoded in the term $`V(x,x^{})`$. We already know how to second quantize $`H_0(x)`$, and so we only need to be able to second quantize the interaction Hamiltonian. We will do this now with the intention of obtaining a typical Hamiltonian for the Bose-Einstein gas. This is done by performing the following transformation

$$\stackrel{~}{V}=𝑑V\mathrm{\Psi }^{}(x,t)\mathrm{\Psi }^{}(x^{},t)V(x,x^{})\mathrm{\Psi }(x,t)\mathrm{\Psi }(x^{},t).$$
(27)
As before, we need to switch from the creation and annihilation operators in the position picture to the creation and annihilation operators in the momentum picture. This is performed in the same way as we did for the non-interacting fields, as in eq. (4). After some manipulations and approximations we can arrive at

$$\stackrel{~}{H}=\underset{p}{}\frac{p^2}{2m}a_p^{}a_p+\underset{p_1,p_2,q}{}\stackrel{~}{V}(q)a_{p_1+q}^{}a_{p_2q}^{}a_{p_2}a_{p_1}.$$
(28)
So, the first term of the Hamiltonian is (as before) a collection of harmonic oscillators and the second term signifies their interaction. The interaction term tells us that two particles of momenta $`p_1`$ and $`p_2`$ collide to produce new particles with momenta $`p_1+q`$ and $`p_2q`$ respectively. The values of the final momenta are such as to obey the momentum conservation law. $`\stackrel{~}{V}(q)`$ is the strength of interaction which is derived from the original potential $`V(x,x^{})`$. The exact form of this is unimportant. In fact, it is sometimes customary to assume that the interaction is independent of the momentum change $`q`$ and to write (after omitting the tilde)

$$H=\underset{p}{}\frac{\mathrm{}p^2}{2m}a_p^{}a_p+V\underset{p_1,p_2,q}{}a_{p_1+q}^{}a_{p_2q}^{}a_{p_2}a_{p_1}.$$
(29)
This becomes the same as the Hamiltonian for a degenerate, but almost ideal Bose gas . The creation and annihilations operators, of course, obey the usual bosonic commutation relations. In order to solve and find the eigenvalues of this Hamiltonian we need to convert it into a bunch of non-interacting harmonic oscillators (or modes). The transformations that achieve this are called the Bogoliubov transformations and are given by

$`a_p`$ $`=`$ $`u_pb_p+v_pb_p^{}`$ (30)
$`a_p^{}`$ $`=`$ $`u_pb_p^{}+v_pb_p`$ (31)
where the coefficients $`u`$ and $`v`$ (functions of the momentum $`p`$) have to be chosen so that the Hamiltonian for $`b`$ modes has no interacting parts (see  for an experimental confirmation of the Bogoliubov transformations in a Bose condensate). The resulting Hamiltonian, written in terms of the $`b`$ modes, is given by

$$H=E_0+\underset{p0}{}E(p)b_p^{}b_p$$
(32)
where $`E_0`$ is, for us an unimportant, constant and $`E(p)`$ is the energy of new modes (all the details, which are irrelevant for the analysis here, can be found in ). Here we have a collection of decoupled (i.e. non-interacting) modes as in the previous section when we discussed the example of two harmonic oscillators. Again, the ground state is an entangled state in terms of the $`a`$ modes and is given by

$$|\psi _0=\underset{p0}{}\frac{1}{u_p}\underset{i=0}{\overset{\mathrm{}}{}}\{\frac{v_p}{u_p}\}^i|n_p=i;n_p=i.$$
(33)
In this state only the $`p`$ and $`p`$ modes are entangled and different momenta are not coupled. The amount of entanglement is again easy to compute to be (c.f. eq. (17))

$$E=\underset{p0}{}\{\frac{\mathrm{ln}(u_p/v_p)^2}{(u_p/v_p)^21}\mathrm{ln}(1(v_p/u_p)^2)\}$$
(34)
This is a sum of entanglements in each individual pair of modes with the opposite momenta. The entanglement is highest when $`u_pv_p`$. This happens in the low momentum regime, more precisely when $`p<<mc`$. In the other regime, when the momenta are high we have that $`u_p1`$ and $`v_p0`$ and so entanglement drops to zero.

A very closely related, but even more intriguing, case of mode transformations generating entanglement can be found in the domain of relativistic quantum field theory. Next we discuss the accelerated vacuum leading to the Unruh effect . The key point in this effect is that the vacuum state for an uniformly moving observer in the flat Minkowski space, which is a continuum of decoupled, disentangled, harmonic oscillators, is actually seen by an accelerated (Rindler) observer to be entangled. The entanglement is created because the creation operator for a particle of momentum $`p`$ for the accelerated observer needs to be expanded as a sum of the creation operations of the uniform observer involving both the $`p`$ and $`p`$ momenta. More precisely, Minkowski quanta (i.e. those detected by an inertial observer) are obtained by quantizing the positive frequency solutions of the Klein-Gordon equation; however, once we change the time parameter which is what happens for the accelerated Rindler observer, we obtain both the set of positive as well as the negative frequency solutions. These, when quantized, lead to completely different set of raising and lowering operators, and therefore a new type of quanta. Interestingly, the transformation of the mode operators between the stationary and the accelerated observer is actually of the Bogoliubov type discussed already in the case of a Bose condensate. The vacuum state of the accelerated observer is, therefore, an entangled superposition of the number states for $`p`$ and $`p`$ of the original uniform field . The interesting observation is that the spectrum of the reduced density matrix (when we trace out the negative momenta for example) is Planckian (just like in the case of two interacting harmonic oscillators) whose temperature is given by $`T=a/2\pi k`$, where $`a`$ is the acceleration of the observer and $`k`$ is the Boltzmann constant. There is, therefore, entanglement between positive and negative modes due to acceleration, and this entanglement is proportional to the amount of acceleration. The amount of entanglement, more precisely, is the same as in eq. (17) with the appropriate value for the temperature.

We would just like to close this section by briefly mentioning an effect closely related to the entanglement in the accelerated vacuum and the corresponding thermal bath. This effect was discovered by Hawking and it explains the radiation from a black hole and the fact that we can associate a temperature and entropy with it . The two effects - accelerated observers and the black hole radiation - are, loosely speaking, related because of Einstein’s equivalence principle, which stipulates that an accelerated observer is the same as an observer in a gravitational field whose strength is that same as the value of acceleration. So, someone sitting close to a black hole for example, should also see a radiation (just like an accelerated observer does) which has a black body spectrum. The temperature of the radiation can be shown to be $`T=\kappa /2\pi k`$, where $`\kappa `$ is the surface gravity. Hawking  showed in $`1974`$ that this exactly is the case. The entanglement here is generated because the incoming modes which contain only positive frequency waves will be converted into the outgoing waves, which will depend also on the negative frequency waves due to the metric change during the collapse of matter to form a black hole (there is a pair creation of this type under more general metric changes in the expanding universe as first discussed by Parker in . He found that particle creation occurs because the positive and negative parts of the field become mixed during the expansion, such that the creation and annihilation operators at one time $`t`$ are a linear (Bogoliubov) combination of those at an earlier time.). So, the situation here can mathematically be described by the same Bogoliubov transformation discussed before . Interestingly, the entropy of a black hole can also be seen to arise due to entanglement between the modes outside the black hole and modes inside the black hole  (this is a static, rather than a dynamic explanation of the entropy of the black hole). There are many issues involved here which have not been resolved, but which are beyond the motivation of the current paper, and I direct the interested reader to consult Wald’s excellent review on this and related subjects .

The above examples illustrate the entanglement properties in the second quantization and its dependence of the choice of modes. A collection of harmonic oscillators that does not interact becomes interacting when we change to a different set of modes. This change may be “passive”, like in the case when we diagonalize the Bose condensate Hamiltonian, or it could be “active”, like when a stationary observer accelerates. But, the question we would really like to ask in this paper is whether entanglement could arise without any interaction between different modes. So far, we have only considered spinless bosons, and the particle statistics has not been crucial in any respect. We would now like to analyze the interplay between the internal and external degrees of freedom in fermions and bosons, which are otherwise non-interacting.

## V Non-interacting Fermions and Bosons

Suppose we have a collection of non-interacting electrons, in a metal, for example. The electrons all move with various different momenta $`p`$ and Pauli’s exclusion principle tells us that there can be only two electrons with the same value of momentum (one with spin, $`s`$, pointing “up” and the other with the spin pointing “down”). The state of these electrons is well approximated by the overall state (the so called Fermi model of a metal):

$$|\varphi _0=\underset{s,p}{}b_s^{}(p)|0.$$
(35)
The fact that these are Fermions (as opposed to Bosons) is reflected in the anticommutation relations which we have to impose between the annihilation (and creation) operators

$$[b_s^{}(p),b_t(q)]_+=\delta _{st}\delta (pq).$$
(36)
Obviously there are no correlations between the momenta of electrons. This means that if we find an electron and measure its momentum to be $`p`$, then if we measure another electron’s momentum $`q`$, its value will be completely independent of the first measurement. That is exactly what the product of creation operators in the above state signifies. However, there is a great deal of entanglement in the above state, we just need to look for it more carefully. Suppose that we measure an electron at the position $`r`$ and another one at the position $`r^{}`$. We will assume, for the sake of simplicity, that we have arbitrarily precise detectors (this is not a point of fundamental importance in the analysis). Suppose further that one electron has spin up, can we infer the direction of the spin of the other electron? The question is, are electrons at different locations correlated in in their spin directions? In order to calculate this, we need to compute the density matrix describing the spin state of two electrons. It is (by definition) given by

$$\rho _{ss^{};tt^{}}=\varphi _0|\psi _t^{}^{}(r^{})\psi _t^{}(r)\psi _s^{}(r^{})\psi _s(r)|\varphi _0$$
(37)
where $`\psi _t^{}(r)`$ creates a particle of spin $`t`$ at the location $`r`$. This is the same as computing the elements $`|sts^{}t^{}|`$ of the two-electron density matrix. Assuming (as we usually do) that

$$\psi _s(r)=\frac{d^3p}{(2\pi )^3}e^{ipr}b_s(p)$$
(38)
we obtain the following form of the density matrix

$`\rho _{ss^{};tt^{}}`$ $`=`$ $`{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{d^3p^{}}{(2\pi )^3}}{\displaystyle \frac{d^3q^{}}{(2\pi )^3}}e^{i(pp^{})r}e^{i(qq^{})r^{}}\times `$ (39)
$`\times `$ $`\varphi _0|b_t^{}^{}(p)b_t^{}(q)b_s^{}(q^{})b_s(p^{})|\varphi _0.`$ (40)
Now we evaluate the last term to be

$`\varphi _0|b_t^{}^{}(p)b_t^{}(q)b_s^{}(q^{})b_s(p^{})|\varphi _0`$ $`=`$ $`(2\pi )^6\delta _{st}\delta _{s^{}t^{}}\delta (pq)\delta (p^{}q^{})`$ (41)
$``$ $`(2\pi )^6\delta _{s^{}t}\delta _{st^{}}\delta (p^{}q)\delta (pq^{}).`$ (42)
This form is intuitively clear as we have both the direct term, i.e. the term where the particle with momentum $`p`$ changes to $`q`$ and $`p^{}`$ to $`q^{}`$, and the exchange term, such that the particle with momentum $`p`$ changes to $`q^{}`$ and $`p^{}`$ to $`q`$. Note that because of the conservation of momenta these have to be the same (there is no interaction to allow the change of momenta or conversion of one particle into two for example). This now enables us to compute the density matrix, by observing that the density of electrons is given by $`n=_0^{p_F}𝑑p^3/(2\pi )^3`$ (where $`p_F`$ is the Fermi momentum),

$$\rho _{ss^{};tt^{}}=n^2(\delta _{st}\delta _{s^{}t^{}}\delta _{s^{}t}\delta _{st^{}}f^2(rr^{})),$$
(43)
where

$$f(rr^{})=_0^{p_F}\frac{d^3p}{(2\pi )^3}e^{ip(rr^{})}=\frac{3j_1(k_F|rr^{}|)}{k_F|rr^{}|}$$
(44)
is the “exchange interaction” term, and $`j_1(k_F|rr^{}|)`$ is the spherical Bessel function with $`k_F=p_F/\mathrm{}`$. This form of $`\rho `$ contains all $`16`$ elements needed for the spin density matrix. This $`4`$ by $`4`$ density matrix is (up to a normalization) given by

$$\rho _{12}=\left(\begin{array}{cccc}1f^2& 0& 0& 0\\ 0& 1& f^2& 0\\ 0& f^2& 1& 0\\ 0& 0& 0& 1f^2\end{array}\right),$$

where the subscript $`12`$ signifies that there are two electrons, and the basis is the usual $`|,|,|,|`$. We first have to decide if this matrix is entangled. We perform the partial transposition and compute the resulting eigenvalues . They are $`1f^2`$ (doubly degenerate), $`1`$ and $`12f^2`$. The first two eigenvalues are positive as $`0f1`$, but the last one can be negative. The condition is that $`f^2>1/2`$, and in this region we have entanglement. Therefore the density matrix is entangled for the region of distances $`0|rr^{}|<r_e`$, where $`r_e`$ is the value such that $`j_1^2(k_Fr_e)=1/2`$ (the value is about $`r_e\pi \mathrm{}/8p_F`$ as expressed in terms of Fermi’s momentum).

We can compute the entanglement of formation  of the spins for a general state $`\rho _{12}`$ (as a function of $`f`$). It is given by

$`E(\rho _{12})`$ $`=`$ $`{\displaystyle \frac{1+\sqrt{1(\frac{2f^21}{2f^2})^2}}{2}}\mathrm{ln}{\displaystyle \frac{1+\sqrt{1(\frac{2f^21}{2f^2})^2}}{2}}`$ (45)
$``$ $`{\displaystyle \frac{1+\sqrt{1+(\frac{2f^21}{2f^2})^2}}{2}}\mathrm{ln}{\displaystyle \frac{1+\sqrt{1+(\frac{2f^21}{2f^2})^2}}{2}}`$ (46)
Entanglement is maximal when $`f=1`$ and it is equal to $`\mathrm{ln}2`$. The resulting state is the spin singlet state $`||`$. Otherwise, as the distance decreases and $`f`$ drops, the more of the triplet state becomes mixed in with the singlet and entanglement drops as a consequence. Beyond the distances of about $`\pi \mathrm{}/8p_F`$, all entanglement disappears. A reasonable estimate for the Fermi momentum in a metal is $`10^{27}kg10^4ms^1`$, which gives a distance of about $`\mathrm{\Delta }r>10^7`$ meters beyond which no two electrons are entangled in spin. Below this distance, every two electrons in a conductor are entangled. If we have a much bigger momentum, such as in a neutron star, then the entanglement distance becomes very small. For a typical neutron star, this distance would be about $`10^{15}`$ meters, which tells us that there is no entanglement there between any two neutrons in a neutron star. Therefore, the higher the Fermi momentum the smaller the “entanglement distance” between the electrons. Note that it is a separate issue whether this entanglement present in the spin degrees of freedom can actually be used to process information, such as in the teleportation protocol. In order to teleport, for example, the electrons need to be separated sufficiently well after being detected. How feasible this is in practice will be a subject of a separate investigation.

The disappearance of entanglement is intuitively expected from the following simple argument. Quantum statistics becomes important only when the spatial wavefunctions of individual fermions have a substantial overlap. The wavefunction spread of fermions is roughly given by the wavelength which can be estimated to be $`\mathrm{}/p_F`$. Therefore beyond this distance we expect particle statistics not to play such a strong role. Let us look at this in the first quantization to make the case even clearer. Suppose that the total (un-normalized) state of two electrons is:

$$|01|\psi _1,\psi _2$$
(47)
where the $`|\psi _1,\psi _2`$ are the spatial wavefunctions of the two electrons, and $`0,1`$ are the values of spin (the first position always belongs to the first electron and the second to the second electron). Suppose that $`\psi _1|\psi _2=ϵ`$. Now this is not an allowed state if the particles are truly indistinguishable and in the case of fermions it needs to be antisymmetrized. The corresponding allowed state is

$$|01|\psi _1,\psi _2|10|\psi _2,\psi _1$$
(48)
Now, this state can be entangled in spin. To check this, we trace out the spatial degree of freedom to obtain the following density matrix (again unnormalised)

$$|0101|+|1010|ϵ^2(|0110|+|1001|)$$
(49)
It can be checked, again using the partial transposition criterion, that this state is entangled, unless $`ϵ=0`$, which means that the two electrons have no spatial overlap. Therefore we see that the quantum statistics plays a crucial role for entanglement.

Classical correlations, on the other hand, are still present even beyond this distance. This is true in the above first quantized model, as well as in the model of many non-interacting fermions. We can quantify them using the mutual information, $`I(\rho _{12})=S(\rho _1)+S(\rho _2)S(\rho _{12})`$, in the state in eq. (43) . Here $`\rho _1`$ and $`\rho _2`$ are the reduced density matrices derived from $`\rho _{12}`$. Computing this for our state of two electrons we obtain

$`I(\rho _{12})`$ $`=`$ $`2\mathrm{ln}2+2{\displaystyle \frac{1f^2}{42f^2}}\mathrm{ln}{\displaystyle \frac{1f^2}{42f^2}}`$ (50)
$`+`$ $`{\displaystyle \frac{1}{42f^2}}\mathrm{ln}{\displaystyle \frac{1}{42f^2}}+{\displaystyle \frac{12f^2}{42f^2}}\mathrm{ln}{\displaystyle \frac{12f^2}{42f^2}}`$ (51)
We see that the classical correlations only vanish when $`f=0`$ as we then have an equal mixture of all four states $`|,,|,,|,,|,`$. Otherwise, spins are always (classically) correlated.

The entanglement in spin when space is traced out is of a similar (but somewhat dual) nature to the spin-to-space entanglement transfer in the double beam-splitter set up in . The entanglement, which in this experiment by Omar el. al. initially exists in the spin degrees of freedom, is, due to the quantum statistical nature of particle interaction, transferred to space. This also is a purely quantum statistical effect as there is no other interaction between the internal and external degrees of freedom (although there is an interaction between the spatial degrees of freedom of different system involved).

Note that the presented method here allows us to extract more complicated information such as correlations between three and more electrons. These would be increasingly more difficult to treat mathematically, but the calculations can be performed in principle. We could, for example, check if genuine multi-electron correlations exist and what are the relevant distances to allow such correlations. Let us briefly summarize the calculation: we have three electrons at positions $`r,r^{}r^{\prime \prime }`$ and they have some values of spin. We would like to extract the spin density matrix, which is now at $`8`$ by $`8`$ square matrix. The form of the matrix can actually be deduced without much effort from the simple “pictorial” representation. We have the following terms

$$\rho _3=\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t& t^{}& t^{\prime \prime }\end{array}\right)\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t& t^{\prime \prime }& t^{}\end{array}\right)\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t^{}& t& t^{\prime \prime }\end{array}\right)\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t^{\prime \prime }& t^{}& t\end{array}\right)+$$

$$\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t^{}& t^{\prime \prime }& t\end{array}\right)+\left(\begin{array}{ccc}s& s^{}& s^{\prime \prime }\\ t^{\prime \prime }& t& t^{}\end{array}\right)$$

There are six terms as there are six possible arrangements; every spin in the density matrix can initially be paired with any of the three spins. The next spin has only two possibilities and the last one has only one remaining choice. The choice of signs should also be clear: the first term is the direct term and therefore positive. The next three terms involve one double exchange and are therefore negative, and the last two terms involve triple exchanges and are therefore positive. The density is, therefore, given by

$`\rho (s,s^{},s^{\prime \prime };t,t^{},t^{\prime \prime })`$ $`=`$ $`\varphi _0|\psi _{t^{\prime \prime }}^{}(r^{\prime \prime })\psi _t^{}^{}(r^{})\psi _t^{}(r)\psi _s(r^{})\psi _s^{}(r)\psi _{s^{\prime \prime }}(r^{\prime \prime })|\varphi _0`$ (52)
$`+`$ $`n^3(\delta _{st}\delta _{s^{}t^{}}\delta s^{\prime \prime }t^{\prime \prime }`$ (53)
$``$ $`f_1f_2\delta _{st}\delta _{s^{}t^{\prime \prime }}\delta s^{}t^{\prime \prime }`$ (54)
$``$ $`f_1f_3\delta _{st^{}}\delta _{s^{}t}\delta s^{\prime \prime }t^{\prime \prime }`$ (55)
$``$ $`f_2f_3\delta _{st^{\prime \prime }}\delta _{s^{}t^{}}\delta s^{\prime \prime }t`$ (56)
$`+`$ $`f_1f_2f_3\delta _{st^{}}\delta _{s^{}t^{\prime \prime }}\delta s^{\prime \prime }t`$ (57)
$`+`$ $`f_1f_2f_3\delta _{st^{\prime \prime }}\delta _{s^{}t}\delta s^{\prime \prime }t^{})`$ (58)
There are three different $`f`$ functions now as we have three arguments $`|rr^{}|`$, $`|rr^{}|`$ and $`|r^{}r^{\prime \prime }|`$. This contains information about all $`64`$ elements. The actual matrix is (for simplicity in notation, I will assume that all the $`f`$ functions are equal)

$$\rho _{123}=\left(\begin{array}{cccccccc}13f^2+2f^3& 0& 0& 0& 0& 0& 0& 0\\ 0& 1f^2& f^2& 0& f^2+f^3& 0& 0& 0\\ 0& f^2& 1f^2& 0& f^2+f^3& 0& 0& 0\\ 0& 0& 0& 1f^2& 0& f^2+f^3& f^2+f^3& 0\\ 0& f^2+f^3& f^2+f^3& 0& 1f^2& 0& 0& 0\\ 0& 0& 0& f^2+f^3& 0& 1f^2& 0& 0\\ 0& 0& 0& f^2+f^3& 0& 0& 1f^2& 0\\ 0& 0& 0& 0& 0& 0& 0& 13f^2+2f^3\end{array}\right),$$

In order to confirm that we obtain the same $`2`$ electron density matrix, we have to trace one electron out, and this we have to do by tracing over both the spin and space. If this is performed, we can easily confirm that the resulting density matrix is the one obtained previously for two electrons.

Now, we would like to know if this state contains three party entanglement in spin. We already know that every two parties can be entangled in the appropriate distance region, for example, when they are very close to each other. What happens when $`r=r^{}=r^{\prime \prime }`$? Then, $`f_i=1`$ and we see that the resulting state is actually not a physical state (all the elements go to zero indicating that this is not a physical state - or that the probability for this to happen vanishes). This is physically clear, as when we have three spins, at least two out of three have to be in the same state and this is not possible as they are all in the same physical location. At the other extreme, when all three electrons are very far away, then the density matrix again becomes maximally mixed. There is no entanglement in this case of any kind, nor do we have classical correlations. We also have a number of other possibilities, such as, for example, that two electrons are close to each other, but the third one is further away from them. We could then imagine that all three are in an entangled state as well as having entanglement between the two close by electrons. Let us take as an example the state when all electrons are very close to each other. Then, to a good approximation, the state is an equal mixture of the four state $`|000`$, $`|111`$, $`|001+|010+|100`$ and $`|011+|110+|101`$. Now, this state clearly contains three party entanglement. This is because every electron is entangled to the other two electrons, indicating that the state cannot be separable in any way. Following this logic therefore we have a three party entanglement providing $`f^2>1/2`$. Can we have a genuine three part entanglement beyond this point even though there is no two party entanglement? This would be something like a GHZ state (where there is no two party entanglement, but all three parties are entangled) and this remains an open question. Generalizing our line of reasoning, however, we can state that if a bunch of $`n`$ electrons is found within the distance of roughly the inverse Fermi wavenumber, then the corresponding total state of their spins would be an entangled state on all levels (two, three-party and so on up to $`n`$-party entanglement). The above treatment also allows us to compute various mutual informations between different parties involved.

Interestingly, an equivalent analysis can be performed for a bunch of non-interacting bosons (with polarization as the internal degree, for example), but we will not see any correlations in polarization. This is because there is no restriction on how many bosons can be in a state with a given momentum and polarization. The most general state of bosons is now given by

$$|\varphi \underset{sp}{}(b_s^{}(p))^{n(p,s)}|0$$
(59)
Where $`n(p,s)`$ is the number of bosons of momentum $`p`$ and polarization $`s`$. The operators obey the following commutation relations

$$[b_s^{}(p),b_t(q)]_{}=\delta _{st}\delta (pq)$$
(60)
It is now intuitive that if we measure a photon of one polarization at $`r`$, there is no reason why this should be correlated to the polarization of another photon at $`r^{}`$. We can compute the same density matrix as in the case of fermions, but the quantity $`\varphi _0|b_t^{}^{}(p)b_t^{}(q)b_s^{}(q^{})b_s(p^{})|\varphi _0`$ will now be independent of the labels $`s,s^{},t,t^{}`$. This means that the resulting $`4`$ by $`4`$ density matrix of polarization will have all the $`16`$ elements equal to each other. This represents a disentangled pure state $`(|H+|V)(|H+|V)`$, where $`H`$ and $`V`$ are the horizontal and vertical polarizations respectively (the polarization plays the same role here as the spin did in the case of fermions). And, so, here we see a fundamental difference between fermions and bosons with respect to the interplay between the entanglement in internal and spatial degrees of freedom. Therefore, there is no correlation between the relative distance between the photons and their polarization. Bosons, in other words, do not obey the Pauli exclusion principle.

Note that in this section we have analyzed the presence of entanglement in the internal degrees of freedom and have not talked about the spatial degrees of freedom. The analysis of long range correlations in fermions and bosons can be found in , although the author there does not distinguish between classical and quantum correlations (and the internal degrees of freedom play no special part). Nevertheless, similar methods could possibly be used in the future to perform a more detailed analysis of more complicated systems. Note also that all our discussions here could be generalized to overall mixed state (here the total states we start with are pure states) as in the statistical mechanics of many body systems . The same analysis can be performed as the one presented above to calculate the amount of entanglement in the internal degrees of freedom. It is to be expected, however, that this additional mixedness can only reduce the existing correlations in an overall pure state computed here. Otherwise, the results would not be fundamentally different to the ones presented here.

## VI Concluding Remarks

We have seen that the second quantization formalism adds some interesting effects to our understanding of entanglement. First and foremost, this is because we have to decide what our subsystems are before we can speak about the amount of entanglement between them. This is not an issue in the first quantized version of quantum mechanics that is used when we deal with a small number of systems, such as $`2`$ or $`3`$ (by definition) distinguishable) qubits. An interesting observation is that the change of modes also affects the amount of entanglement present. A change of modes is, of course, a “global change” relating the old modes to the new ones and this is why entanglement can be affected in the first place - “local changes” on their own would not be able to generate entanglement (We have to be a bit more careful here as we have to restrict ourselves to speaking about modes and excitations of the modes. We cannot really speak about local operations on “particles”, since they are identical and whatever operation we perform on one of them, we, in fact, have to perform on all of them). The change of modes produces an interaction between the new modes which is what gives rise to entanglement. In section $`3`$, we presented the most general formalism for mode changes and showed how the amount of entanglement varies with the change of modes. However, we have also been able to show that we do not really need an interaction to produce entanglement. We can instead rely on the effects of quantum particle statistics to produce entanglement. This effect was illustrated by looking at the Fermi sea model of a conductor, where we had a collection of non-interacting electrons (neither in spin nor in space) occupying states up to the highest Fermi level (assumed arbitrary). Here we found that two electrons are entangled in the spin degrees of freedom providing that their distance is (roughly) below the inverse Fermi momentum (measured in units of $`\mathrm{}=1`$). Our method can be used to compute the density matrix for any number of subsystems. We have also pointed out a crucial difference between non-interacting fermions and non-interacting bosons, in that non-interacting fermions can still be found to be entangled in the internal degrees of freedom, while bosons cannot. It would be interesting to investigate the further possibility of using these “non-interacting” correlations to perform information processing . It would also be interesting to extend our discussion to more complicated systems composed of fermions and bosons.

Acknowledgments. I would like to thank A. J. Fisher, J. R. Gittings, N. Paunković, J. Pachos and G. Castagnoli on many interesting and helpful discussions on the nature of the link between entanglement and particle statistics. S. J. van Enk and J. Schliemann are acknowledged for pointing out some related references to this work. My research on this subject is supported by European Community, Engineering and Physical Sciences Research Council, Hewlett-Packard and Elsag spa.


format: text

text:
# Transport phenomena and microscopic structure in partially miscible binary fluids: A simulation study of the symmetrical Lennard-Jones mixture

## 1 Introduction

While the interplay between static structure and dynamical properties of simple fluids is a problem intensively studied both by analytic theory  and simulations , and also extensions to rather complex fluids such as polymers are already given , somewhat less attention has been paid to simple liquid mixtures: while static properties of liquid-liquid miscibility have been analyzed in some detail both by analytic theory (see e.g. ) and simulations (see e.g. ), dynamic properties of fluid mixtures have been investigated by computer simulation only occasionally (e.g. ). In most of these studies, the precise phase behavior of the used model is not known, and also very small system sizes were used, and hence the proper interpretation of the results is difficult.

The present work aims to contribute towards filling this gap by studying both phase behavior and static structure in conjunction with dynamic correlations and transport phenomena. We address the simplest possible case, a symmetric binary Lennard-Jones mixture with a miscibility gap. Extensions of our approach to asymmetric Lennard-Jones mixtures and models suitable for a realistic description of binary mixtures of liquid metals (such as ) will be given in future work. We note that a specific asymmetric binary Lennard-Jones mixture with 80% A and 20% B particles is a very popular model for the study of the dynamics of glass-forming fluids . However, the static phase diagram of that model system is not yet known.

The investigation of a symmetric mixture is of particular interest, because one can clearly distinguish in this case long–ranged correlations in the vicinity of the critical point from the different local ordering of A- and B- particles that drives the unmixing transition well below the critical temperature $`T_\mathrm{c}`$: At the critical point a mixture of 50% A- and 50% B-particles is formed, i.e. $`x_\mathrm{B}=1x_\mathrm{A}=N_\mathrm{B}/(N_\mathrm{A}+N_\mathrm{B})=0.5`$ ($`N_\mathrm{A}`$ and $`N_\mathrm{B}`$ being the number of A and B particles in the simulation box, respectively), and due to the symmetry of the model static A-A and B-B correlations are identical in the latter case. Well below $`T_\mathrm{c}`$ at coexistence one of the species, say B, is the minority species and a B particle prefers a larger number of neighboring B particles than in the one-phase region at the same concentration of B particles. We shall discuss below how this interplay of long-ranged correlations in the vicinity of $`T_\mathrm{c}`$ and local ordering far below $`T_\mathrm{c}`$ along the coexistence states affects the dynamics of the mixture.

For a symmetric mixture the phase diagram and other static properties are most easily computed with Monte Carlo (MC) methods using the semi-grandcanonical ensemble rather than Molecular Dynamics (MD) methods. This type of MC method has been used with great success to establish binary alloy phase diagrams in the solid phase, both for lattice  and off-lattice  models. It was also used to study lattice models of polymers , and, in conjunction with finite size scaling methods, to investigate the critical behavior of symmetric binary Lennard-Jones fluids . It is also suitable for asymmetric models of fluid mixtures . In this work we combine semi-grandcanonical MC simulations with microcanonical MD runs (i.e. at constant energy $`E`$ and constant volume $`V`$). The former simulations yield well-equilibrated configurations at concentrations $`x_\mathrm{B}=1x_\mathrm{A}`$ and temperatures $`T`$ along the two-phase coexistence curve of the system and thus, these configurations can be used as initial states for the MD to study the structure and dynamics at coexistence and also in the vicinity of coexistence by changing the temperature or the composition.

The rest of the Paper is organized as follows. In the next section we shall briefly comment on the model and the used simulation techniques. A precise definition of the quantities that are computed is given in Sec. 3. Then we present the results for the structure (Sec. 4) and dynamics (Sec. 5) of our model that have been obtained from MD runs. Finally, Sec. 5 presents our conclusions and gives an outlook on further work.

## 2 Model and simulation methods

We consider a binary mixture of point-like particles interacting with Lennard-Jones potentials ($`\alpha ,\beta `$ A,B),

$$u(r_{ij})=4\epsilon _{\alpha \beta }\left[\left(\frac{\sigma _{\alpha \beta }}{r_{ij}}\right)^{12}\left(\frac{\sigma _{\alpha \beta }}{r_{ij}}\right)^6\right],r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|,$$
(1)
where $`\{\stackrel{}{r}_i\}`$ are the positions of the particles, and the Lennard-Jones parameters $`\epsilon _{\alpha \beta }`$ are chosen as follows

$$\sigma _{\mathrm{AA}}=\sigma _{\mathrm{BB}}=\sigma _{\mathrm{AB}}=\sigma ,\epsilon _{\mathrm{AA}}=\epsilon _{\mathrm{BB}}=\epsilon ,\delta \epsilon _{\mathrm{AB}}/\epsilon .$$
(2)
We measure all lengths in units of $`\sigma 1`$, and we also choose Boltzmann’s constant $`k_B1`$. The temperature $`T`$ will be measured in units of $`\epsilon 1`$, and (in the MD part) the masses of the two particle species are also chosen the same, $`m_\mathrm{A}=m_\mathrm{B}=1`$. The potential is truncated and shifted at $`r_{ij}=2.5\sigma `$.

The same model has been studied earlier by Wilding  (note that in this work the potential was also truncated at $`r_{ij}=2.5\sigma `$ but not shifted, and thus minor differences in the properties between our model and that of Ref.  must be expected). In the latter work, also the density of the binary fluid was varied over a wide range, $`0\rho 0.7`$, and it was found for $`\delta =0.7`$, that the $`\lambda `$-line of critical points for fluid-fluid phase separation ends at a critical end point at a density of about $`\rho \sigma ^30.6`$ at the coexistence curve of liquid-gas phase separation. At sufficiently smaller values of $`\delta `$ than 0.7 the $`\lambda `$ line moves towards the liquid-gas critical point where it merges into a tricritical point and one may expect that this happens at densities below $`\rho \sigma ^30.6`$ . In the present work, we are not interested in the liquid-gas transition of the model at all, and hence work at a much higher density, namely $`\rho =1.0`$. For this density we stay in the fluid phase for the whole temperature range of interest, i.e. $`1T1.8`$. Choosing $`\delta =0.5`$, the critical temperature of unmixing is $`T_c1.638\pm 0.005`$ at this density, see Fig. 2 below.

The equilibration of the model is done as follows. We start from random configurations with an equal number of A and B particles. The particles are then displaced by means of a standard Monte Carlo scheme in the NVT ensemble, i.e. at constant total particle number $`N`$, constant volume $`V`$, and constant temperature $`T`$, rejecting or accepting a trial move according to the standard Metropolis criterion. Thereby, the length of a trial displacement is a randomly chosen number in the interval $`[\sigma /20,\sigma /20]`$. After 100000 MC steps the semi-grandcanonical MC is switched on: At the end of each displacement step an attempted identity switch (B $``$ A or A $``$ B) of $`N/10`$ randomly chosen particles is made, which is again accepted or rejected according to a standard Metropolis criterion, where both the energy change $`\mathrm{\Delta }E`$ and the change of the chemical potential $`\pm \mathrm{\Delta }\mu `$ needs to be taken into account in the Boltzmann factor . Each configuration after this Monte Carlo move (whether the change was accepted or not) is counted as a state for the Monte Carlo averaging in the semi-grandcanonical ensemble, which is defined by the independent variables $`N`$, $`V`$, $`T`$, and the chemical potential differences $`\mathrm{\Delta }\mu `$. Note that the relative concentration $`x_\mathrm{B}=N_\mathrm{B}/(N_\mathrm{A}+N_\mathrm{B})=N_\mathrm{B}/N`$ is a fluctuating quantity and its average value is an output rather than an input to the simulation. Only for the choice $`\mathrm{\Delta }\mu =0`$ the symmetry of the model between A and B dictates that for $`TT_c`$ we must have $`x_\mathrm{B}=1/2`$, while for $`T<T_c`$ we obtain for $`\mathrm{\Delta }\mu =0`$ states at the two-phase coexistence curve. This is illustrated in Fig. 1, where the concentration distribution function $`P(x_\mathrm{A})`$ is plotted for three temperatures. Note that we have done 5 independent runs, each with a length of 400000 MC steps in the semi-grandcanonical ensemble, to calculate the distribution functions in Fig. 1 whereby we started the averaging after 100000 MC steps in each run. In principle, for a finite system there is no symmetry breaking possible in full thermal equilibrium, and hence one should observe always two peaks, one at $`x_\mathrm{A}^{\text{coex}(2)}`$ and the other one at $`x_\mathrm{A}^{\text{coex}(1)}=1x_\mathrm{A}^{\text{coex}(2)}`$, corresponding to the two branches of the two-phase coexistence curve. In practice, these two branches are separated at low temperatures by a huge free energy barrier in phase space. Therefore, as can be seen in Fig. 1, at $`T=1.2`$ and $`T=1.4`$ only the A-rich phase is observed, and no transition to the B-rich side of the coexistence curve has happened yet at all. Only for $`T=1.6`$, i.e. rather close to $`T_c`$ (which in Fig. 1 is evident from the large width of the peaks), transitions back and forth between the two coexisting phases have occurred (although there were not enough of those transitions to make $`P(x_\mathrm{A})`$ strictly symmetric around $`x_\mathrm{A}^{\text{crit}}=1/2`$).

This “ergodicity breaking” due to finite observation time in the context of first-order phase transitions is a well-known phenomenon  and does not hamper our analysis at all. From data such as shown in Fig. 1 we hence obtain, when we normalize the part of $`P(x_\mathrm{A})`$ for $`x_\mathrm{A}>1/2`$ to $`P(x_\mathrm{A})𝑑x_\mathrm{A}=1`$, both estimates for the location of the coexistence curve and the “concentration susceptibility” $`\chi `$,

$$x_\mathrm{A}^{\text{coex}(2)}=\underset{1/2}{\overset{1}{}}x_\mathrm{A}P(x_\mathrm{A})dx_\mathrm{A},k_BT\chi =N(\underset{1/2}{\overset{1}{}}x_\mathrm{A}^2P(x_\mathrm{A})dx_\mathrm{A}[x_\mathrm{A}^{\text{coex}(2)}]^2).$$
(3)
In order to investigate the dynamics of well-equilibrated states precisely at the coexistence curve, we store a number of $`n=10`$ independent configurations at each of the 5 independent runs which have a number $`N_\mathrm{A}`$ of A particles such that the resulting ratio $`x_\mathrm{A}`$ is as close to $`x_\mathrm{A}^{\text{coex}(2)}`$ as possible (note that $`x_\mathrm{A}`$ is “quantized”, because $`N_A`$, $`N`$ are integers, but for large $`N`$ the discreteness of $`x_A`$, which is $`1/N`$, is smaller than the statistical error with which $`x_\mathrm{A}^{\text{coex}(2)}`$ can be estimated, and therefore this effect does not matter). With these configurations we can realize microcanonical MD runs right at the coexistence curve. Note that we have used in this work a standard MD code applying the velocity Verlet algorithm with a time step $`\delta t=0.01`$ {in units of the time $`t_0=(m\sigma ^2/48\epsilon )^{1/2}`$}.

When we work in the one phase region we just heat up the configurations that we have gotten from the semi-grandcanonical Monte-Carlo at coexistence and we equilibrate the system by MD at the desired temperature $`T_\mathrm{f}`$. The temperature is then kept constant by coupling the system to a stochastic heat bath, i.e. every 50 MD steps the velocities are replaced by new ones that are chosen randomly from a Boltzmann distribution corresponding to the temperature $`T_\mathrm{f}`$.

We have described these simple and almost straightforward considerations  in such detail here, because the literature is still full of examples where the authors are not aware of the pitfalls presented to equilibration by the presence of slowly relaxing variables. In the present problem, the conservation of the concentration (in the ensemble where both $`N_\mathrm{A}`$ and $`N_\mathrm{B}`$ are held fixed) together with the presence of a large correlation length of concentration fluctuations (which for $`N\mathrm{}`$ diverges at the critical point and causes the critical slowing down ) may lead to a very slow relaxation. Hydrodynamic slowing down is eliminated by the choice of an ensemble where the slow variable is not conserved, in our case the semi-grandcanonical ensemble. Critical slowing down cannot be avoided (apart from the special case of Ising models on lattices, for which cluster algorithms have been put forward  that eliminate also critical slowing down to a large extent). In the present work, we hence do not attempt to study the immediate vicinity of the critical point.

Fig. 2 shows the coexistence curve as estimated from our semi-grandcanonical method. One can see clearly that for $`T1.4`$ finite size effects on the coexistence curve are quite negligible, if $`N400`$. For $`T=1.5`$, data for $`N=400`$ overestimate the width of the two phase region very slightly, while data for $`N=`$ 800, 1600 and 3200 agree within statistical errors. For $`T=1.6`$, however, data for both $`N=400`$ and $`N=800`$ overestimate the width of the coexistence curve slightly, while the data for both $`N=1600`$ and $`N=3200`$ still agree. For this reason, we have chosen $`N=1600`$ as the standard particle number, for which all analysis of static and dynamic correlations have been made (Secs. 3, 4). We have estimated the critical temperature $`T_c`$ from power law fits of the form

$$f(x_\mathrm{B})=1/2\pm x_\mathrm{B}^{\text{coex}}=\widehat{B}(1T/T_c)^\beta $$
(4)
with the critical exponent $`\beta =0.33`$ which corresponds to the universality class of the three-dimensional Ising model. $`\widehat{B}`$ and $`T_c`$ have been used as adjustable parameters. For $`N=400`$ the fit yields $`T_c1.666\pm 0.005`$, while for all larger values of $`N`$ the result is $`T_c1.638\pm 0.005`$ (see Fig. 2). We emphasize at this point, however, that we deliberately do not address the problem of analyzing the critical behavior of transport coefficients in our study: for such a purpose one would need to work with system sizes $`N`$ that are several orders of magnitude larger, and extremely accurate data in the temperature range $`1.6T1.7`$ would be needed. This interesting problem would require substantially more computer resources than were available to us. In addition, we would have to locate $`T_c`$ much more accurately with a finite size scaling analysis . Even though we do not study critical behavior, our choice of $`N`$ is much larger than most of the choices made in previous work on the dynamics of mixtures (e.g. Vogelsang and Hoheisel  work with $`N=108`$ and $`N=256`$ particles, and most data of Asta et al.  refer to $`N=500`$ particles).

## 3 Static and dynamic quantities

We now proceed to define the quantities that are computed in our simulations. With respect to static quantities, we have monitored the standard pairwise radial distribution functions $`g_{\alpha \beta }(r)`$ between the different pairs $`\{\alpha ,\beta (`$A,B$`)\}`$ of particles ,

$$g_{\alpha \beta }(r)=\frac{N}{\rho N_\alpha N_\beta }\underset{i=1}{\overset{N_\alpha }{}}\underset{j=1}{\overset{N_\beta }{}}{}_{}{}^{}\delta (r|\stackrel{}{r}_i(0)\stackrel{}{r}_j(t)|).$$
(5)
The prime in the second sum means that $`i=j`$ has to be left out if $`\alpha =\beta `$.

In terms of the $`g_{\alpha \beta }(r)`$, partial structure functions $`S_{\alpha \beta }(q)`$ are defined as follows :

$$S_{\alpha \beta }(q)=x_\alpha \delta _{\alpha \beta }+x_\alpha x_\beta \rho \underset{0}{\overset{\mathrm{}}{}}g_{\alpha \beta }(r)\frac{\mathrm{sin}(qr)}{qr}\mathrm{\hspace{0.17em}\hspace{0.17em}4}\pi r^2𝑑r.$$
(6)
with $`x_\alpha =N_\alpha /N`$ ($`\alpha `$ A, B). In order to compute the partial structure factors we used the formula

$$S_{\alpha \beta }(\stackrel{}{q})=\frac{1}{N}\underset{k=1}{\overset{N_\alpha }{}}\underset{l=1}{\overset{N_\beta }{}}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r}_{kl}),$$
(7)
and performed the angular integration numerically to obtain $`S_{\alpha \beta }(q)`$.

It turns out that it is more useful to form the following linear combinations of the radial distribution functions $`g_{\alpha \beta }`$ and the partial structure factors $`S_{\alpha \beta }`$ which describe correlations in number density ($`g_{nn}`$ and $`S_{nn}`$), concentration correlations ($`g_{cc}`$ and $`S_{cc}`$), and cross correlations between number density and concentration ($`g_{nc}`$ and $`S_{nc}`$) . The real space functions are given by

$`g_{nn}(r)`$ $`=`$ $`x_\mathrm{A}^2g_{\mathrm{AA}}(r)+2x_\mathrm{A}x_\mathrm{B}g_{\mathrm{AB}}(q)+x_\mathrm{B}^2g_{\mathrm{BB}}(r),`$ (8)
$`g_{cc}(r)`$ $`=`$ $`x_\mathrm{A}^2x_\mathrm{B}^2\left[g_{\mathrm{AA}}(r)+g_{\mathrm{BB}}(r)2g_{\mathrm{AB}}(r)\right],`$ (9)
$`g_{nc}(r)`$ $`=`$ $`x_\mathrm{A}x_\mathrm{B}\left[x_\mathrm{A}g_{\mathrm{AA}}(r)x_\mathrm{B}g_{\mathrm{BB}}(r)+(x_\mathrm{B}x_\mathrm{A})g_{\mathrm{AB}}(r)\right],`$ (10)
and the functions in reciprocal space by

$`S_{nn}(q)`$ $`=`$ $`S_{\mathrm{AA}}(q)+2S_{\mathrm{AB}}(q)+S_{\mathrm{BB}}(q),`$ (11)
$`S_{cc}(q)`$ $`=`$ $`x_\mathrm{B}^2S_{\mathrm{AA}}(q)+x_\mathrm{A}^2S_{\mathrm{BB}}(q)2x_\mathrm{A}x_\mathrm{B}S_{\mathrm{AB}}(q),`$ (12)
$`S_{nc}(q)`$ $`=`$ $`x_\mathrm{B}S_{\mathrm{AA}}(q)x_\mathrm{A}S_{\mathrm{BB}}(q)+(x_\mathrm{B}x_\mathrm{A})S_{\mathrm{AB}}(q).`$ (13)
We will see below that strong concentration fluctuations are observed if one approaches the critical point of unmixing (and this happens not only in the immediate vicinity of the critical point). One can try to describe $`S_{cc}(q)`$ at small $`q`$ by the Ornstein-Zernike form,

$$S_{cc}(q)=\frac{k_BT\chi }{1+\xi ^2q^2},$$
(14)
where the “susceptibility” $`\chi `$ is given by Eq. (3) and $`\xi `$ has the meaning of a static correlation length.

Turning to dynamic quantities, the most straightforward quantity to consider is the mean square displacement of tagged particles

$$g_\mathrm{A}(t)=[\stackrel{}{r}_{i,\mathrm{A}}(0)\stackrel{}{r}_{i,\mathrm{A}}(t)]^2,g_\mathrm{B}(t)=[\stackrel{}{r}_{j,\mathrm{B}}(0)\stackrel{}{r}_{j,\mathrm{B}}(t)]^2,$$
(15)
where it is understood that the average $`\mathrm{}`$ includes an average over all particles of type A or B, respectively. Another quantity which also monitors the motion of individual particles are the incoherent intermediate structure functions $`F_s^{(\mathrm{A})}(q,t)`$ and $`F_s^{(\mathrm{B})}(q,t)`$,

$`F_s^{(\mathrm{A})}(q,t)`$ $`=`$ $`{\displaystyle \frac{1}{N_\mathrm{A}}}{\displaystyle \underset{i\mathrm{A}}{}}\mathrm{exp}\{i\stackrel{}{q}[\stackrel{}{r}_i(0)\stackrel{}{r}_i(t)]\},`$ (16)
$`F_s^{(\mathrm{B})}(q,t)`$ $`=`$ $`{\displaystyle \frac{1}{N_\mathrm{B}}}{\displaystyle \underset{j\mathrm{B}}{}}\mathrm{exp}\{i\stackrel{}{q}[\stackrel{}{r}_j(0)\stackrel{}{r}_j(t)]\},`$ (17)
where the sum runs over all A particles {Eq. (16)} or B particles {Eq. (17)}, respectively. Note that due to the isotropy of the fluid, Eqs. (16), (17) cannot depend on the direction of the scattering vector $`\stackrel{}{q}`$. When the displacements of the particles are Gaussian distributed, Eqs. (16), (17) reduce to

$$F_s^{(\mathrm{A})}(q,t)=\mathrm{exp}\left[\frac{1}{6}q^2g_\mathrm{A}(t)\right],F_s^{(\mathrm{B})}(q,t)=\mathrm{exp}\left[\frac{1}{6}q^2g_\mathrm{B}(t)\right].$$
(18)
Noting further the Einstein relations

$$g_\mathrm{A}(t)=6D_\mathrm{A}t,g_\mathrm{B}(t)=6D_\mathrm{B}t,t\mathrm{},$$
(19)
where $`D_\mathrm{A}`$ and $`D_\mathrm{B}`$ are self-diffusion constants of the particles, we expect that the asymptotic decay of $`F_s^{(\mathrm{A})}(q,t)`$ and $`F_s^{(\mathrm{B})}(q,t)`$ can be described by a simple exponential variation with time,

$$F_s^{(\mathrm{A})}(q,t)\mathrm{exp}[t/\tau _\mathrm{A}(q)],F_s^{(\mathrm{B})}(q,t)\mathrm{exp}[t/\tau _\mathrm{B}(q)],$$
(20)
and Eqs. (18)-(20) then imply

$$\underset{q0}{lim}[\tau _\mathrm{A}(q)q^2]^1=D_\mathrm{A},\underset{q0}{lim}[\tau _\mathrm{B}(q)q^2]^1=D_\mathrm{B}.$$
(21)
In addition, we introduce partial coherent intermediate structure functions $`S_{\alpha \beta }(q,t)`$ defined as $`[\alpha ,\beta `$ (A,B)\]

$$S_{\alpha \beta }(q,t)=\frac{1}{N}\underset{\genfrac{}{}{0pt}{}{i=1}{i\alpha }}{\overset{N_\alpha }{}}\underset{\genfrac{}{}{0pt}{}{j=1}{j\beta }}{\overset{N_\beta }{}}\mathrm{exp}\{i\stackrel{}{q}[\stackrel{}{r}_i(0)\stackrel{}{r}_j(t)]\},$$
(22)
and from these functions we can define number density and concentration correlation as well as cross correlation functions $`\{S_{nn}(q,t),S_{cc}(q,t)`$ and $`S_{nc}(q,t)\}`$ by forming exactly the same type of linear combinations as written in Eqs. (11)-(13), but with the $`S_{\alpha \beta }(q,t)`$ rather than their static counterpart $`S_{\alpha \beta }(q)`$. In order to compare the time correlation functions $`S_{\alpha \beta }(q,t)`$ for different values of $`q`$ it is convenient to consider the normalized functions

$$F_{\alpha \beta }(q)=\frac{S_{\alpha \beta }(q,t)}{S_{\alpha \beta }(q)}.$$
(23)
The normalized functions $`F_{nn}(q,t)`$, $`F_{cc}(q,t)`$, and $`F_{nc}(q,t)`$ are defined in a similar way by normalizing $`S_{nn}(q,t)`$, $`S_{cc}(q,t)`$ and $`S_{nc}(q,t)`$ by the corresponding static correlation function.

Further quantities of interest are the instantaneous values of the components of the pressure tensor ($`x,y,z`$ denote the Cartesian components, $`\stackrel{}{v}_i`$ is the velocity of the $`i`$’th particle)

$$\sigma _{xy}(t)=\underset{i=1}{\overset{N}{}}\left[m_iv_{ix}v_{iy}+\frac{1}{2}\underset{j(i)}{}|x_ix_j|F_y(|\stackrel{}{r}_i\stackrel{}{r}_j|)\right],$$
(24)
$`\stackrel{}{F}`$ being the force acting between particles $`i,j`$. From the autocorrelation function of $`\sigma _{xy}(t)`$ one gets the shear viscosity $`\eta `$ of the fluid by a Green-Kubo formula ,

$$\eta =\frac{1}{Vk_BT}\underset{0}{\overset{\mathrm{}}{}}𝑑t\sigma _{xy}(0)\sigma _{xy}(t).$$
(25)
Finally we mention the Green-Kubo formula for the interdiffusion constant (note that we have chosen the masses of all particles equal to each other here)

$$D_{\text{int}}=\frac{1}{NS_{cc}(0)}\underset{0}{\overset{\mathrm{}}{}}𝑑tJ_x^{\mathrm{int}}(0)J_x^{\mathrm{int}}(t),$$
(26)
where the current $`\stackrel{}{J}^{\mathrm{int}}(t)`$ is defined as follows:

$$\stackrel{}{J}^{\mathrm{int}}(t)=x_\mathrm{B}\underset{i=1}{\overset{N_\mathrm{A}}{}}\stackrel{}{v}_i(t)x_\mathrm{A}\underset{i=1}{\overset{N_\mathrm{B}}{}}\stackrel{}{v}_i(t).$$
(27)
One can easily calculate $`D_{\text{int}}`$ for an ideal mixture \[such a mixture is formed if the particles are labeled but are otherwise identical and thus the entropy of mixing is equal to $`Nk_B(x_\mathrm{A}\mathrm{log}x_\mathrm{A}+x_\mathrm{B}\mathrm{log}x_\mathrm{B})`$\]. In this case cross correlations in the autocorrelation function for $`J_x^{\mathrm{int}}`$ in Eq. (26) vanish and the interdiffusion constant can be expressed by the self-diffusion constants ,

$$D_{\text{int}}=x_\mathrm{B}D_\mathrm{A}+x_\mathrm{A}D_\mathrm{B}.$$
(28)
It has been shown by MD simulations that Eq. (28) is a good approximation for mixtures of Lennard-Jones fluids in the one-phase region . However, for ideal mixtures on a rigid lattice (where A, B-atoms hop with jump rates $`\sigma _\mathrm{A}`$, $`\sigma _\mathrm{B}`$ to vacant sites, assuming that a small number of vacancies is present) Eq. (28) does not hold .

## 4 Static properties of the symmetrical binary fluid mixture

In this section we analyze the structure of the symmetric Lennard-Jones mixture. Results are shown along the coexistence line and along the different paths in the one-phase region that are indicated in Fig. 2.

In Fig. 3 we present the pair correlation functions $`g_{\alpha \beta }(r)`$ for the three states at temperatures $`T=1.2`$, 1.4, and 1.6 at the coexistence curve (cf. Fig. 2). These data show a typical normal fluid behavior in all cases, as expected. Although it would be hard to recognize from these partial radial distribution functions that one approaches the critical point of unmixing in a binary fluid, the comparison of the different correlations shows nontrivial features. At each temperature, the first peak in $`g_{\mathrm{AB}}(r)`$ is at a slightly smaller distance than in $`g_{\mathrm{AA}}(r)`$ and $`g_{\mathrm{BB}}(r)`$. Moreover, the second peak in $`g_{\mathrm{BB}}(r)`$ is shifted towards smaller distances compared to the other two functions and this effect is more pronounced the lower the temperature is, i.e. the farther one is away from the critical point. The explanation of this feature is simple: The B particles are the minority species at the considered state points and thus, at small concentrations $`x_\mathrm{B}`$ it is very likely that one finds an A particle between two next-nearest B particles forming a B-A-B sequence. Since the distance between nearest A-B neighbors is smaller than that between nearest A-A and B-B neighbors one may expect that a large amount of such B-A-B sequences, where the distance between next-nearest B particles is smaller than in a B-B-B sequence, leads to the shift in the second peak of $`g_{\mathrm{BB}}(r)`$. Certainly, the difference between $`g_{\mathrm{AA}}`$ and $`g_{\mathrm{BB}}`$ becomes less pronounced the closer one approaches $`x_\mathrm{B}=0.5`$ (and thus the critical point) where, due to the symmetry of our model, $`g_{\mathrm{AA}}`$ and $`g_{\mathrm{BB}}`$ are identical.

The temperature dependence of $`g_{\mathrm{AA}}(r)`$ and $`g_{\mathrm{BB}}(r)`$ in the one-phase region is shown in Fig. 4 for the constant concentration $`x_\mathrm{B}=0.10375`$ in the temperature range $`1.4T1.8`$ ($`T=1.4`$ corresponds to a state on the coexistence curve). We see that the temperature dependence is very weak for $`g_{\mathrm{AA}}(r)`$ whereas we observe a significant increase of the amplitude of the first peak in $`g_{\mathrm{BB}}(r)`$ decreasing the temperature towards the coexistence line. The latter effect is even more pronounced in the coordination number function $`z_{\alpha \alpha }(r)`$ that is shown in the insets of Fig. 4. This function is defined by

$$z_{\alpha \alpha }(r)=x_\alpha _0^r𝑑r^{}\mathrm{\hspace{0.33em}4}\pi r^{\mathrm{\hspace{0.17em}2}}g_{\alpha \alpha }(r^{})\alpha [\mathrm{A},\mathrm{B}]$$
(29)
which gives the number of particles of type $`\alpha `$ surrounding a particle of type $`\alpha `$ within a distance $`r^{}<r`$. In the considered temperature range $`z_{\mathrm{BB}}(r)`$ changes at $`r=1.4`$ (i.e. around the location of the first minimum) from $`z_{\mathrm{BB}}=1.37`$ for $`T=1.8`$ to $`z_{\mathrm{BB}}=1.55`$ for $`T=1.4`$. By a closer inspection of $`z_{\mathrm{AA}}(r)`$ one finds that the changes in this quantity are of the same order. Thus there are only minor changes in the local order if one approaches the coexistence line from above.

We now proceed to discuss the behavior of the structure factors $`S_{nn}(q)`$, $`S_{cc}(q)`$, and $`S_{nc}(q)`$ {Eqs. (11)-(13)} along the coexistence line and for the latter two quantities also the corresponding functions in real space {Eqs. (9)-(10)}, see Fig. 5. The number density structure factor $`S_{nn}(q)`$ behaves exactly as expected for any dense simple liquid, and in particular $`S_{nn}(q0)`$ is very small, as expected for liquids that are almost incompressible. The height of the first peak is around 3, as is typical for fluids at temperatures somewhat higher than the melting temperature, and there is a pronounced second peak, due to the rather regular close packing of atoms in a dense fluid.

In contrast, $`S_{cc}(q)`$ at $`T=1.2`$ is very small and structureless throughout: in symmetrical binary mixtures, there is very little coupling between density and concentration fluctuations. For $`T=1.4`$, we see already some enhancement of $`S_{cc}(q)`$ at small $`q`$, reflecting the growing concentration fluctuations. For $`T=1.6`$, at small $`q`$ a rather dramatic growth of $`S_{cc}(q0)`$ signals the proximity of the critical point. For small $`q`$ we have fitted the data by the Ornstein-Zernike form Eq. (14) (bold lines in Fig. 5b), where we have used the correlation length $`\xi `$ as an adjustable parameter and we have estimated the susceptibility $`\chi `$ by means of Eq. (3) from the distribution function $`P(x_\mathrm{A})`$ (the values of $`S_{cc}(0)=k_BT\chi `$ are shown as crosses in Fig. 5b). From the fits we recognize that the data for small $`q`$ can be described by Eq. (14) within the statistical errors (note that there is only one adjustable parameter, the correlation length $`\xi `$). From the fits we find a growing correlation length $`\xi `$ with increasing temperature towards $`T_c`$. This behavior can be also inferred from $`g_{cc}(r)`$ (see inset of Fig. 5b) at $`T=1.6`$ which has a rather pronounced tail for $`r>2`$. Also remarkable in $`g_{cc}(r)`$ is the negative peak around $`r0.9`$ which stems from the fact that these distances are avoided by nearest A-A and B-B neighbors but are typical for nearest A-B neighbors.

Particularly interesting is the structure factor $`S_{nc}(q)`$. In a mixture with $`\mathrm{\Delta }\mu =0`$ at $`T>T_c`$, where then $`x_\mathrm{A}=x_\mathrm{B}=1/2`$, we expect that $`S_{nc}(q)0`$ due to the perfect symmetry between A and B. Along the coexistence curve, the spontaneous symmetry breaking between A and B invalidates this argument, and hence a nonzero $`S_{nc}(q)`$ is possible. As we see in Fig. 5c this is obviously the case. The origin of the oscillations in $`S_{nc}(q)`$ can be most easily understood by means of its counterpart in real space, $`g_{nc}(r)`$ (see inset of Fig. 5c). The most pronounced feature in the latter quantity is a peak with negative intensity around $`r_10.9`$. This feature is again due to the fact that the first peak in $`g_{\mathrm{AB}}(r)`$ is slightly shifted to smaller distances compared to the first peak in $`g_{\mathrm{AA}}`$ and in $`g_{\mathrm{BB}}`$ such that the intensity of the latter two functions in the interval $`0.85<r<0.9`$ is essentially zero whereas $`g_{\mathrm{AB}}`$ is rapidly increasing for the latter distances. This leads to the occurrence of negative correlations in $`g_{nc}`$ and thus to the peak around $`r_1`$. We emphasize that the scale for $`S_{nc}(q)`$ and $`g_{nc}(r)`$ is clearly very much smaller than for the other structure factors and pair correlation functions, respectively.

Similar data have been taken for two other paths in the phase diagram of Fig. 2, namely one at constant temperature, $`T=1.5`$, varying $`\mathrm{\Delta }\mu `$ (and hence $`x_\mathrm{B}`$), and one at constant $`x_B=0.10375`$. It turns out that the radial distribution functions $`g_{\alpha \beta }(r)`$ always look similar to those that are shown in Fig. 3, and hence are not shown here. The same statement applies to $`S_{nn}(q)`$, which exhibits only a weak dependence on either temperature or composition. Here we hence show only those quantities which have a more pronounced dependence on our control parameters, namely $`S_{cc}(q)`$ and $`S_{nc}(q)`$, Figs. (6), (7). From Fig. 6a we see that the enhancement of $`S_{cc}(q0)`$ persists to rather high temperatures, and also $`S_{nc}(q)`$ is only weakly temperature-dependent in the temperature region shown. In contrast, when we change the composition at $`T=1.5`$, we find that $`S_{cc}(q)`$ quickly loses all structure when $`x_\mathrm{B}0`$ (Fig. 7a), and a related trend is observed in $`S_{nc}(q)`$.

## 5 Dynamic properties of the symmetrical binary fluid mixture

The detailed description of the structure of our Lennard-Jones model that we have done so far is necessary to understand its dynamic properties. The next two sections are devoted to the results for time-dependent correlation functions and transport coefficients. In Sec. 5.1 we analyze the dynamic properties of our system along the fluid-fluid coexistence line whereas in Sec. 5.2 dynamic quantities are discussed that have been obtained along the path at $`x_\mathrm{B}=0.10375`$ which starts in the one-phase region and ends at the coexistence line at $`T=1.4`$ (see Fig. 2).

### 5.1 The dynamics at phase coexistence

We start by showing the incoherent intermediate structure function for the A particles, $`F_s^{(\mathrm{A})}(q,t)`$ {Eq. (16)}, see Fig. 8 (data for B-particles look very similar and therefore are not shown). Note that for $`q7`$, where the peaks of the static structure factors $`S_{\alpha \beta }(q)`$ and $`S_{nn}(q)`$ occur, $`F_s^\alpha (q,t)`$, $`\alpha `$ (A,B), decays on the time scale of about $`t=20`$ and this decay is certainly even faster for $`q=10`$. Moreover, the temperature dependence of the $`F_s^\alpha (q,t)`$ is weak in the temperature range under consideration. Such a short structural relaxation time as well as a weak temperature dependence indeed are expected for any ordinary fluid, which has no tendency to glass formation .

On the other hand, for small $`q`$ the decay of $`F_s^{(\mathrm{A})}(q,t)`$, $`F_s^{(\mathrm{B})}(q,t)`$ is much slower, and this slowing down simply reflects the diffusive behavior predicted in Eqs. (18)-(21). Fig. 9 shows that this interpretation in fact is nicely consistent with the data in that, in agreement with Eq. (30), the product $`\tau _\alpha q^2`$ approaches the inverse diffusion constant $`D_\alpha `$ for $`q0`$ (the values for $`D_\alpha ^1`$ are indicated as arrows in Fig. 9). One can further see that the B particles diffuse a bit faster than the A particles. Of course, we expect that this dynamic asymmetry must vanish when we are at the critical point, and indeed, as one can infer from Fig. 9, the ratio $`D_\mathrm{B}/D_\mathrm{A}`$ decreases to one if the critical point is approached from below along the coexistence line. The difference in the diffusion constants has a simple reason: In the A rich phase the B particles find mainly A particles as nearest neighbors, and since the interaction between particles of different species is weaker than for equal particles it is easier for the B particles to escape from an A rich neighborhood than from a B rich one. Thus, the effective activation energy for the diffusion of B particles is smaller in an A rich than in an B rich environment. Note that the behavior of $`\tau _\alpha q^2`$ at finite $`q`$ is very different from liquids near a glass transition: In such liquids one would see a peak in $`\tau _\alpha q^2`$ around a $`q`$ value that corresponds to the location of the structure factor maximum (see Ref.  and references therein).

Particularly interesting is the behavior of the collective structure functions (Figs. 1011). In $`F_{\mathrm{AA}}(q,t)`$ (Fig. 10) we observe again a very fast decay for $`q`$ values around the structure factor maximum $`q=7`$ and at $`q=0.992`$, i.e. at small $`q`$, the decay occurs on a two orders of magnitude larger time scale. For the latter value of $`q`$ one can clearly recognize that several relaxation processes contribute, since the decay emerges in two steps. Note that $`F_{\mathrm{AB}}(q,t)`$ and $`F_{\mathrm{BB}}(q,t)`$ exhibit a similar behavior and so are not shown here. In contrast to the partial structure functions, in the concentration-concentration correlation function $`F_{cc}(q,t)`$ only the slow decay is found at $`q=0.992`$ (Fig. 11). Surprisingly, this is also the case at $`T=1.2`$ although, at this temperature, there is almost no structure in the corresponding static function $`S_{cc}(q)`$. As can be also inferred from Fig. 11 the function $`F_{nn}(q,t)`$ shows a very fast decay to zero at $`q=0.992`$ (accompanied by oscillations due to acoustic modes), whereas it has a slightly slower relaxation to zero than $`F_{cc}(q,t)`$ at $`q=7`$, i.e. the location of the structure factor maximum.

The relaxation times for the collective correlation functions are of the same order as those detected from the self-correlation functions: disappointingly, there is no clear indication of critical slowing down when we approach $`T_c`$ from below, since in the temperature region studied the overall Arrhenius-like increase of the relaxation time as $`T`$ is lowered dominates (the Arrhenius-like behavior is demonstrated below in the case of the diffusion constants).

### 5.2 The dynamics from the one phase region to coexistence

The qualitative character of the relaxation functions $`F_s^{(\mathrm{A})}(q,t)`$, $`F_s^{(\mathrm{B})}(q,t)`$ and $`F_{\mathrm{AA}}(q,t)`$, $`F_{cc}(q,t)`$ in the one phase region is rather similar to the behavior found along the coexistence curve. Therefore, we do not show these functions in any detail here, but proceed immediately to a counterpart of Fig. 9 in Fig. 12, where now the relaxation times $`\tau _\mathrm{A}(q)q^2`$ and $`\tau _\mathrm{B}(q)q^2`$ extracted from the incoherent intermediate scattering functions at constant $`x_\mathrm{B}=0.10375`$ are shown as a function of wavenumber $`q`$. Again one finds that the minority species diffuses a bit faster than the majority. The reason for the larger discrepancies of $`D_\alpha ^1`$ and $`\tau _\alpha (q)q^2`$ for $`q0`$ in the case of the B particles is of course due to the fact that the statistics is worse for the B particles which is the minority species.

This difference in the diffusion constants is emphasized in the Arrhenius plot, Fig. 13, where the logarithm of $`D_\alpha `$ is plotted vs. inverse temperature. It is seen that a simple thermally activated behavior (which would show up as straight lines in this plot) holds only over very restricted temperature regimes which might be related to the fact that there are small but significant changes in the local structure when approaching the coexistence line (see especially Fig. 4). Also plotted in Fig. 13 are the inverse shear viscosity $`\eta ^1`$ and the interdiffusion constant $`D_{\mathrm{int}}`$. For an accurate computation of these quantities we have undertaken a much larger effort than for the single particle quantitities: At each temperature we averaged over 50 independent runs of 400000 MD steps to calculate $`\eta `$ and $`D_{\mathrm{int}}`$ by means of Eq. (25) and Eq. (26), respectively. We show also in Fig. 13 the interdiffusion constant according to Eq. (28) that is valid for an ideal mixture. Already for $`T<2.5`$ $`D_{\mathrm{int}}`$ decouples from the ideal mixture approximation and thus from the self-diffusion constants, and $`D_{\mathrm{int}}`$ becomes significantly smaller than the $`D_\alpha `$ if one approaches the coexistence curve.

It is now interesting to check for the Stokes-Einstein relation, which relates to the self-diffusion constant of a diffusing spherical particle of diameter $`d`$ in a fluid of viscosity $`\eta `$ as

$$D=\frac{k_BT}{2\pi \eta d}.$$
(30)
The factor 2 in the denominator of Eq. (30) corresponds to the assumption of slip boundary conditions on the surface of the diffusing particle (for stick boundary conditions a factor 3 instead of 2 appears) .

Using the data for $`\eta `$, $`D_\mathrm{A}`$, $`D_\mathrm{B}`$ and invoking Eq. (30) one obtains the corresponding effective Stokes-Einstein diameters $`d_\mathrm{A}`$, $`d_\mathrm{B}`$ shown in Fig. 14. The data for $`d_\alpha `$ show that despite the strong concentration fluctuations in the system the Stokes-Einstein relation is a good approximation. The finding that the values for $`d_\mathrm{A}`$ and $`d_\mathrm{B}`$ are different is reasonable since the distance for nearest A-B neighbors is slightly smaller than that for A-A or B-B neighbors and the B particles as the minority species are mostly surrounded by A particles (see Fig. 4); thus, the effective hydrodynamic diameter for the B particles is slightly smaller than that for the A particles. Note that with the assumption of slip boundary conditions in Eq. (30) we obtain reasonable values around 1 for $`d_\mathrm{A}`$ and $`d_\mathrm{B}`$ whereas for stick boundary conditions these values are rather small (around 0.6).

## 6 Concluding remarks

The present paper considers as a simple model system for a binary fluid mixture with a miscibility gap a dense Lennard-Jones mixture with fully symmetrical interactions, $`\sigma _{\mathrm{AA}}=\sigma _{\mathrm{BB}}=\sigma _{\mathrm{AB}}=\sigma `$, $`\epsilon _{\mathrm{AA}}=\epsilon _{\mathrm{BB}}=\epsilon =1`$, $`\delta =\epsilon _{AB}\epsilon =0.5`$, and also the masses of the two types of particles are chosen the same, $`m_A=m_B=m`$. The aim of the present work was to carry out a feasibility study, where by combination of the semi-grandcanonical Monte Carlo technique with Molecular Dynamics simulations both static properties, including the phase diagram, and dynamic correlation functions and associated transport coefficients are obtained simultaneously. It is shown that outside of the critical region of the mixture rather small system sizes (such as a total number of $`N=1600`$ particles) already suffice to obtain quantitatively reliable results, while in the critical region (which could not yet be studied here) much larger sizes clearly are mandatory. Such large sizes would also be needed for the very interesting problem of spinodal decomposition of the binary fluid mixture after a quench from the one-phase region into the unstable region of the phase diagram (Fig. 2). Since our study is the first study of a fluid mixture where both the phase diagram and the relevant transport coefficients have been determined simultaneously, such an extension to far from equilibrium dynamics clearly would be very interesting.

Our studies of transport coefficients are distinct from earlier work on similar models by the fact the we know precisely where in the phase diagram the considered state points are for which transport coefficients have been determined. This fact clearly facilitates the proper interpretation of the results. A particular interesting finding concerns the asymmetric composition in our otherwise fully symmetric model.

In further work we also plan to study systematically Onsager coefficients to test the concepts of phenomenological irreversible thermodynamics for the present model and to check “mixing rules” for the interdiffusion coefficient. Moreover, we will present extensions to more realistic models of real materials which always exhibit asymmetry in the interactions, to allow also a comparison to suitable experiments.

Acknowledgments: The present research was supported by the Deutsche Forschungsgemeinschaft (DFG) under grant N<sup>o</sup> Bi314/18 (SPP 1120). One of the authors (J. H.) acknowledges the support of the DFG under grant N<sup>o</sup> HO 2231/2-1.


text:
# 1 Contour map of 𝑧_𝐹 which satisfies Δ𝑡_{𝑧=3.62}=1.3 Gyr as a function of Ωₘ₀ and 𝐻₀.

Probing the star formation epoch with Wilkinson Microwave Anisotropy Probe (WMAP) data and high-redshift QSOs

Jun’Ichi Yokoyama Department of Earth and Space Science, Graduate School of Science,
Osaka University, Toyonaka 560-0043, Japan

## Abstract

Using the values of the cosmological parameters obtained from Wilkinson Microwave Anisotropy Probe (WMAP) and the cosmological chemical clock observed by high-redshift quasars, it is shown that star formation should have started by the redshift $`z_F>9.5`$, and that $`z_F`$ is sensitive only to $`\mathrm{\Omega }_{\mathrm{m0}}h^2`$ in flat $`\mathrm{\Lambda }`$CDM cosmology.

The first year data release of Wilkinson Microwave Anisotropy Probe (WMAP) has opened a new era of high precision cosmology (Bennet et al. 2003). As for the values of the cosmological parameters it has basically confirmed that of the concordance model with significantly smaller error bars than before. Together with the observation of nearly scale-invariant adiabatic fluctuations, the inflation-based flat $`\mathrm{\Lambda }`$CDM model can be regarded as the primary candidate of the standard evolution model of our Universe.

Ever since Hubble’s discovery of cosmic expansion (Hubble 1929), who concluded that the cosmic expansion rate was as large as $`H_0500`$km/s/Mpc, we encountered the cosmic age problem from time to time in the history of modern cosmology. The introduction of vacuum-like energy density with negative pressure such as cosmological constant $`\mathrm{\Lambda }`$ as inferred from the magnitude-redshift relation of high-redshift type Ia supernovae (Perlmutter et al. 1999; Riess et al. 1998) has reconciled a relatively large Hubble parameter $`h=H_0/(100`$km/s/Mpc)$`0.7`$ (Freedmann et al. 2001) with a fairly large cosmic age $`t_013`$Gyr. Indeed the cosmic age obtained by WMAP observations supplemented by other cosmological observations, $`t_0=13.7\pm 0.2`$Gyr with $`h=0.71_{0.03}^{+0.04}`$, $`\mathrm{\Omega }_{\mathrm{m0}}h^2=0.135_{0.009}^{+0.008}`$, $`\mathrm{\Omega }_{\mathrm{v0}}=\mathrm{\Omega }_{\mathrm{tot0}}\mathrm{\Omega }_{\mathrm{m0}}`$, and $`\mathrm{\Omega }_{\mathrm{tot0}}=1.02\pm 0.02`$ are in good agreement with astrophysical estimate of globular cluster age and cosmological nuclear chronology (Spergel et al. 2003). Here $`\mathrm{\Omega }_{\mathrm{v0}}`$ denotes the current value of the vacuum-like energy density or dark energy in unit of the critical density.

Thus as far as the current Universe is concerned we have indeed the concordance of the cosmic age and various cosmological parameters at hand. However, the situation may change drasitically if we turn our attention to the high-redshift Universe. Indeed with the improved accuracy of current values of the cosmological parameters, we can calculate the cosmic age at redshift $`z`$, $`t(z)`$, more accurately. Comparing $`t(z)`$ with other measures of cosmic age at high-redshift we can obtain important cosmological and astrophysical information and constraints.

As for the latter we consider observation of Fe/Mg abundance ratio probed by high-redshift quasars. This quantity can serve as a cosmic chemical clock due to the difference of the origin (Hamann & Ferland 1993). That is, Fe is predominantly produced by type Ia supernovae whose precursor has a lifetime of $``$ 1 Gyr or larger, while Mg and other alpha elements are produced from type II supernovae with much shorter time scale. Thus the ratio Fe/Mg at the redshift $`z`$ indicates the time elapsed between the star formation epoch $`z_F`$ and $`z`$.

Previously this measure was mainly used to constrain the values of cosmological parameters with particular emphasis on whether cosmological constant is nonvanishing (Yoshii, Tsujimoto, & Kawara 1998). But at the new era of precision cosmology with WMAP and other sophisticated observational tools we can extract useful information on $`z_F`$ and chemical evolution models of the Universe as argued below.

Although the observational constraint of WMAP on the spatial curvature or the total energy density of the Universe still has a two percent uncertainty, $`\mathrm{\Omega }_{\mathrm{tot0}}=1.02\pm 0.02`$, we assume $`\mathrm{\Omega }_{\mathrm{tot0}}`$ is equal to unity with much higher accuracy, because we have good reasons to believe in standard inflation as argued above, which predicts $`|\mathrm{\Omega }_{\mathrm{tot0}}1|<Q5\times 10^6`$ with $`Q`$ being the quadrupole anisotropy of cosmic microwave background radiation (Kashlinsky, Tkachev & Frieman 1994). Then, in the spatially flat Robsertson-Walker spacetime with the scale factor $`a(t)`$, cosmic age at redshift $`z`$ is given by

$`t\left(z\right)`$ $`=`$ $`{\displaystyle \frac{da}{\dot{a}}}={\displaystyle \frac{1}{H_0}}{\displaystyle _z^{\mathrm{}}}{\displaystyle \frac{dz}{(1+z)E(z)}},`$ (1)
$`E(z)`$ $`=`$ $`\left[\mathrm{\Omega }_{\mathrm{m0}}(1+z)^3+\mathrm{\Omega }_{\mathrm{v0}}f(z)\right]^{1\mathrm{/}2},`$ (2)
where $`f(z)`$ is a function which specifies the evolution of dark energy component. If it has a constant index of equation of state, $`w_v\rho _v/p_v`$, where $`\rho `$ and $`p`$ denote energy density and pressure respectively, one finds $`f(z)=(1+z)^{3(1+w_v)}`$. The pure cosmological constant corresponds to $`w_v=1`$. Below we mainly consider this case because it is theoretically well-motivated (Yokoyama 2002a,b) and observationally supported among many other candidates of dark energy (Yokoyama 2003). In the above expression we have neglected radiation component, which is a good approximation in the range of redshifts we consider.

Several authors have calculated time evolution of Fe/Mg ratio using a chemical evolution model. For example, Yoshii, Tsujimoto & Kawara (1998) calculated its evolution based on the model of Yoshii, Tsujimoto & Nomoto (1996) with the following three basic ingredients. The first one is the star formation rate $`C(t)`$ parametrized as $`C(t)=\nu _k\left[f_g(t)\right]^k`$ where $`f_g(t)`$ is the gas fraction and $`\nu _k`$ and $`k`$ are parameters. As the standard model they choose $`k=1`$ and $`\nu _{k=1}=7.6\mathrm{Gyr}^1`$ as inferred from metal abundance in quasar host galaxies (Hamann & Ferland 1993). The second is the initial stellar mass function (IMF) which is assumed to be time-invariant with a spectrum $`\varphi (m)dmm^xdm`$ for $`0.05M_{}m50M_{}`$ with $`x=1.35`$ based on Tsujimoto et al. (1997). The third and the most important ingredient is the fraction, $`A`$, and the lifetime, $`t_{Ia}`$, of the progenitors that eventually produce SNe Ia, and their time spread. They adopt $`A=0.055`$ for $`m=38M_{}`$ and $`A=0`$ outside this mass range as a standard value. $`t_{Ia}`$ is obtained from the break in \[O/Fe\] observed at \[Fe/H\]$`1`$ in the solar neighborhood (Yoshii, Tsujimoto & Nomoto 1996) and taken as $`t_{Ia}=1.5`$ Gyr, while its spread function, $`g(t_{Ia})`$, is modeled as a power-law $`g(t_{Ia})t_{Ia}^\gamma `$. In practice, they take $`\gamma =0`$ and $`g(t_{Ia})0`$ for $`t_{Ia}=13`$ Gyr, again to match the observational features of \[O/Fe\] and \[Fe/H\].

As for the observation, Yoshii, Tsujimoto and Kawara used UV-optical spectrum of QSO B1422+231 at $`z=3.62`$ and found FeII(UV+opt)/MgII $`\lambda \lambda `$ 2796, 2804 flux ratio to be $`12.2\pm 3.9`$. Assuming this ratio traces the abundance ratio Fe/Mg, they obtained \[Mg/Fe\]$`=0.61_{0.300.12}^{+0.12+0.16}`$, which, according to their standard evolution model, implies that the time elapsed from $`z_F`$ to $`z=3.62`$ is given by

$$\mathrm{\Delta }t_{z=3.62}t(z=3.62)t(z_F)=1.50\mathrm{Gyr}.$$
(3)
They have also considered possible modification of their standard model and concluded that the lower bound is 1.3 Gyr, namely, $`\mathrm{\Delta }t_{z=3.62}1.3`$ Gyr, which may be realized for a higher value of $`A0.08`$. As is seen in Table 1 of Yoshii, Tsujimoto and Kawara (1998), modification of other model parameters tends to increase $`\mathrm{\Delta }t`$.

This quasar is by no means the only QSO whose flux ratio FeII/MgII has been observed. Thompson, Hall and Elston (1999) analyzed spectra of 12 high-redshift quasars to probe evolution of FeII/MgII ratio. As a result they found no decline in the iron abundance even at $`z=4.47`$ and concluded that about 1Gyr had passed from $`z_F`$ to $`z=4.47`$. Dietrich et al. (2002), on the other hand, observed six high-redshift QSOs with $`z=33103.488`$ and basically confirmed the conclusion of Yoshii, Tsujimoto and Kawara (1998), namely, $`\mathrm{\Delta }t_{z3.4}1.5`$ Gyr.

Comparing these values with the formula (2) calculated with the improved values of cosmological parameters obtained with WMAP, we can find constraint on the epoch of star formation $`z_F`$. Specifically we consider the three possible tests on $`z_F`$ as argued above, to draw contours of $`z_F`$ corresponding to $`\mathrm{\Delta }t_{z=3.62}=1.3`$ Gyr, $`\mathrm{\Delta }t_{z=3.62}=1.5`$ Gyr and $`\mathrm{\Delta }t_{z=4.77}=1.0`$ Gyr.

Figures 1-3 depict these contour on an $`\mathrm{\Omega }_{\mathrm{m0}}H_0`$ plane. As is seen there the condition $`\mathrm{\Delta }t_{z=3.62}=1.5`$ Gyr is more stringent than that $`\mathrm{\Delta }t_{z=4.77}=1.0`$ Gyr. For the most likely value of the cosmological parameters obtained by WMAP, $`h=0.71`$ and $`\mathrm{\Omega }_{\mathrm{m0}}=0.268`$ under the assumption that $`\mathrm{\Omega }_{\mathrm{v0}}=1\mathrm{\Omega }_{\mathrm{m0}}`$, we find $`z_F=10.13`$ from $`\mathrm{\Delta }t_{z=3.62}=1.3`$ Gyr, $`z_F=14.98`$ from $`\mathrm{\Delta }t_{z=3.62}=1.5`$ Gyr, and $`z_F=11.91`$ from $`\mathrm{\Delta }t_{z=4.77}=1.0`$ Gyr. Thus significant star formation before $`z_F>10`$ is inferred. In principle we could yield even more stringent constraint than $`z_F=14.98`$, if we applied Yoshii et al’s model (1998) to Thompson et al’s data (1999), but we remain conservative here, allowing the room for uncertainty in chemical evolution models.

More interesting are contour maps depicted on a $`\mathrm{\Omega }_{\mathrm{m0}}h^2H_0`$ plane, which are shown in figs. 4-6. As is seen there in the range of cosmological parameters under consideration $`z_F`$ is insensitive to $`H_0`$ and depends only on $`\mathrm{\Omega }_{\mathrm{m0}}h^2`$. Thus for $`\mathrm{\Omega }_{\mathrm{m0}}h^2=0.135_{0.009}^{+0.008}`$, we obtain $`z_F=10.1_{0.6}^{+0.7}`$ from $`\mathrm{\Delta }t_{z=3.62}=1.3`$ Gyr, $`z_F=15.0_{1.7}^{+1.9}`$ from $`\mathrm{\Delta }t_{z=3.62}=1.5`$ Gyr, and $`z_F=11.9\pm 0.7`$ from $`\mathrm{\Delta }t_{z=4.77}=1.0`$ Gyr. We may therefore conclude $`z_F>9.5`$ in the flat $`\mathrm{\Lambda }`$CDM Universe with the cosmological parameters obtained by combination of WMAP and other measures.

Finally we briefly mention the case dark energy has a softer equation of state $`w>1`$ than the cosmological constant. For a constant coefficient $`w`$, WMAP has obtained a constraint $`w<0.78`$ (Spergel et al. 2003). Dark energy with a softer equation of state tends to require star formation at higher redshift, although the change is modest. For example, for $`w=0.8`$, $`h=0.71`$, and $`\mathrm{\Omega }_{\mathrm{m0}}=0.268`$ we find $`z_F=10.36,`$ 15.57, and 12.10, respectively instead of $`z_F=10.13,`$ 14.98, and 11.91 for the case of the cosmological constant as discussed above.

In conclusion we have pointed out that in the new era of precision cosmology we can obtain some useful information on the evolution of the high-redshift Universe from the cosmological chemical clock rather than constraining current values of cosmological parameters with them. We have shown that the epoch of star formation $`z_F`$ should satisfy $`z_F>9.5`$ and that it is sensitive only to $`\mathrm{\Omega }_{\mathrm{m0}}h^2`$ in the context of the standard flat $`\mathrm{\Lambda }`$CDM cosmology. This constraint would bring about interesting implication to formation of galaxies that host quasars.

This work was partially supported by the JSPS Grant-in-Aid for Scientific Research No. 13640285.

Bennett, C.L., Halpern, M., Hinshaw, G., Jarosik, N., Kogut, A., Limon, M., Meyer, S.S., Page, L., Spergel, D.N. et al., 2003, astro-ph/0302207

Dietrich, M., Appenzeller, I., Vestergaard, M., and Wagner, S.J., 2002, ApJ, 564, 581

Freedman, W.L., Madore, B.F., Gibson, B.K., Ferrarese, L., Kelson, D.D., Sakai, S., Mould, J.R., Kennicutt, R.C. et al., 2001, ApJ, 553, 47

Hamann, F. and Ferland, G., 1993, ApJ, 418, 11

Hubble, E.P, 1929, Proc. Nat. Acad. Sci, 15, 168

Kashlinsky, A., Tkachev, I.I., Frieman, J., 1994, Phys. Rev. Lett., 73, 1582

Perlmutter, S., Aldering, G., Goldhaber, G., Knop, R.A., Nugent, P., Castro, P.G., Deustua, S., Fabbro, S., Goobar, A., et al., 1999, ApJ, 517, 565

Riess, A.G., Nugent, P.E., Gilliland, R.L., Schmidt, B.P., Tonry, J., Dickinson, M., Thompson, R.I., Budavari, T.S., et al., 1998 AJ, 116, 1009

Spergel, D.N., Verde, L., Peiris, H.V., Komatsu, E., Nolta, M.R., Bennett, C.L., Halpern, M., Hinshaw, G., et al., 2003, astro-ph/0302209

Thompson, K.L., Hill, G.J., Elston, R., 1999, ApJ, 515, 487

Tsujimoto, T., Yoshii, Y., Nomoto, K., Matteucci, F., Thielemann, F.-K., & Hashimoto, M., 1997, ApJ, 483, 228

Yoshii, Y., Tsujimoto, T., and Nomoto, K., 1996, ApJ, 462, 266

Yoshii, Y., Tsujimoto, T., and Kawara, K., 1998, ApJL, 507, L113

Yokoyama, J. 2002a, Phys. Rev. Lett., 88, 151302

Yokoyama, J. 2002b, Int. J. Mod. Phys., D11, 1603

Yokoyama, J. 2003, In Proc. 12th workshop on general relativity and gravitation, eds. Y. Eriguchi and M. Shibata. (University of Tokyo)