text:
# Scaling solutions on a brane

## 1 Introduction

Cosmological evolution of a brane world has become a matter of intense investigation after the famous paper of L. Randall and R. Sundrum . Due to additional terms, which appear in effective Einstein equations on a brane, the brane cosmology has some important differences from the standard cosmological scenario.

The field equations on the brane are

$$G_{\mu \nu }=\kappa ^2T_{\mu \nu }+\stackrel{~}{\kappa }^4S_{\mu \nu }_{\mu \nu }.$$
(1)
Here $`\stackrel{~}{\kappa }=m_5^{3/2}`$ is the fundamental 5-dimensional gravitational constant, $`\kappa =m_4^1`$ is the effective 4-dimensional gravitation constant on the brane. Bulk corrections to the Einstein equations on the brane are of two forms: there are quadratic energy-momentum corrections via the tensor $`S_{\mu \nu }`$ and nonlocal effects from the free gravitational field in the bulk $`_{\mu \nu }`$. The nonlocal corrections can be decomposed into a scalar, vector and tensor parts . On FRW brane there is only a scalar part $`𝒰`$ which decays during the expansion of the Universe as $`𝒰a^4`$ where $`a`$ is the scale factor. This means that this term behaves exactly as a radiation fluid, and for this reason it is often called as ”dark radiation” despite its pure geometrical nature. It should be noted, however, that ”dark radiation” can be either positive or negative, and in the latter case it gives rise to some cosmological solutions, impossible in the standard cosmology (for example, the collapse of a flat brane Universe ).

More important differences between standard and brane cosmologies appear at early stage of the Universe evolution, where quadratic corrections proportional to $`S_{\mu \nu }`$ dominate over the standard terms proportional to $`T_{\mu \nu }`$. The rate of expansion of the Universe changes with respect to matter density, this leads to interesting peculiarities of the brane scenario.

The case of matter source in the form of perfect fluid with the equation of state $`p=(\gamma 1)\rho `$ have been extensively investigated during last three years, see, for example, .

The next step is to consider a scalar field on a brane. Even in the standard cosmology the dynamics of Universe filled by a scalar field is much more complicated in comparison to the perfect fluid case . The effective $`\gamma `$ of the field can vary in the whole possible interval $`\gamma (0,2)`$. This property, if the scalar field potential is not too steep, leads to existence of inflation and, in the case of positive spatial curvature, to the possibility of bouncing solutions and dynamical chaos. However, if the potential is steep enough, the situation changes and this transition occurs for potentials as steep as exponential: $`V(\varphi )\mathrm{exp}\lambda \varphi `$ with large enough $`\lambda `$ can not provide inflation  and bounces . It already have been noticed that in brane cosmology we have a similar picture but with inverse square potential $`V(\varphi )=B\varphi ^2`$ playing the role of exponential potentials in the standard cosmology .

In the present paper we consider a brane Universe filled by a scalar field and a perfect fluid. In the ordinary cosmology exponential potentials in this situation have again a particular property – in this case a tracker solution (the energy density of a scalar field tracks exactly the energy density of an ordinary matter) exists . Not surprising that a tracker solution on a brane have been found for the inverse square potential. It was done by K.-I. Maeda  in the particular case of ordinary matter in the form of radiation fluid ($`\gamma =4/3`$). We present a general form of brane tracker solution for an arbitrary $`\gamma `$. We also study other power-law potentials which give us brane scaling solutions – solutions with scalar field energy density scales as some power of a scale factor, when the perfect fluid, being the dominant component, has an energy density scaling as a possibly different power (we use the terminology of ). Such potentials are often used in modern quintessence models  to explain observable acceleration of our Universe.

Our analysis can be useful for understanding evolution of a quintessence field before the radiation-dominated stage, for example, during inflation. Though for the full treatment of this problem it is necessary to consider a more complicated case of a two-field dynamics (quintessence and inflaton fields), we can use the approximation of a perfect fluid if $`\gamma `$ of the matter varies slowly enough.

In the present paper we consider only a flat isotropic brane. Shear and curvature terms can modify the scaling solutions, in the standard cosmology it was shown in . A two-component (a scalar field and a perfect fluid) brane dynamics with nonzero shear and curvature has already become a matter of investigations . We leave the problem of scaling solutions on an anisotropic brane to a future work.

## 2 Tracker solutions

### 2.1 Basic equations

We start with system of equations describing a scalar field on a flat isotropic brane . There are the Raychaudhuri equation

$$\dot{H}=H^2\frac{\kappa ^2}{6}(\rho +3p)\frac{\kappa ^2\rho }{6\lambda }(2\rho +3p)\frac{2𝒰}{\kappa ^2\lambda },$$
(2)
the Friedmann equation

$$H^2=\kappa ^2\frac{\rho }{3}+\kappa ^2\frac{\rho ^2}{6\lambda }+\frac{2𝒰}{\kappa ^2\lambda },$$
(3)
the Klein-Gordon equation for the scalar field

$$\ddot{\varphi }+3H\dot{\varphi }+V^{}(\varphi )=0,$$
(4)
and the evolution equation for the ”dark radiation”

$$\dot{𝒰}+\frac{4}{3}H𝒰=0.$$
(5)
Here $`\lambda `$ is the brane tension. The overall pressure and energy density in (2) – (3) are the sums

$$p=p_m+p_\varphi ,\rho =\rho _m+\rho _\varphi ,$$

where the perfect fluid equation of state $`p_m=(\gamma _m1)\rho _m`$, $`\gamma _m=\mathrm{const}`$ and the scalar field equation of state $`p_\varphi =(\gamma _\varphi 1)\rho _\varphi `$ with $`\gamma _\varphi `$ changing with time. The pressure $`p_\varphi `$ and the energy density $`\rho _\varphi `$ of the scalar field are given by

$$p_\varphi =\frac{\dot{\varphi }^2}{2}V(\varphi ),\rho _\varphi =\frac{\dot{\varphi }^2}{2}+V(\varphi ).$$

Using Eq. (4) we can write the dynamical equation for $`\gamma _\varphi `$ as

$$\dot{\gamma }_\varphi =3H\gamma _\varphi (\gamma _\varphi 2)2\gamma _\varphi \frac{V^{}}{\dot{\varphi }}.$$
(6)
In this section we study the particular potential $`V(\varphi )=B\varphi ^2`$, for this case the Eq.(6) becomes

$$\dot{\gamma }_\varphi =3H\gamma _\varphi (\gamma _\varphi 2)+\frac{\sqrt{2\gamma _\varphi }}{\sqrt{B}}\rho _\varphi (2\gamma _\varphi )^{3/2}.$$
(7)
We would like to notice an interesting difference between scalar field dynamics in ordinary cosmology and on a brane. In standard cosmology a constant coefficient in field potential can be absorbed by rescaling scale factor and time without changing a general dynamical picture. In brane cosmology it is not the case, and coefficient $`B`$ plays an important role, as we will see later.

After introducing the deceleration parameter $`q`$ as

$$\dot{H}=(1+q)H^2,$$

dimensionless variables

$$\mathrm{\Omega }_{\rho m}=\frac{\kappa ^2\rho _m}{3H^2},\mathrm{\Omega }_{\rho \varphi }=\frac{\kappa ^2\rho _\varphi }{3H^2},\mathrm{\Omega }_{\lambda m}=\frac{\kappa ^2\rho _m^2}{6\lambda H^2},\mathrm{\Omega }_{\lambda \varphi }=\frac{\kappa ^2\rho _\varphi ^2}{6\lambda H^2},U=\frac{2𝒰}{\lambda \kappa ^2H^2}$$

and a new time variable $`\tau `$ such that $`\frac{dt}{d\tau }=\frac{1}{H}`$ we can rewrite our system in the form

$$\begin{array}{c}q=\frac{1}{2}\mathrm{\Omega }_{\rho m}(3\gamma _m2)+\frac{1}{2}\mathrm{\Omega }_{\rho \varphi }(3\gamma _\varphi 2)+\mathrm{\Omega }_{\lambda m}(3\gamma _m1)+\hfill \\ \hfill \mathrm{\Omega }_{\lambda \varphi }(3\gamma _\varphi 1)+\sqrt{\mathrm{\Omega }_{\lambda m}\mathrm{\Omega }_{\lambda \varphi }}(3(\gamma _m+\gamma _\varphi )2),\end{array}$$
(8)
$$1=\mathrm{\Omega }_{\rho m}+\mathrm{\Omega }_{\rho \varphi }+U+\mathrm{\Omega }_{\lambda m}+\mathrm{\Omega }_{\lambda \varphi }+2\sqrt{\mathrm{\Omega }_{\lambda m}\mathrm{\Omega }_{\lambda \varphi }},$$
(9)
$$\mathrm{\Omega }_{\rho m}^{}=\mathrm{\Omega }_{\rho m}[2(1+q)3\gamma _m],\mathrm{\Omega }_{\rho \varphi }^{}=\mathrm{\Omega }_{\rho \varphi }[2(1+q)3\gamma _\varphi ],$$
(10)
$$\mathrm{\Omega }_{\lambda m}^{}=\mathrm{\Omega }_{\lambda m}[2(1+q)6\gamma _m],\mathrm{\Omega }_{\lambda \varphi }^{}=\mathrm{\Omega }_{\lambda \varphi }[2(1+q)6\gamma _\varphi ],$$
(11)
$$U^{}=2U(q1),$$
(12)
$$\gamma _\varphi ^{}=3\gamma _\varphi (\gamma _\varphi 2)+\frac{\sqrt{72\gamma _\varphi \mathrm{\Omega }_{\lambda \varphi }}}{\sqrt{B}}m_5^3(2\gamma _\varphi )^{3/2}.$$
(13)
Here the prime denoted derivative with respect to $`\tau `$.

Note, that in the case of several matter sources nonlinear crossing terms appear in the r.h.s. of Eqs. (8) – (9).

In this paper we study only a pure nonstandard regime when the contribution from linear terms proportional to $`T_{\mu \nu }`$ in (1) is negligible in comparison to the quadratic energy-momentum terms proportional to $`S_{\mu \nu }`$. Therefore we put $`\mathrm{\Omega }_{\rho m}`$ and $`\mathrm{\Omega }_{\rho \varphi }`$ to zero. For the potential $`V(\varphi )=B\varphi ^2`$ the Hubble parameter $`H`$ decouples (in standard regime this property belongs to an exponential potential ) and we have a compact phase space unless a ”dark energy” density $`U<0`$.

### 2.2 Equilibrium points

Now we list equilibrium points for the system (8) – (13) which can be stable in future direction with corresponding results of their stability analysis. To simplify the formulae we use the dimensionless coefficient $`bB/m_5^6`$. Remember, that $`\mathrm{\Omega }_{\rho m}`$ and $`\mathrm{\Omega }_{\rho \varphi }`$ are set to zero. As for the full system (8) – (13) the equilibrium points obtained correspond to an intermediate brane behavior, when some asymptotic regime have already been established, but the brane is still driven by the quadratic terms in (1). First, consider the case of $`U=0`$.

$$\gamma _\varphi =\frac{2}{1+b/8},\mathrm{\Omega }_{\lambda m}=0,\mathrm{\Omega }_{\lambda \varphi }=1.$$

Stable for $`\gamma _m>\frac{2}{1+b/8}`$.

$$\gamma _\varphi =\gamma _m,\mathrm{\Omega }_{\lambda \varphi }=\frac{b\gamma _m}{8(2\gamma _m)},\mathrm{\Omega }_{\lambda m}=1+\frac{b\gamma _m}{8(2\gamma _m)}\frac{\sqrt{b\gamma _m}}{\sqrt{2(2\gamma _m)}}.$$

Exists and stable for $`\gamma _m<\frac{2}{1+b/8}`$.

The point I is the regime of scalar field dominance. We see that depending on $`b`$ the scalar field density can mimic an ordinary matter in the whole interval of possible values $`\gamma _\varphi (0,2)`$. Naturally, this point is stable for $`\gamma _m>\gamma _\varphi `$. The point II represents the brane tracker solution with $`\mathrm{\Omega }_{\lambda \varphi }/\mathrm{\Omega }_{\lambda m}=\mathrm{const}`$, which was obtained in  for the particular case of $`\gamma _m=4/3`$.

The case of $`U>0`$ leads to the following changes. The point I disappears for $`b>16`$, the point II becomes unstable for $`\gamma _m>2/3`$. In these regions a new stable point appears:

$$\gamma _\varphi =2/3,\mathrm{\Omega }_{\lambda \varphi }=\frac{b}{16},U=1\frac{b}{16},\mathrm{\Omega }_{\lambda m}=0.$$

Exists for $`b<16`$, stable for $`\gamma _\varphi >2/3`$.

This is an analog of tracker solution when the ”dark energy” plays the role of matter. The ratio $`\mathrm{\Omega }_{\lambda \varphi }/U`$ remains constant. Of course, for $`2/3<\gamma _m<4/3`$, this situation (”dark energy” dominates the ordinary matter) changes after the low-energy stage begins and $`\mathrm{\Omega }_{\rho m}`$ becomes important. Zones of stability for the future direction in the plane $`(b,\gamma _m)`$ for the case of $`U>0`$ are plotted in Fig.1.

If $`U`$ is negative, the phase space becomes incompact. Physically this means that a recollapse of the brane Universe is possible. Our numerical studies indicate that in zones of stability of points I and II the recollapse occurs for large enough initial $`|U|`$, while for small $`|U|`$ the Universe tends asymptotically to corresponding stable points I or II with $`U=0`$. On the contrary, in the region III there are no stable equilibrium points, corresponding to initial negative $`U`$. We have confirmed numerically, that all trajectories from the region $`b<16,\gamma _m>2/3,U<0`$ recollaps. It is worth to remind, however, that all this results are obtained for a pure nonstandard regime. When the full system of equations (8) – (13) is considered, some trajectories may still not recollaps until transition to the standard regime occurs, and, if $`2/3<\gamma _m<4/3`$ the recollapse will never occur.

## 3 General power-law potentials

### 3.1 Scaling solutions

In this section we describe solutions with scalar field energy density scales exactly as some power (in general, different from the fluids one) of the scale factor. The case of scalar field dominance have been described in , therefore we concentrate our effort to the fluid dominance case. Our work generalize  to fluids with $`\gamma 4/3`$.

It is worth to note that these solutions are not future attractor in the sense of the previous section even if we do not consider the transition to the standard regime. Indeed, the dynamical evolution may break down the fluid dominance. Nevertheless, it is an interesting example of transient behavior. This behavior is followed by a regime of a scalar field dominance. Late-time scalar field dominance in the standard regime leads to inflation  and provides a possible explanation of accelerated expansion of our present Universe.

We will follow the method of , where a similar class of solutions have been found in the standard cosmology. Suppose that the scalar field energy density behaves as $`\rho _\varphi a^n`$, where $`a`$ is the scale factor. Since the Klein-Gordon equation is not modified on a brane, this property, as in the standard case, requires that the ratio of a scalar field kinetic energy density and a total scalar field energy density remains constant :

$$\frac{\dot{\varphi }^2/2}{\rho _\varphi }=\frac{n}{6}.$$
(14)
The perfect fluid with equation of state $`p=(\gamma 1)\rho `$ has the energy density $`\rho a^m`$ where $`m=3\gamma `$. In the pure nonstandard brane regime the fluid dominance leads to

$$at^{1/m}.$$
(15)
All our considerations remain valid also for the case of the ”dark energy” dominance (in this case, of course, $`𝒰>0`$), because the ”dark energy” behaves as a ”normal” matter with $`m=4`$. Due to (15), Eq.(4) becomes

$$\ddot{\varphi }=\frac{6}{m}\frac{1}{t}\dot{\varphi }\frac{dV}{d\varphi }.$$
(16)
Eq. (14) gives us

$$\dot{\varphi }t^{n/2m}.$$

Integrating this equation we obtain

$$\varphi =At^{1n/2m}.$$
(17)
Substituting (17) into (16) and solving for $`V(\varphi )`$ we find the potentials which allow the scaling behavior in the nonstandard brane regime

$$V(\varphi )=A^{2\alpha }\frac{2(3(\alpha 2)\alpha m)}{\alpha m(\alpha 2)^2}\varphi ^\alpha ,$$
(18)
where

$$\alpha =\frac{2n}{n2m}.$$
(19)
It means that, as in the standard cosmology, scaling solutions exist for power-law potentials $`V(\varphi )\varphi ^\alpha `$. If $`\alpha `$ is given, the scalar field energy density is proportional to $`a^n`$, where

$$n=\frac{2\alpha }{\alpha 2}m.$$
(20)
Though this expression differs only for the factor $`2`$ in the numerator from its analog in the standard cosmology , the properties of brane scaling solutions are rather different. Let us begin with the case of negative $`\alpha `$. In the standard scenario the energy density of a scalar field with $`\alpha <0`$ always falls slower than the density of a perfect fluid while in the brane cosmology it falls slower for $`\alpha >2`$ and faster for $`\alpha <2`$. Second, the positivity of the potential in (18) requires that $`\alpha (3m)<6`$, which indicates that for $`m>3`$ (a matter with a positive pressure) and

$$\alpha <\frac{6}{3m}$$
(21)
the scaling solution does not exist. On the contrary, in the standard cosmology the scaling solution exists for an arbitrary negative $`\alpha `$ independently on $`m`$.

Consider now the case of positive $`\alpha `$. The positivity of (18) leads to $`\alpha (3m)>6`$. It means that for $`m>3`$ scaling solutions do not exist for any $`\alpha `$. In particular, they do not exist for the case of $`\gamma =4/3`$, investigated in . On the other side, for $`m<3`$ the situation is similar to the standard cosmology – scaling solution exists for sufficiently steep potential with $`\alpha `$ satisfying in the brane case

$$\alpha >\frac{6}{3m}.$$
(22)
### 3.2 Stability analysis

We have to analyze the stability of the scaling solutions (17) with a power-law potential (18). To do this we use new variables

$$\tau =e^t,u(\tau )=\frac{\phi (\tau )}{\phi _0(\tau )},$$
(23)
where $`\phi _0(\tau )`$ is the exact solution given by (17) Then the equation for scalar field becomes

$$u^{\prime \prime }+u^{}\left(\frac{2+\alpha }{2\alpha }+\frac{3}{m}\right)+\frac{2}{\alpha 2}\left(\frac{\alpha }{\alpha 2}\frac{3}{m}\right)\left(uu^{\alpha 1}\right),$$
(24)
where the prime denotes the derivative with respect to $`\tau `$. This second-order differential equation can be divided into the autonomous system

$$\begin{array}{c}u^{}=p,\hfill \\ p^{}=\frac{2}{\alpha 2}\left(\frac{3}{m}\frac{\alpha }{\alpha 2}\right)\left(uu^{\alpha 1}\right)\left(\frac{2+\alpha }{2\alpha }+\frac{3}{m}\right)p.\hfill \end{array}$$
(25)
Scaling solution corresponds to a critical point $`(u,p)=(1,0)`$. Linearizing (25) about this point we find the eigenvalues of these coupled equations

$$\lambda _{1,2}=\frac{\alpha +2}{2(\alpha 2)}\frac{3}{2m}\pm \sqrt{\left(\frac{\alpha +2}{2(\alpha 2)}\frac{3}{2m}\right)^2+\frac{2\alpha }{\alpha 2}\frac{6}{m}}.$$
(26)
Then the condition for stability is the negativity of the real part of both eigenvalues.

For $`\alpha <0`$ after analysis we find that scaling solution is a stable attractor for all region of existence (21). And for $`\alpha >0`$ there is an additional inequality

$$\alpha >\frac{2m+6}{3m},$$
(27)
which excludes a narrow domain from (22). For example, for $`V\varphi ^4`$ the scaling solution is stable if $`m<1`$, or, equivalently, $`\gamma _m<1/3`$, this is exactly the condition for a brane inflation driven by an ordinary matter. For $`V\varphi ^6`$ we have $`m<3/2`$ (or $`\gamma _m<1/2`$).

Zones of stability of the scaling solutions on the plane $`(\alpha ,m)`$ are plotted in Fig.2.

### 3.3 Other solutions

In the end of this section we briefly describe solutions, existing in the regions of the plane $`(\alpha ,m)`$ where scaling solutions do not exist.

For positive $`\alpha `$ with even $`n`$ we have infinite damping oscillations of the scalar field, similar to the standard case. Some asymptotic formulae for potentials $`V=\frac{1}{2}m^2\varphi ^2`$ and $`\frac{1}{4}\lambda \varphi ^4`$ have been given in .

For negative $`\alpha `$ the other type of solution exists. In this kinetic-term dominated solution the contribution of the potential for the energy density of the scalar field is negligible in comparison to the contribution of the kinetic term. Its explicit form is

$$\varphi =\varphi _0\left(\frac{t}{t_0}\right)^{\frac{m3}{m}},$$
(28)
where $`\varphi _0`$ and $`t_0`$ are integration constants. This solution does not exist for $`m3`$. In this case the scaling solution is stable for all negative $`\alpha `$. For $`m=4`$ we recover the known behavior $`\varphi t^{1/4}`$ .

## 4 Conclusions

We have investigated a scalar field dynamics with a power-law potential on a brane filled by a perfect fluid with equation of state $`p_m=(\gamma _m1)\rho _m`$ in a regime when quadratic energy-momentum correction dominates. The tracker solution when the ratio of scalar field and fluid energy densities remains constant is found for the scalar field potential $`V(\varphi )=B\varphi ^2`$. Two other future asymptotic for this case are described and zones of their stability in the plane $`(B,\gamma _m)`$ are found.

We also describe scalar field evolution with $`V(\varphi )\varphi ^\alpha `$ in a brane Universe, dominated by perfect fluid energy density. We have found scaling solution, analogous to known in the standard cosmology. The main differences between standard and brane cosmology are:

* For negative $`\alpha `$ in the standard cosmology the scaling solutions exist for arbitrary $`\gamma _m`$ and $`\alpha `$, whereas in the brane cosmology they do not exist for $`\gamma _m>1`$ and large enough $`|\alpha |`$, which do not satisfy (21).
* For positive $`\alpha `$ in the standard cosmology the scaling solutions exist for any $`\gamma _m`$ if $`\alpha `$ exceeds some value, depending on $`\gamma _m`$. In the brane cosmology the scaling solutions exist only for $`\gamma _m<1`$ and large enough $`\alpha `$, satisfying (27). They are absent for $`\gamma _m>1`$.

It is also interesting, that in the case of negative $`\alpha `$ the scalar field energy density may fall more rapidly or more slowly than the fluid density, depending on the sign of $`\alpha +2`$, whereas in the standard scaling solution with negative $`\alpha `$ the ratio of scalar field and fluid energy densities always increases. This was discovered in the particular case of radiation fluid in , where it was argued that this fact may be important for quintessence scenario. In our paper we confirm that this property do not depend on the state equation of the perfect fluid.

## Acknowledgments

The work was supported by RFBR grants Ns. 00-15-96699 and 02-02-16817.


format: text

text:
# Atomic ordering effects and bond indices in quaternary systems Si-(Be,Al)-O-N

## I Introduction

Atomic ordering effects (AOE) have been discovered in a wide range of non-stoichiometric and multi-component alloys, solid solutions (SS) and compounds. The influence they produce on the properties is great and comparable with the effect of altering the chemical composition of materials gus .

In the majority of theoretical simulations of multi-component materials, atomic ordering effects are studied using integral energy characteristics (for example, free, total or cohesion energies) of a limited number of atomic configurations of a system set a priory order . Generally, this approach is applied to binary or ternary compounds, whereas for more complex systems its application is confined by computer power due to rapidly increasing sizes of cells to be considered and the number of atomic configurations.

In the present report, we describe a new method for studying AOE based on comparative analysis of detached types of interatomic bonds in multi-component systems. In the framework of this method, optimal atomic configurations are chosen by minimizing unfavorable interatomic bonds.

The method is illustrated using quaternary solid solutions Si6-xAlxOxN8-x ($`\beta `$-SiAlONs) considered earlier in ref. okatov2001  as an example. It is also applied to predict atomic ordering effects in the hypothetical $`\beta `$-Si3N4-based solid solutions Si6-xBexO2xN8-2x.

## II The method for modeling atomic ordering effects: $`\beta `$-SiAlONs

As is known, the stability of a thermodynamic system corresponds to the free energy or, at zero temperature, total energy minimum condition. In terms of the chemical bonding theory, this condition implies the optimal type of bands filling, when all bonding states are occupied and all antibonding states are vacant.

In the proposed method, AOE are determined by searching for the multi-component system configurations having the maximum number of the most favorable interatomic bonds. To characterize the detached type of paired bonds, crystal orbital overlap populations (COOP) and the numbers of filled antibonding states (FAS) may be used. In this work, these values are calculated employing the Mulliken analysis by the tight-binding band method with matrix elements parameterization according to the Hückel theory huck . Other techniques of bond characterization can be used as well.

The application of this approach can be demonstrated on the four-component SS in the Si-Al-O-N system with the basic structure $`\beta `$-Si3N4 (the so-called $`\beta `$-SiAlONs of the formal stoichiometry Si6-xAlxOxN8-x, where $`x=04.2`$).

A previous investigation of $`\beta `$-SiAlONs by the tight-binding band method okatov2001  showed that Al and O atoms form extended quasi-one-dimensional (1D) structures constituted by 12-atomic rings — the so-called “impurity channels”, fig. 1. Subsequent studies by the DFT cluster ryz  and the band structure full-potential LMTO methods flmto  showed that ordered structures exhibited the best stability. Besides, their chemical bonding and electronic properties were examined.

A detailed analysis of bond indices in $`\beta `$-SiAlONs demonstrated okatov2001  that the most unfavorable bonds are Si-O, which contain the maximum number of FAS. The antibonding states of Al-O bonds are less populated. Si-N and Al-N bonds do not contain any FAS, while COOP(Si-N$`=0.58`$ e)$`>`$COOP(Al-N$`=0.24`$ e). Hence, a bonds hierarchy taking into account their stabilizing effects can be represented by the so-called sequence of bonds preferences:

$$P_{\mathrm{Si}\mathrm{O}}<P_{\mathrm{Al}\mathrm{O}}<P_{\mathrm{Al}\mathrm{N}}<P_{\mathrm{Si}\mathrm{N}},$$
(1)
where $`P_{\mathrm{A}\mathrm{B}}`$ reflects the preference of A-B-type bonds presence in the system.

Thus, in order to determine the superstructure of a multi-component material, it is necessary to choose such distributions of atoms over lattice sites, in which the number of bonds is minimized consecutively in accordance with (1).

For example, for $`\beta `$-SiAlONs at the first step the superstructures with the minimum number of Si-O bonds are chosen from all possible ones according to (1). Among them we choose the superstructures with the minimum number of Al-O bonds etc. As a result, in a system of a given composition, the superstructures are found, which possess the same number of the same types of bonds, but differ only by their positions in the supercell.

To determine the atomic order type, the following algorithm is applied. The SS Si6-xAlxOxN8-x of a certain composition modelled by a supercell is considered. In this work we make use of the 126-atomic supercells (3x3x1) analogous to those applied in okatov2001 . For example, the composition of the SS with $`x=2.0`$ (Si4Al2O2N6) is described by the Si36Al18O18N54 cell. The following items are fulfilled:

Atoms (Si,Al,N,O) are randomly distributed over the supercell sites so that the result does not depend on initial conditions.

A random substitution, which is the only permitted one for this system (for example Si$``$Al, N$``$O), is performed, while all other types of substitution (Si$``$O, N$``$Al) are not considered.

If the correspondence to (1) improves (the number of favorable bonds increases and the number of unfavorable bonds3 decreases), the step is accepted, otherwise it is rejected.

Items (ii) and (iii) are repeated until the number of steps denied in succession does not exceed $`N_1`$.

Items (i-iv) are repeated until the number of structures denied in succession does not exceed $`N_2`$. This cycle is necessary since the result of cycle (iv) depends on the initial position of atoms and on the first accepted substitution (item iii) defining the subsequent shape of the superstructure.

The choice of $`N_1`$, $`N_2`$ values is determined by the supercell size and composition. In the present work they range from 10000 to 100000 and from 10 to 400 respectively.

An application of the above algorithm for $`\beta `$-SiAlONs with the same cells and compositions as in okatov2001  shows that the resulting superstructures shapes (1D impurity channels, fig. 1) are in complete agreement with those obtained in a series of complicated calculations of total energies. The formation of impurity channels in $`\beta `$-SiAlONs leads to the disappearance of the most unfavorable Si-O bonds, while the number of Al-O bonds is minimal.

The evident drawback of the aforesaid scheme concerns condition (1), which reflects only the effect of paired bonds A-B. It should not be excluded that in some cases the formation of AnBm clusters is more preferable in spite of low preference of A-B bonds. This restriction may be removed by using a sequence of clusters preferences instead of a sequence of bonds preferences. In respect to $`\beta `$-SiAlONs, the above restriction does not permit differentiating between distant and adjacent positions of impurity channels.

## III Electronic structure and atomic ordering effects in the $`\beta `$-Si3N4-Be-O system

### III.1 Bond indices and chemical composition

At the first stage of modeling the hypothetical SS in the quaternary system $`\beta `$-Si3N4-Be-O, we analyzed the electronic properties and bond indices of pure $`\beta `$-Si3N4 and that containing Be and O impurities, fig. 2. It is seen that in the binary nitride $`\beta `$-Si3N4 (cell Si54N72, valence electrons concentration (VEC) equals 576 e/cell), all antibonding states are vacant, while all bonding states are completely occupied.

O$``$N or Be$``$Si substitutions (cells Si54ON71 and Si53BeN72 respectively) reduce the stability of the $`\beta `$-Si3N4:O and $`\beta `$-Si3N4:Be systems. In the first case this is due to partial filling of antibonding states when VEC increases (577 e/cell), while in the second case — due to a depopulation some bonding states when VEC decreases (574 e/cell). A “mutual compensation” of the destabilizing effects of single impurities Be, O can be achieved by choosing a composition, at which the VEC remains 576 e/cell, i.e. when the impurities are introduced in the proportion (Be+2O). In this case the type of bands filling (as compared with the basic $`\beta `$-Si3N4) does not change and the number of filled antibonding and empty bonding states is minimal, fig. 2. Thus, the formal stoichiometry of the SS in the $`\beta `$-Si3N4-Be-O system is Si6-xBexO2xN8-2x.

### III.2 Short-range atomic ordering

When heterovalent substitutions occur in $`\beta `$-Si3N4, the adjacent arrangement of donour (O) and acceptor (Be) atoms (i.e. the formation of {BeO2} complexes) is more probable, otherwise the charge screening effect hinders the above mentioned compensation of bonds filling.

### III.3 Long-range atomic ordering

Analysis of separate bonds in the Si6-xBexO2xN8-2x system shows that Si-O and Be-O bonds consist of filled antibonding states, and their number in Si-O is greater than in Be-O (fig. 2). Determining atomic ordering effects by the proposed algorithm, we search for atomic configurations satisfying the condition of type (1):

$$P_{\mathrm{Si}\mathrm{O}}<P_{\mathrm{Be}\mathrm{O}}<P_{\mathrm{Be}\mathrm{N}}<P_{\mathrm{Si}\mathrm{N}},$$
(2)
In other words, for the considered series of solid solutions of variable concentrations Si6-xBexO2xN8-2x ($`x00.5`$, supercells Si53BeO2N70 – Si24Be30O60N12), we determine the configurations of modelled cells, in which the number of less advantageous bonds Si-O and Be-O is successively minimized.

The results of the simulation are depicted in fig. 3 and table 1. Their major difference from $`\beta `$-Si6-xAlxOxN8-x is the presence of Si-O bonds, the number of which grows monotonously with the concentration of Be and O. This indicates that Si6-xBexO2xN8-2x is less stable in comparison with Si6-xAlxOxN8-x. When $`x`$ is small (cells Si53BeO2N70 – Si51Be3O6N66), the system contains Be-N bonds, and the number of “unfavorable” Si-O bonds grows with $`x`$ most rapidly. The noted peculiarities are due to a small number of Be and O atoms in the supercell and its limited size. In real crystals this situation is not possible.

The smallest cluster, which does not contain any Be-N bonds, is made up of four {BeO2} complexes (supercell Si50Be4O8N64, fig. 3). As the number of Be, O atoms grows further, the number of Si-O bonds increases linearly and equales the double number of {BeO2} complexes in the cell. This fact shows that the system has no definite order because, for example, one cluster of 8 complexes {BeO2} contains the same set of bonds as two separate clusters each containing 4 {BeO2} complexes (fig. 3). As a result, atomic configurations at a definite supercell composition may have the form of extended, compact or multiple clusters constituted by {BeO2} complexes, fig. 3.

The analysis of bonds allows us to conclude that the formation of a whole compact cluster with the maximal size is most probable. In this case, the “depopulation” of FAS occurs analogously as in precipitated phases of concentration polytypes in the AlN-Si-O system refr ; metal , where antibonding states are present only in “interlayer” bonds Al-O and are absent within the “block of defects”.

In Si6-xBexO2xN8-2x, an increase in the size of the “impurity cluster” brings about a growth of the number of “intracluster” Be-O bonds reducing the number of FAS. For example, COOP spectra of “intracluster” and “frontier” Be-O bonds of the Si46Be8O16N56 supercell are depicted in fig. 4. It can be seen that intracluster bonds do not possess any antibonding states. This indicates that the formation of a compact cluster of the largest size is more advantageous. The growth of a cluster in Si6-xBexO2xN8-2x is not limited meaning the formation of a heterophase state of the system ($`\beta `$-Si3N4/BeO).

Thus, the proposed semiempirical approach based on quantum-chemical analysis of separate interatomic bonds makes it possible to determine the most probable effects of short-range (clusterization) and long-range (superstructures) atomic ordering in multi-component solid solutions and to find out the factors responsible for their stability.

Analysis of the differences between Si6-xBexO2xN8-2x and Si6-xAlxOxN8-x solid solutions shows that the impossibility of formation of homogeneous SS in the $`\beta `$-Si3N4-Be-O system is due to the excess of oxygen atoms with regard to those in $`\beta `$-SiAlONs. This leads to the presence of Si-O bonds in the system, which does not allow the impurities to form compact structures of a limited size like impurity channels. Therefore it can be suggested that SS in the $`\beta `$-Si3N4-Be-O system may be stabilized by altering the stoichiometry of the compounds (composition Si6-xBexOxN8-x) or by partial replacement of Si atoms by elements of the IV-th group, for example, Zr (Si6-xBex/2Zrx/2OxN8-x).

Undoubtedly, the detailed analysis of the potentially stable structures, which were suggested in the framework of the used approach, implies further investigations of their energy states by ab-initio methods for the conformation of the presented results. At present, the above studies of Si6-xBexO2xN8-2x, Si6-xBexOxN8-x and Si6-xBex/2Zrx/2OxN8-x SS are being performed by the full-potential LMTO method.

## IV Conclusions

In this paper we proposed a new microscopic semi-empirical method for describing atomic ordering effects in multi-component systems. This method is based on the analysis of interatomic bonds of separate types and includes a procedure of minimizing the number of unfavorable bonds.

The method is used to study AOE in solid solutions of variable compositions: Si6-xAlxOxN8-x ($`\beta `$-SiAlONs) and hypothetical Si6-xBexO2xN8-2x in the $`\beta `$-Si3N4-Be-O system. It is established that in $`\beta `$-SiAlONs Al, O atoms form separate quasi-one-dimensional “impurity channels”. On the contrary, in Si6-xBexO2xN8-2x Be, O atoms trend to form a whole compact cluster. This implies that in the $`\beta `$-Si3N4-Be-O system homogeneous solid solutions cannot be formed under equilibrium conditions.

Potential ways of stabilizing the SS in the $`\beta `$-Si3N4-Be-O system are supposed to be a variation of its stoichiometry (reduction in the oxygen concentration) or partial replacement of silicon by other atoms from the IV-th group, for example, zirconium.

###### Acknowledgements.

This work was supported by the Russian Foundation for Basic Research, grant # 01-03-96515 (Ural).


text:
# A semiclassically entangled puzzle

## 1. introduction

When Einstein, Podolsky and Rosen  explored the “strange” behaviour of entangled states, their main goal was to set forth an argument concerning the ontology of quantum mechanics. Only quite later, did Bell  show that some of the discussion had a definite empirical bearing, namely, on matters related to the nonlocal character of quantum mechanics. After the empirical confirmation of this latter , emphasis has been placed on the “peaceful coexistence”  of this kind of nonlocality and the partial order of spacetime events derived from the principle of relativity .

In terms of the Wigner  function $`W`$, the cartesian Weyl  representation of the density operator $`\widehat{\rho }`$ as a real function on euclidean phase space, this question is addressed in the following way: let $`x_1`$ and $`x_2`$ denote the phase space points of two subsystems whose total system is in an entangled eigenstate of the hamiltonian $`\widehat{H}_1+\widehat{H}_2`$ , with corresponding Wigner function $`W_0(x_1,x_2)`$. Also, let $`\widehat{H}_2^{}`$ be a new hamiltonian, whose Weyl symbol $`H_2^{}`$ generates a flow $`x_2(0)x_2(t)`$, and let $`W_t`$ denote the evolution of $`W_0`$ effected by this new hamiltonian. Further: let $`W_0^1(x_1)W_0(x_1,x_2)𝑑x_2`$ and $`W_t^1(x_1)W_t(x_1,x_2)𝑑x_2`$ be the corresponding reduced Wigner functions, which encode all the empirical information of subsystem $`1`$ alone, and denote their difference by $`\delta W^1(x_1,t)`$.

In this context, conservation of the relativistic partial order of events is implied by

(1)
$$\delta W^10,$$

which is always true in full quantum mechanics $`(\mathrm{}=1)`$, due to its linear structure . In a classical setting $`(\mathrm{}=0)`$, on the other hand, a sufficient condition for the validity of (1) is given by

(2)
$$W_0(x_1,x_2(0))W_t(x_1,x_2(t))=W_0(x_1,x_2(0)).$$

Years ago, Heller  noticed that the generic expression corresponding to (2) is not always valid at the semiclassical level $`(\mathrm{}0^+)`$, a fact which has recently been reinterpreted in a more geometrical way . Accordingly, for a generic semiclassical Wigner function $`𝒲`$, (2) is generally true only if $`H_2^{}`$ is quadractic. For a nonlinear interaction, (2) is valid only for points $`x_2`$ close to the Lagrangian leaf corresponding to the classical limit of $`W`$, whenever its semiclassical expression $`𝒲`$ is an oscillatory function of $`x_2`$, as is usually the case for pure states .

The purpose of this letter is to state the following puzzle: at the semiclassical level, it is not obvious whether (1) is valid in general. To see this puzzle, we focus on a particularly simple semiclassically entangled state, namely the (anti)symmetric state of two identical $`1`$-dimensional harmonic oscilators, under a small nonlinear interaction $`H_2^{}`$. We apply the results of  to ask whether, in a suitable limit, (1) might not hold. Finally, we emphatically stress the highly speculative nature of this letter.

## 2. the setup

We study the semiclassical dynamics of an eigenstate of $`\widehat{H}=\widehat{H}_1+\widehat{H}_2`$ , where $`\widehat{H}_i=(\widehat{p}_i^2+q_i^2)/2`$, and of $`\widehat{P}_{12}:12`$. This state is: $`\sqrt{2}|\mathrm{\Psi }=|e_m,1|e_n,2\pm |e_n,1|e_m,2`$ , where the $`\pm 1`$ refers to the eigenvalue of $`\widehat{P}_{12}`$ and where $`e_{m+n}=e_m+e_n`$ is the eigenvalue of $`\widehat{H}`$ , with $`e_m=(m+1/2)\mathrm{}`$ , $`e_n=(n+1/2)\mathrm{}`$ and where the energy difference is $`\lambda =e_me_n=l\mathrm{}`$, with $`mn=l>0`$.

The corresponding Wigner function factors as $`2W_{m+n}(x_1,x_2)=W_m(x_1)W_n(x_2)+W_n(x_1)W_m(x_2)\pm M_m^n(x_1)M_n^m(x_2)\pm M_n^m(x_1)M_m^n(x_2)`$ , where $`M_m^n`$ is the Moyal  function, or cartesian Weyl representation of the transition operator $`|e_me_n|`$, satisfying $`M_m^n=\overline{M_n^m}`$. We write this function in polar $`\{x=(r,\theta )\}`$ coordinares but emphasize that we are not considering the action-angle Weyl representation, which is not equivalent . In taking the semiclassical limit 111There are different phenemenological ways of approaching the classical limit. The one considered here, large quantum numbers or small effective Planck’s constant, should not be identified with the limit of many degrees of freedom, or large complexity, for which most present theories of decoherence apply., we formally let $`\mathrm{}0^+`$ . This is equivalent to letting $`m,n\mathrm{}`$ , while keeping $`e_m,e_n<\mathrm{}`$ , or to varying an “effective” $`\mathrm{}`$.

Our analysis is based on the semiclassical expression for the Moyal function which is not corrected by uniform approximations on the caustics. Ours is a particularly simple case of the general expression carefully treated in  and is given by (we use curly letters to indicate semiclassical expressions):

(3)
$$_m^n(x)=\frac{e^{il\theta }cos(\varphi _{m,n}(r)/\mathrm{}\pi /4)}{\sqrt{\pi ^3\mathrm{}/2}D_{m,n}(r)},$$

where $`\varphi _{m,n}(r)`$ is half of the symplectic area between the circle with radius $`\sqrt{(2m+1)\mathrm{}}`$ centered at the origin and the circle with radius $`\sqrt{(2n+1)\mathrm{}}`$ centered at $`(2r,\theta )`$ , and $`D_{m,n}(r)=|\dot{x}_m^+\dot{x}_n^{}|^{1/2}`$, where $`\dot{x}_m^+`$ , $`\dot{x}_n^{}`$ are the phase space velocity vectors of the intersecting circles. These intersections are real as long as $`r_{}<r<r_+`$, where $`r_\pm =(\sqrt{(2m+1)\mathrm{}}\pm \sqrt{(2n+1)\mathrm{}})/2`$ , and the denominator goes to zero when the circles become tangent, i.e. $`r=r_\pm `$ , which are the caustic lines. Outside this ring, $`_m^n0`$ , in this crude expression. Naturally, the corresponding formula for $`𝒲_m`$ is obtained from (3) with $`m=n`$. By performing uniform approximations as the intersection points coalesce, $`_m^n`$ can be expressed in terms of Airy functions which do not blow up and decay exponentially outside the ring. Still using (3), the sum of the cross terms in $`𝒲_{m+n}`$ is expressed as

(4)
$$Re\{_m^n(x_1)_n^m(x_2)\}=\frac{cos(l(\theta _1\theta _2))cos(\varphi _{m,n}(r_1)/\mathrm{}\pi /4)cos(\varphi _{m,n}(r_2)/\mathrm{}\pi /4)}{(\pi ^3\mathrm{}/2)D_{m,n}(r_1)D_{m,n}(r_2)}.$$

Of course, $`Re\{_m^n(x_1)_n^m(x_2)\}𝑑x_2=0`$ and so $`𝒲_{m+n}^1(𝒲_m+𝒲_n)/2`$ , as it should be.

## 3. the puzzle

The seeming puzzle to be laid forth below follows from the fact that $`_m^n`$ is an oscillatory function, that is, (3) has an essential singularity at $`\mathrm{}=0`$. This is the reason why (2) fails to hold when $`\mathrm{}0^+`$ as can be seen by directly plugging (3) into the Moyal bracket. Its successive derivatives then bring negative powers of $`\mathrm{}`$ which alter the usual $`\mathrm{}`$-expansion obtained for smooth functions. Recently , a simple geometrical prescription for a more correct semiclassical propagation of such functions was obtained, in agreement with another old result derived by Berry and Balazs  according to which one should propagate the whole family of curves, or lagrangian submanifolds, associated with a quantum state and then re-evaluate the semiclassical expressions accordingly.

The new result obtained in  somehow simplifies this procedure by showing that, equivalently, one can obtain a more correct semiclassical propagation of an oscillatory function like $`_m^n(x)`$ by classically propagating the tips of the chord centered on $`x`$ instead of classically propagating the argument $`x`$ itself, as in the Liouville case (2). Of course, Lioville propagation is still semiclassically correct for smooth functions, when all chords are null, or if the flow is linear, as discussed in more detail in . Since $`_m^n(x)`$ is not smooth, we might expect surprises when under nonlinear interactions. And these might be related, in some way or another, to the old strange fact that for generic pure-state Wigner functions under nonquadractic hamiltonians, if the interaction is first performed in the full quantum regime and then the classical limit is taken, the final result generally appears to be rather different from first taking the classical limit and then interacting.

Therefore, let us now suppose that subsystem $`2`$ interacts nonlinearly for an interval of time $`t`$. For simplicity, we state the puzzle when the interaction is either cubic $`H_2^{}=\alpha q_2^3/3`$ or quartic $`H_2^{}=\alpha q_2^4/4`$ and denote by $`ϵ=\alpha t`$ the strength of this interaction. After rewriting (4) as a sum of oscillatory functions of $`\varphi _{m,n}^\pm (r_2,\theta _2)=\varphi _{m,n}(r_2)\pm \lambda \theta _2`$ , it is not too difficult to see that, within the stationary phase regime $`(\mathrm{}0^+)`$, the effect of a small interaction is, in a first approximation, the equivalent to having only phase shifts $`\varphi _{m,n}^\pm (r_2,\theta _2)\varphi _{m,n}^\pm (r_2,\theta _2)+\delta \varphi _{m,n}^\pm (r_2,\theta _2)`$.

More precisely, to get these phase shifts we first apply equation (16) in  to the chord $`\xi _{m,n}^\pm =J(\varphi _{m,n}^\pm (x_2)/x_2)`$, centered on $`x_2`$, and get the “midpoint phase difference” $`\stackrel{~}{\varphi }_{m,n}^\pm (\stackrel{~}{x}_2^t)\varphi _{m,n}^\pm (x_2)`$, where $`\stackrel{~}{x}_2^t`$ is the midpoint of the hamiltonian flow of the tips of the chord $`\xi `$. Then we note that, for small $`ϵ`$, the total phase difference acquired for the hamiltonian flow of $`x_2`$ is approximately half of $`\stackrel{~}{\varphi }_{m,n}^\pm (\stackrel{~}{x}_2^t)\varphi _{m,n}^\pm (x_2)`$. The amplitude being approximately covariant under the hamiltonian flow, for small $`ϵ`$, this determines the looked-for phase shifts as

$$\delta \varphi _{m,n}^\pm (r_2,\theta _2)(2ϵ/3)\{sin(\theta _2)[e_{m+n}(\lambda /2r_2)^2r_2^2]^{1/2}\pm cos(\theta _2)[\lambda /2r_2]\}^3,$$

for the cubic case, while for the quartic case we get

$$\delta \varphi _{m,n}^\pm (r_2,\theta _2)(2ϵ)r_2cos(\theta _2)\{sin(\theta _2)[e_{m+n}(\lambda /2r_2)^2r_2^2]^{1/2}\pm cos(\theta _2)[\lambda /2r_2]\}^3.$$

The next step is to compute $`Re\{_m^n(x_1)_n^m(x_2)\}𝑑x_2`$ and check if it vanishes, as before. But note that the phase shifts have broken the symmetry of $`Re\{_m^n(x_1)_n^m(x_2)\}`$ with respect to $`\theta _2`$ and it is not obvious now whether the integral vanishes or not. In fact, we must now perform the new integration $`Re\{_m^n(x_1)_n^m(x_2)\}𝑑x_2`$ by stationary phase and it is not difficult to see that the critical points of $`\varphi _{m,n}^\pm (r_2,\theta _2)+\delta \varphi _{m,n}^\pm (r_2,\theta _2)`$ are generally complex. This means that we must generally perform a 2-dimensional steepest descent evaluation of the integral, carefully deforming the whole ring into the complex 2-plane, counting all the significant contributions, adding them all up… This is a rather delicate computation, but it is not unreasonable to suppose that, instead of zero, we should at least get an exponentially small result , as in a “tunneling” effect 222 For a much simpler analogy, consider the integral $``$ $`=`$ $`_0^{2\pi }cos(l\theta )𝑑\theta `$ . Writing $`l\theta =\varphi (\theta )/\mathrm{}`$ , where $`\varphi (\theta )=\lambda \theta `$, then under the phase shift $`\varphi \stackrel{~}{\varphi }=\varphi +\delta \varphi `$, where $`\delta \varphi =ϵ\theta ^3/3`$ , we have that $`\stackrel{~}{}`$ is asymptotically $`(\mathrm{}0^+)`$ mapped into the Airy integral, that is, $`\stackrel{~}{}`$ $``$ $`\pi (\mathrm{}/ϵ)^{1/3}𝒜_i(y)`$ , where $`y=\lambda ϵ^{1/3}\mathrm{}^{2/3}`$..

This being an asymptotic analysis, however, we still need to make sure whether the result of this stationary phase computation implies that $`\delta 𝒲^1`$ is really different from zero, to lowest order of $`\mathrm{}`$ (remember that (3) is only the first term of an $`\mathrm{}`$-expansion), in some region of the positive $`(\lambda ,ϵ,\mathrm{})`$-space close to the $`\mathrm{}=0`$ plane. Also in this region, we must analyse the behaviour of $`𝒲_m(x_2)𝑑x_2`$ and $`𝒲_n(x_2)𝑑x_2`$ under similar phase shifts to find out how unitary such a propagation of $`𝒲^1`$ actually is (note that $`Re\{_m^n(x_1)_n^m(x_2)\}𝑑x_1𝑑x_2`$ vanishes). Then we must compare all the above results and further analyse and interpret them (for instance, in possible relation to Stokes phenomenon).

Finally, to properly account for stronger interactions ($`ϵ`$ not so small), closer resonance ($`\lambda `$ small), or to actually get more precise estimates for $`\delta 𝒲^1`$ (any $`ϵ`$ or $`\lambda `$), we must leave the crude expression (3) and proceed much more carefully with the uniform approximations (this is the really difficult part), since non-Liouville propagation of the amplitude is so much more relevant near caustics.

It is very important to emphasize that any possible nontrivial answer for $`\delta 𝒲^1`$ would only be an asymptotic one. Such a possibility seems due to the $`\theta _2`$-symmetry breaking of the whole semiclassical Wigner function $`𝒲_{m+n}(x_1,x_2)`$, effected by the nonlinear interaction with $`H_2^{}`$ , and the consequences of this symmetry breaking should become more explicit as the amplitude corrections are included. But, actually, this is precisely the very hard crucial point that needs to be fully investigated. In other words, it is quite likely that only a careful (and very difficult) analysis based on more accurate semiclassical approximations can produce clear answers, even because the accuracy of the chord expression (3), from which we started, is in fact worse than exponential.

Then, if such a possibility is indeed correct, we could conjecture there would be an optimal way of approaching the classical limit, in the positive $`(\lambda ,ϵ,\mathrm{})`$-space, leading to a maximal $`\delta 𝒲^1`$. Actually, combined with the singular nature of the semiclassical limit, the presence of so many independent parameters is what seems to open this problem in such a manner. And for more degrees of freedom, issues from the classical (non)integrability of a subsystem, with respect to the nonlinear interaction, might well play a further significant role. However, it is clear that nothing in the previous discussion proves or even suggests that a nontrivial $`\delta 𝒲^1`$ can probably happen. In this letter, we have merely pointed out to an apparent possibility that this somehow might come to pass.

## 4. speculations

Final words on this puzzle from the phenomenological, empirical and philosophical points of view.

Phenomenologically, how could it ever be possible for (1) to fail at all, since it is valid in both the classical and the full quantum regimes 333Examples of phenomena that appear in one or the other but not in both regimes are easier to come by. On the other hand, the non-geometrical phenomenon of the glory  can be seen as a phenomenon that vanishes both “classically” (ray regime) and “fully quantum” (long wave regime) but not “semiclassically” (short wave asymptotics). ? Well, we could argue that (1) is valid in both regimes for different reasons. In the classical case, (2) translates the fact that no quantum correlations have “survived”. But this can be seen as an effect of the correlations becoming “random”, which is best expressed by the oscillations in the Wigner function becoming so high that they average to zero everywhere but on the classical leaf. In the full quantum case, by contrast, the correlations are fully structured, in other words, they have a definite “form”, which is usually expressed by the concept of a definite normalized vector in the Hilbert space. Such a vector can only be preserved by linear transformations, from which (1) follows. Therefore, the puzzle ammounts to a question of whether, at the semiclassical level, the linear structure of quantum mechanics could feel some disturbing effects of a nonlinear classical dynamics and, if so, how they could show up 444In a certain way, the two important factors for a Bell’s telephone, namely the distinction between pure and mixed states and between linear and nonlinear evolution, show up in the semiclassical propagation of Wigner functions..

Empirically, if a suitable limit where $`\delta 𝒲^10`$ could ever be found, would this mean the possibility of faster-than-light-communication signals? Here, we should stress that these sub-hypothetical “falico” signals would be fundamentally threshold phenomena and so a deeper understanding of this threshold between the quantum and the classical regimes could be required. In other words, the very existence of stable semiclassically entangled states would be at question. On the other hand, even if granted their existence for some special systems, in a suitable limit of nontrivial $`\delta 𝒲^1`$, it could still be a formidable, pehaps nearly impossible task to actually measure such signals, even though some direct measurements of Wigner functions (if ever needed) have already been claimed .

Philosophically, what would be at stake here is, among other things, the principle of relativity. But in this respect, it seems interesting to point out the usually overlooked fact that, as originally formulated , the relativistic partial order presupposes events which are correlated by an exchange of energy. In other words, which are correlated by “material” causes. Such is not the case for the quantum correlations treated here, which could pehaps be described as “formal” correlations. Maybe, these two kinds of correlations being essentially distinct, their principles need not be the same.

Acknowledgements: Speculations supported by CNPq via a postdoctoral fellowship to UC Berkeley. I thank Alan Weinstein, Jair Koiller, Robert Littlejohn, A. Ozorio de Almeida and Michael Berry for some pertinent, though not necessarily positive, comments (the responsability for these speculations is entirely mine). This letter is not under consideration for publication in any other way, or elsewhere.


text:
# Optical implementation of Deutsch-Jozsa and Bernstein-Vazirani quantum algorithms in eight dimensions

## Abstract

We report on a fiber-optics implementation of the Deutsch-Jozsa and Bernstein-Vazirani quantum algorithms for 8-point functions. The measured visibility of the 8-path interferometer is about 97.5%. Potential applications of our setup to quantum communication or cryptographic protocols using several qubits are discussed.

The last decade has seen the emergence of the field of quantum information processing. A particularly promising application is the concept of quantum algorithms, which allow certain problems such as factorizationShor  or searchingGrover  to be solved much faster than on a classical computer. Another algorithm which we will be interested in here is Deutsch’s algorithmDeutsch , the first quantum algorithm ever discovered, which was later generalized by Deutsch and JozsaDeutschJozsa  (DJ). The DJ algorithm discriminates between a constant or a balanced $`N`$-point binary function using one single quantum query, while a classical algorithm requires $`𝒪(N)`$ classical queries. It was later on adapted by Bernstein and Vazirani (BV) for efficiently querying a quantum databaseBV ; BS .

In the present paper, we report on a fiber-optics implementation of the DJ algorithm using standard telecom optical components and a single-photon detector. The DJ algorithm has already been implemented using NMRnmr  (see also nmr2  for a NMR implementation of the BV algorithm), table-top opticsDJoptical  (optical demonstrations of other quantum algorithms also include Grover’s algorithmKwiat ; Amsterdam ), molecular statesmolecular , and very recently using an ion traptrap . However, our setup separates from these realizations (especially that of DJoptical ) on several major aspects. First, it relies on guided optics components, which makes it unnecessary to perform a precise alignment, and it is made robust against phase fluctuations by use of an autocompensation technique. Second, although it relies on linear optics, our realization is relatively efficient in terms of used optical resources compared to standard linear optical implementations of quantum computation. The central idea of such implementations consists in representing the basis states of a $`N`$-dimensional Hilbert space by $`N`$ optical paths so that unitary transformations are obtained by chaining linear optics components that make these paths interfereReck ; CerfKwiat . Such implementations seem, however, to be inherently inefficient since the space requirement (the number of optical components) and the time requirement both grow exponentially with the number $`n`$ of qubits (with $`N=2^n`$) fn . In contrast, in our setup, the number of components is kept linear in $`n`$, while the time needed still grows exponentially. Note that any implementation of an algorithm involving an arbitrary $`2^n`$-point function (also called oracle) does in any case require exponential resources to simulate this function. Therefore, the linear optical implementation of quantum algorithms involving oracles can reasonably be made as efficient as any other implementation in this respect. For all these reasons, our experimental demonstration works with a 8-point (3-qubit) function and might probably be extended even further without fundamental difficulty, while today’s largest size optical demonstrator of the DJ algorithm involves a 4-point functionDJoptical .

Let us start by recalling the principle of the DJ algorithm. At the core of the algorithm is the oracle which computes a function $`f(x)`$, where $`x\{0,1\}^n`$ is an $`n`$ bit string, and $`f\{0,1\}`$ is a single bit. The DJ problem is to discriminate whether $`f`$ is a constant or balanced function, while querying the oracle as few times as possible. A balanced function is such that the number of $`x`$’s on which $`f(x)=0`$ is equal to the number of $`x`$’s on which $`f(x)=1`$. Classically, $`2^{n1}+1`$ queries are necessary in the worst case, whereas the DJ algorithm requires a single query as we shall see. In this algorithm, $`n`$ qubits are used, and the basis of the Hilbert space is chosen as $`|x=|x_1x_2\mathrm{}x_n`$ where $`x_i\{0,1\}`$. The quantum oracle carries out the transformation

$$|x|y\stackrel{\mathrm{oracle}}{}|x|yf(x),$$
(1)
where $`|y`$ is an ancilla qubit. By choosing $`|y=(|0|1)/\sqrt{2}`$, the action of the oracle simplifies into

$$|x\stackrel{\mathrm{oracle}}{}(1)^{f(x)}|x.$$
(2)
since $`|y`$ then remains unchanged.

The DJ algorithm starts with the system in the state $`|0=|00\mathrm{}0`$. Next, a Hadamard transform $`H`$ is applied independently on each of the $`n`$ qubit. Using the definition $`H|x=2^{n/2}_{z\{0,1\}^n}(1)^{xz}|z`$, where $`xz=_ix_iz_imod2`$ is the inner product of two $`n`$-bit strings, we see that the Hadamard transform acting on the initial state simply yields a uniform superposition of all states: $`H|0=2^{n/2}_x|x`$. This state is then sent through the oracle whereupon it becomes

$$2^{n/2}\underset{x}{}(1)^{f(x)}|x.$$
(3)
The superposition principle allows the oracle to be queried on all input values in parallel. A second Hadamard transform is then carried out to obtain the state

$$2^n\underset{x,z}{}(1)^{xz+f(x)}|z$$
(4)
which is finally measured in the $`z`$ basis. One easily deduces from Eq. (4) that when $`f`$ is constant, the probability of measuring $`|0`$ is one. In contrast, when $`f`$ is balanced, this probability is always zero, so the DJ algorithm can distinguish with certainty between these two classes of functions by querying the oracle a single time. The BV variant of this algorithm is also based on the transformation leading to Eq. (4). Suppose that the oracle is restricted to be of the form $`f_j(x)=xj`$ where $`j\{0,1\}^n`$ is an arbitrary $`n`$-bit string. The aim is to find the bit string $`j`$ with as few queries as possible. Classically one needs at least $`n`$ queries since each query provides one independent bit of information at most about $`j`$. Quantum mechanically, a single query suffices since using $`f_j(x)`$ in Eq. (4) shows that the measurement outcome is $`z=j`$ with probability.

Our all optical fiber (standard SMF-28) setup is illustrated in Fig. 1. Initially, a 3 ns light pulse is produced by a laser diode at $`1,55\mu `$m, attenuated by an optical attenuator (Agilent 8156A), and then is processed through three unbalanced Mach-Zehnder (M-Z) interferometers with path length differences $`\mathrm{\Delta }_l`$ ($`l=1,2,3`$) obeying $`\mathrm{\Delta }_3=2\mathrm{\Delta }_2=4\mathrm{\Delta }_1=15`$ ns. Each MZ interferometer doubles the number of pulses so that, at the coupler C6, we get eight equally spaced pulses. This corresponds to the action of the Hadamard transform on the three input qubits in state $`|0`$. The pulses are then reflected by a Faraday mirror and, on their way back, are modulated by a phase modulator (Trilink) commanded by a pattern generator (Agilent 81110A) which selectively puts a phase shift of $`0`$ or $`\pi `$ on each pulse according to the 8-point function \[see Eq. (3)\]. The pulses then pass back through the three M-Z interferometers, thereby realizing a second triple Hadamard transform, and are sent via a circulator to a single-photon detector (id Quantique id200), which completes the DJ or BV algorithm. The additional delay lines $`L_1`$, $`L_2`$, and $`L_3`$, which obey $`L_3>2L_2>4L_1>8\mathrm{\Delta }_3`$, ensure that the different outputs of the DJ or BV algorithm all reach the detector at different times. The delay lines $`L_2`$ and $`L_3`$ also contain an isolator so that the pulses are not transmitted on their way to the mirror. The photodetector was gated during 5 ns around the arrival time of each pulse by the pattern generator. The output of the detector was registered using a time to digital delay converter (ACAM-GP1) connected to a computer. All delays $`\mathrm{\Delta }_i`$ and $`L_i`$ were chosen to be integer multiples of $`\mathrm{\Delta }_1`$ within 0,2 ns. All electronic components were triggered by a pulse generator (Standford Research Inc. DG535). In order to maximize the visibility, polarization controllers were introduced in the long arm of each M-Z interferometer and in front of the polarization-sensitive phase modulator. Once optimized, the setup was stable for days.

This implementation of the DJ and BV algorithms differs from an earlier optical implementation of the DJ algorithmDJoptical  in several important aspects. First, we run the algorithm for $`n=3`$ qubits and, more importantly, we measure all 8 outcomes, which makes it possible to realize the BV algorithm as well (the previous implementation DJoptical  works with $`n=2`$ qubits and only measures the outcome $`z=0`$). Second, we operate at telecom wavelengths in optical fibers using a setup closely inspired from the “plug-and-play” quantum cryptographic system developed by Gisin and collaborators (see e.g. plugandplay ). For this reason, the present setup in a slightly modified version can be adapted to implement quantum cryptography using higher dimensional systemshighdQcrypto  or to illustrate quantum communication complexity protocolsBCW  over distances of a few kilometers. These potential applications will be discussed below. Third, the used resources are quite different in the two implementations: when scaled to a large $`n`$, the implementation of DJoptical  requires exponential time and exponentially many optical elements, whereas our implementation also requires exponential time but only a linear number of optical elements. This is because the $`n`$ qubits are realized as $`2^n`$ separate optical paths in DJoptical , whereas in our case they are represented as $`2^n`$ light pulses traveling in a single optical fiber, extending naturally the “time-bin” realization of qubits used in plugandplay . The small number of optical elements in our setup therefore implies that it is relatively easy to increase the number of qubits $`n`$ while keeping the optical setup stable. As we shall see, a disadvantage of our setup is that the Hadamard transform can only be implemented with a probability of success of $`1/2`$. Since $`2n`$ Hadamard transforms are needed for the DJ algorithm with $`n`$ qubits, the resulting attenuation is $`2^{2n}`$.

Let us now prove that our optical setup indeed realizes the DJ and BV algorithm. The quantum state describing the eight pulses at coupler C6 can be written as

$`|\psi {\displaystyle \underset{x=000}{\overset{111}{}}}\mathrm{exp}\left[i{\displaystyle \underset{l=1}{\overset{3}{}}}(kx_l\mathrm{\Delta }_l+\pi x_l)\right]|{\displaystyle \underset{l=1}{\overset{3}{}}}x_l\mathrm{\Delta }_l`$ (5)
where $`x`$ stands for $`(x_1,x_2,x_3)`$ and $`|p`$ denotes a pulse located at position $`p`$. The three bits $`x_1,x_2,x_3=0,1`$ label whether the pulse took the short ($`x=0`$) or the long ($`x=1`$) path through each interferometer. The factor $`\mathrm{exp}[ik_lx_l\mathrm{\Delta }_l]`$ with $`k`$ being the wave number takes into account the phase difference between a pulse traveling along the short or long paths of the interferometers. The factor $`\mathrm{exp}[i_l\pi x_l]`$ takes into account the phase accumulated at the couplers: if the pulse takes the short path, it is transmitted at two couplers, whereas, if it takes the long path it is reflected twice. After reflection at the Faraday mirror and phase modulation, the pulses cross again the three M-Z interferometers and reach the photodetector in the state

$`|\psi {\displaystyle \underset{x,y,z}{}}(1)^{f(x)}\mathrm{exp}\left[i{\displaystyle \underset{l=1}{\overset{3}{}}}\left(k(x_l+y_l)\mathrm{\Delta }_l+\pi (x_l+y_ly_l(z_l+z_{l+1})+{\displaystyle \frac{z_l+z_{l+1}}{2}})\right)\right]|{\displaystyle \underset{l=1}{\overset{3}{}}}((x_l+y_l)\mathrm{\Delta }_l+z_lL_l)`$ (6)
where the bits $`y_1,y_2,y_3`$ equal $`0`$ or $`1`$ according to whether the pulse passed through the short or the long path of each of the interferometers on its way back, and the bits $`z_1,z_2,z_3`$ equal $`0`$ or $`1`$ according to whether or not the pulse exited each of the interferometer in the path containing the delay lines $`L_1`$, $`L_2`$ or $`L_3`$. Note that we put $`z_4=0`$. We have again taken into account the phases induced by transmission or reflexion at the couplers C1-C6. The final state contains 120 pulses, but we are only interested in the eight pulses such that $`x_1+y_1=x_2+y_2=x_3+y_3=1`$, which are those that exhibit 8-path interference. The other pulses are filtered out in the computer analysis (they correspond to different time bins). The final state then becomes

$`|\psi {\displaystyle \underset{z}{}}(i)^{z_1}(1)^{z_2+z_3}{\displaystyle \underset{x}{}}(1)^{f(x_1,x_2,x_3)+x_1(z_1+z_2)+x_2(z_2+z_3)+x_3z_3}|\mathrm{\Delta }_1+\mathrm{\Delta }_2+\mathrm{\Delta }_3+z_1L_1+z_2L_2+z_3L_3`$ (7)
By relabeling the time bins according to the substitution $`z_1z_1+z_2+z_3`$ mod 2, $`z_2z_2+z_3`$ mod 2, and $`z_3z_3`$, this equation coincides (up to irrelevant phases and an overall normalization factor) with Eq. (4) with the 8 logical states $`|z`$ identified as specific time bins.

The setup thus realizes the DJ or BV algorithm, the main difference with an ideal algorithm being an extra attenuation by a factor of $`2^7`$. A factor $`2^3`$ originates from the couplers C2, C4, and C6 because each time a pulse passes through these couplers it only has a probability $`1/2`$ of exiting by the right path. Otherwise, it is absorbed by the isolators I1 and I2 or by the unconnected fiber pigtail at coupler C6. Another factor $`2^3`$ is due to the filtering out of the 112 pulses produced on the way back that do not correspond to 8-path interferences. The remaining factor $`2^1`$ is due to the coupler C7. This overall loss of 21 dB could be remedied by replacing the couplers C2, C4, C6, and C7 by optical switches which would direct the light pulses along the appropriate path. High speed, low loss optical switches are not available commercially at present, so we had to use couplers in the present experiment. However, we emphasize that this is a technological rather than a fundamental limitation.

In order to characterize the performances of our setup, we considered the $`2^n`$ oracles of the form $`f_j(x)=xj`$ and $`\overline{f}_j(x)=xj+1\mathrm{mod}2`$ (i.e., the oracles in the BV algorithm and their complements). For each oracle $`f_j`$ (or $`\overline{f}_j`$), we ran the algorithm 500,000 times and registered the number of counts in time bin $`z`$, denoted as $`N_j(z)`$ (or $`\overline{N}_j(z)`$). The algorithm gives constructive interference in the time bin $`z=j`$ for the oracle $`f_j`$ or $`\overline{f}_j`$, and destructive interference elsewhere. We then computed

$$V_j(z)=\frac{1}{2}\left(\frac{N_z(z)N_j(z)}{N_z(z)+N_j(z)}+\frac{\overline{N}_z(z)N_j(z)}{\overline{N}_z(z)+N_j(z)}\right)$$
(8)
for each pair of oracles with $`jz`$, and calculated the visibility $`V(z)`$ in time bin $`z`$ by taking the average of $`V_j(z)`$ over all values of $`j`$. The measured visibilities $`V(z)`$ for 2 and 3 qubits are shown in Table I. Remarkably, they remain relatively high when going from 2 to 3 qubits in spite of the fact that 8 path interferences are involved. This is because the path differences are automatically compensated and only $`n+1`$ polarizations must be adjusted. It should therefore be relatively easy to go beyond $`n=3`$ without significantly decreasing the visibilities.

Because of the attenuation in our setup along with the quantum efficiency ($`10.5`$ %) and the dark count probability ($`10^4`$ ns<sup>-1</sup>) of our detector, the signal-to-noise ratio was not high enough to perform these visibility measurements in the single-photon regime. For this reason, in our experiment, approximatively 20 photons per pulse entered the phase modulator (oracle) in 4 dimensions, and approximatively 50 in 8 dimensions. However, minimal modifications should allow us to decrease the number of photons while keeping the signal-to-noise ratio constant. In particular, by reducing the pulse length or using two detectors instead of the coupler C7, it should be possible to operate in the single-photon regime in 4 (and possibly 8) dimensions. Moreover, as already mentioned, the Hadamard transform could be rendered deterministic by using fast low-loss optical switches instead of couplers, which would strongly reduce the losses.

The present experiment can be extended in several ways. For example, one could implement the distributed Deutsch-Jozsa problemBCW  where two parties, which receive each as input a $`2^n`$-bit string (denoted as $`f`$ and $`g`$), must decide whether $`f=g`$ or $`f`$ differs from $`g`$ in exactly $`2^{n1}`$ bits (they are promised that only one of these two cases can occur). The value of $`n`$ for which the gap between the classical and quantum algorithms sets in is unknown, but recent results suggest that it could be as early as $`n=4`$ B . The distributed Deutsch-Jozsa algorithm could be realized with a slight modification of our setup in which the two parties, separated by by a few kilometers of optical fiber, would use each a phase modulator. This quantum protocol can also be easily adapted to multidimensional quantum cryptography, where two parties randomly choose their patterns $`f`$ and $`g`$. By publicly revealing part of $`f`$ and $`g`$, the parties can use the correlations between the measurement outcomes to establish a secret key. A detailed analysis shows that for n = 2, this exactly coincides with the 4-dimensional cryptosystem based on 2 mutually conjugate bases, which has been shown to present advantages over quantum cryptography in two dimensionsmultidimcrypto . A final potential application of this setup is to test quantum non locality using the entanglement-based Deutsch-Jozsa correlationsBCT . An entangled state of $`2^n`$ time bins must be produced, for example using the source RMZG , and each party must then carry out phase modulation and a Hadamard transform (which we have demonstrated are easy to realize on time bin entangled photons). The correlations between the chosen phases and the measured time of arrival of the photons at each side should exhibit quantum non locality. It has recently been shown that these correlations are non local for $`n=4`$ B , and that they exhibit exponentially strong resistance to detector inefficiency for large $`n`$ M , which means here that the 3 dB losses at each Hadamard transform can in principle be tolerated.

In summary, we have demonstrated, by implementing the DJ and BV algorithms for 3 qubits, a simple and robust method for manipulating multidimensional quantum information encoded in time bins in optical fibers. We anticipate that our method will have wide applicability for quantum information processing and quantum communication using higher dimensional systems.

It is a pleasure to thank N. Gisin, W. Tittel, G. Ribordy, H. Zbinden, and all the members of GAP optique for very helpful discussions and for providing the single-photon detector. We acknowledge financial support from the Communauté Française de Belgique under ARC 00/05-251, from the IUAP programme of the Belgian government under grant V-18, and from the EU under project EQUIP (IST-1999-11053).