text:
# Neutrino oscillations with three flavors in matter of varying density

## I Introduction

In previous papers , we have given analytic expressions for the three flavor neutrino oscillation evolution operator and the transition probabilities in presence of constant matter densities expressed in the vacuum mixing matrix elements and the neutrino energies or masses, i.e., incorporating the so-called Mikheyev–Smirnov–Wolfenstein (MSW) effect . Here we will discuss the application of this to realistic matter density variations in a “semi-classical” approximation based on our previous results. This allows a simple and efficient calculation of neutrino oscillations in media of varying densities. We compare this approximate formula with a numerical simulation in a multi-step model. We will as before assume that the $`CP`$ phase $`\delta `$ is equal to zero. Thus, the neutrino mixing matrix is real. The semi-classical approximation for three neutrino flavors, we believe, is a unique part of our investigation.

Previous work on models for three flavor neutrino oscillations in matter for constant matter density includes works of Barger et al. , Kim and Sze , and Zaglauer and Schwarzer . Approximate solutions for three flavor neutrino oscillations in matter have been presented by Kuo and Pantaleone  and Joshipura and Murthy . Approximate treatments have also been done by Toshev and Petcov . D’Olivo and Oteo have made contributions by using an approximative Magnus expansion for the time evolution operator . Extensive numerical investigations for matter enhanced three neutrino oscillations have been made by Fogli et al. . Studies of neutrino oscillations in Earth has been performed by several authors .

Neutrino oscillations for matter with linearly varying density have been treated by Petcov  and Lehmann et al. . Osland and Wu  have also solved the case for exponentially varying density. Matter enhanced two flavor neutrino oscillations with an arbitrary monotonic matter density profile have been studied by Balantekin and Beacom  using a uniform semi-classical approximation. See also Fishbane et al.  for two flavor neutrino oscillations in matter of varying density.

## II The evolution operator in presence of matter

Let the flavor state basis and mass eigenstate basis be denoted by $`_f\{|\nu _\alpha \}_{\alpha =e,\mu ,\tau }`$ and $`_m\{|\nu _a\}_{a=1}^3`$, respectively. Then, the flavor states $`|\nu _\alpha _f`$ can be obtained as a superpositions of the mass eigenstates $`|\nu _a_m`$, or vice versa. The bases $`_f`$ and $`_m`$ are of course just two different representations of the same Hilbert space $``$.

In the present analysis, we will use the plane wave approximation to describe neutrino oscillations. In this approximation, a neutrino flavor state $`|\nu _\alpha `$ is a linear combination of neutrino mass eigenstates $`|\nu _a`$’s such that

$$|\nu _\alpha =\underset{a=1}{\overset{3}{}}U_{\alpha a}^{}|\nu _a,$$
(1)
where $`\alpha =e,\mu ,\tau `$. In what follows, we will use the short-hand notations $`|\alpha |\nu _\alpha `$ and $`|a|\nu _a`$ for the flavor states and the mass eigenstates, respectively.

The components of a state $`\psi `$ in flavor basis and mass basis, respectively, are related to each other by

$$\psi _f=U\psi _m,$$
(2)
where

$$\psi _f\left(\psi _\alpha \right)\left(\begin{array}{c}\psi _e\\ \psi _\mu \\ \psi _\tau \end{array}\right)_f\text{and}\psi _m\left(\psi _a\right)\left(\begin{array}{c}\psi _1\\ \psi _2\\ \psi _3\end{array}\right)_m.$$

A convenient parameterization for $`U=U(\theta _1,\theta _2,\theta _3)`$ is given by

$$U=\left(\begin{array}{ccc}C_2C_3& S_3C_2& S_2\\ S_3C_1S_1S_2C_3& C_1C_3S_1S_2S_3& S_1C_2\\ S_1S_3S_2C_1C_3& S_1C_3S_2S_3C_1& C_1C_2\end{array}\right),$$
(3)
where $`S_i\mathrm{sin}\theta _i`$ and $`C_i\mathrm{cos}\theta _i`$ for $`i=1,2,3`$. This is the standard representation of the neutrino mixing matrix. The quantities $`\theta _i`$, where $`i=1,2,3`$, are the vacuum mixing angles. Since we have put the $`CP`$ phase equal to zero in the neutrino mixing matrix, this means that $`U_{\alpha a}^{}=U_{\alpha a}`$ for $`\alpha =e,\mu ,\tau `$ and $`a=1,2,3`$.

In mass basis, the Hamiltonian $``$ for the propagation of the neutrinos in vacuum is diagonal and given by

$$H_m=\left(\begin{array}{ccc}E_1& 0& 0\\ 0& E_2& 0\\ 0& 0& E_3\end{array}\right),$$
(4)
where $`E_a=\sqrt{m_a^2+𝐩^2}`$, $`a=1,2,3`$, are the energies of the neutrino mass eigenstates $`|a`$, $`a=1,2,3`$ with masses $`m_a`$, $`a=1,2,3`$. We will assume the three-momentum $`𝐩`$ to be the same for all mass eigenstates.

When neutrinos propagate in ordinary matter, there is an additional term in the Hamiltonian $``$ coming from the presence of electrons in matter . This term, the potential term, is diagonal in flavor basis and is given by

$$V_f=A\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)AK_f,$$
(5)
where

$$AA(r)=\pm \sqrt{2}G_FN_e(r)\pm \frac{1}{\sqrt{2}}G_F\frac{1}{m_N}\rho (r)$$

is the matter density parameter and $`K_f`$ is the projector in flavor basis on the electron neutrinos. Here $`G_F`$ is the Fermi weak coupling constant, $`N_e`$ is the electron density, $`m_N`$ is the nucleon mass, and $`\rho `$ is the matter density. The sign of the matter density parameter depends on weather we deal with neutrinos ($`+`$) or antineutrinos ($``$). In mass basis, this piece of the Hamiltonian is $`V_m=U^1V_fU`$, where $`U`$ is again the neutrino mixing matrix.

In the case when the neutrinos propagate through matter, as here, the Hamiltonian is not diagonal in either the mass basis or the flavor basis, and we have to calculate the evolution operator $`U_f(t)`$ or $`U_f(L)e^{i_fL}=Ue^{i_mL}U^1`$ if we set $`t=L`$ ($`L`$ is the traveling (propagation) path length of the neutrinos.).

To do so it is convenient to introduce the traceless real symmetric matrix $`T`$ defined by $`T_m(\mathrm{tr}_m)I/3`$. The trace of the Hamiltonian in mass basis $`_mH_m+U^1V_fU`$ is $`\mathrm{tr}_m=E_1+E_2+E_3+A`$, and the matrix $`T`$ can then be written as

$$T=(T_{ab})=\left(\begin{array}{ccc}AU_{e1}^2\frac{1}{3}A+\frac{1}{3}\left(E_{12}+E_{13}\right)& AU_{e1}U_{e2}& AU_{e1}U_{e3}\\ AU_{e1}U_{e2}& AU_{e2}^2\frac{1}{3}A+\frac{1}{3}\left(E_{21}+E_{23}\right)& AU_{e2}U_{e3}\\ AU_{e1}U_{e3}& AU_{e2}U_{e3}& AU_{e3}^2\frac{1}{3}A+\frac{1}{3}\left(E_{31}+E_{32}\right)\end{array}\right),$$
(6)
where $`E_{ab}E_aE_b`$. Of the six antisymmetric quantities $`E_{ab}`$, where $`a,b=1,2,3`$ and $`ab`$, only two are linearly independent, since the $`E_{ab}`$’s fulfill the relations $`E_{ba}=E_{ab}`$ and $`E_{12}+E_{23}+E_{31}=0`$.Later, we will use the usual (vacuum) mass squared differences $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{\Delta }m_{32}^2`$, instead of $`E_{21}`$ and $`E_{32}`$, which are related to each other by $`\mathrm{\Delta }m_{21}^2=2E_\nu E_{21}`$ and $`\mathrm{\Delta }m_{32}^22E_\nu E_{32}`$, where $`E_\nu `$ is the neutrino energy. This means that the evolution operator in mass basis can be written as

$$U_m(L)e^{i_mL}=\varphi e^{iLT}=\varphi \underset{a=1}{\overset{3}{}}e^{iL\lambda _a}\frac{1}{3\lambda _a^2+c_1}\left[(\lambda _a^2+c_1)I+\lambda _aT+T^2\right],$$
(7)
where $`\varphi e^{iL(\mathrm{tr}_m)I/3}`$, $`\lambda _a`$, $`a=1,2,3`$, are the eigenvalues of the matrix $`T`$, and $`I`$ is the $`3\times 3`$ identity matrix.§§§Using Cayley–Hamilton’s theorem, the exponential of a matrix $`M`$, $`e^M=_{n=0}^{\mathrm{}}\frac{1}{n!}M^n`$ (infinite series), can be written as $`e^M=_{n=0}^{N1}a_nM^n`$ (finite series, $`N`$ is the dimension of $`M`$), where $`a_n`$ ($`n=0,1,\mathrm{},N1`$) are some coefficients to be determined. In this case, since $`N=3`$ (the dimension of $`T`$ is three), this means that we have $`e^{iLT}=a_0Iia_1LTa_2L^2T^2`$, i.e., there are no higher power terms of $`T`$ in $`e^{iLT}`$ than that of order two. The coefficients $`c_0`$, $`c_1`$, and $`c_2`$ are all real and the eigenvalues $`\lambda _a`$, $`a=1,2,3`$, can be expressed in closed form in terms of these .

The evolution operator for the neutrinos in flavor basis is thus given by

$$U_f(L)=e^{i_fL}=Ue^{i_mL}U^1=\varphi \underset{a=1}{\overset{3}{}}e^{iL\lambda _a}\frac{1}{3\lambda _a^2+c_1}\left[(\lambda _a^2+c_1)I+\lambda _a\stackrel{~}{T}+\stackrel{~}{T}^2\right],$$
(8)
where $`\stackrel{~}{T}UTU^1`$. Equation (8) is our final expression for $`U_f(L)`$.

Since $`_f=U_mU^1`$, it is clear that $`\stackrel{~}{T}=_f(\mathrm{tr}_f)I/3=_f(\mathrm{tr}_m)I/3`$ due to the invariance of the trace under transformation of $`U`$. In fact, the characteristic equation is also invariant under transformation of $`U`$ and therefore so are the coefficients $`c_0`$, $`c_1`$, $`c_2`$, and the eigenvalues $`\lambda _1`$, $`\lambda _2`$, $`\lambda _3`$. However, the expression for $`_f`$ is much more complicated than that for $`_m`$, which is the reason why we work with $`_m`$ instead of $`_f`$.

The formula (8) expresses the time (or $`L`$) evolution directly in terms of the mass squared differences and the vacuum mixing angles without introducing any auxiliary matter mixing angles. By dividing the density variation in small, approximately constant segments, and using this formula repeatedly in each segment, we can numerically study neutrino oscillations in matter with varying density. We will use this method as a standard test for the semi-classical approximation to the evolution operator that we study below.

## III The semi-classical approximation

The formula given in Eq. (8) is the evolution operator in flavor basis for constant matter density. To handle the case of varying matter density, let us divide the distance $`L`$ from the source to the detector into $`N`$ equidistant parts and put the index $`k`$ on the eigenvalues $`\lambda _a^k`$, where $`a=1,2,3`$ denote the three mass eigenstates. For any matter density profile $`\rho (r)`$ we first make this profile discrete and introduce $`\rho _k`$ as the matter density in the interval $`r_{k1}rr_k`$, where $`k`$ varies from $`1`$ to $`N`$ with $`r_0=0`$ and $`r_N=L`$. The length of each segment is then $`\mathrm{\Delta }r_k=r_kr_{k1}=L/N`$. The evolution operator in mass basis from $`0`$ to $`L`$ can then be written as the ordered product

$$U_m(L)=U_m(r_Nr_{N1})U_m(r_{N1}r_{N2})\mathrm{}U_m(r_2r_1)U_m(r_1r_0).$$
(9)
Note that the order of the $`U_m(r_kr_{k1})`$’s are important, since these operators do not in general commute.

When $`N`$ is large, each step is small and the exponent in the evolution operator $`U_m(r_kr_{k1})=U_m(\mathrm{\Delta }r_k)`$ is small. We can then approximate this operator with

$$U_m(\mathrm{\Delta }r_k)e^{i\mathrm{\Delta }r_k_m^k},$$
(10)
where $`_m^kH_m+A_kK_m`$ and $`A_k\rho _k`$. Inserting this into Eq. (9) gives

$$U_m(L)e^{i\mathrm{\Delta }r_N_m^N}e^{i\mathrm{\Delta }r_{N1}_m^{N1}}\mathrm{}e^{i\mathrm{\Delta }r_2_m^2}e^{i\mathrm{\Delta }r_1_m^1}.$$
(11)
Since the $`_m^k`$’s do not commute, the higher order terms have to be calculated with the time-ordering (here rather $`r`$-ordering) operator. However, here we will at first be satisfied with the lowest order result, which we call the semi-classical approximation. In this approximation, we retain only the terms proportional to $`\mathrm{\Delta }r_k=L/N`$, and thus, neglect the noncommutativity of the $`_m^k`$’s for different $`k`$’s. We can thus write

$$U_m(L)e^{i_{k=1}^N\frac{L}{N}_m^k}.$$
(12)
In the limit $`N\mathrm{}`$, this gives the integral formula

$$U_m(L)=e^{i_0^L_m(r)𝑑r}=\varphi (L)e^{i_0^LT(r)𝑑r},$$
(13)
where $`T(r)`$ is the traceless part of the Hamiltonian corresponding to the electron density at position $`r`$ between $`0`$ and $`L`$ and $`\varphi (L)`$ is the phase factor coming from the trace.

For further discussion it is often convenient to retain the original form of the Hamiltonian and to use $`_m=H_m+AU^1K_fU`$ rather than $`T`$. Thus, when $`A=A(r)`$, we obtain

$$_0^L_m(r)𝑑r=L(H_m+\overline{A}(L)K_m),$$
(14)
where

$$\overline{A}(L)\frac{1}{L}_0^LA(r)𝑑r$$

is the average matter density along the baseline $`L`$ and

$$K_mU^1K_fU,$$

which means that the evolution operator in mass basis can be written as

$$U_m(L)=e^{iL\left(H_m+\overline{A}(L)K_m\right)}\overline{\varphi }e^{iL\overline{T}},$$
(15)
where $`\overline{\varphi }e^{iL(\mathrm{tr}\overline{}_m)/3}`$, $`\overline{T}\overline{}_m(\mathrm{tr}\overline{}_m)I/3`$, and $`\overline{}_mH_m+\overline{A}(L)K_m`$.

For this case we can thus use the previous expression for $`T=_m(\mathrm{tr}_m)I/3`$ by simply replacing $`A`$ with $`\overline{A}(r)`$ and then pass to flavor basis by using the $`U`$ transformation, i.e., $`\stackrel{~}{T}=UTU^1`$.

Thus, for any $`L`$ we can use the spectral decomposition theorem and we find that

$$U_f(L)=\varphi (L)\underset{a=1}{\overset{3}{}}e^{iL\lambda _a(L)}P_a(L),$$
(16)
where $`\lambda _a(L)`$ is the $`a`$th eigenvalue of $`T(L)`$ (or $`\stackrel{~}{T}(L)`$) and

$$P_a(L)=\frac{1}{3\lambda _a^2(L)+c_1(L)}\left[(\lambda _a^2(L)+c_1(L))I+\lambda _a(L)\stackrel{~}{T}(L)+\stackrel{~}{T}^2(L)\right]$$
(17)
is the projection operator. Everything here is of course $`L`$-dependent, since the operator is $`L`$-dependent and therefore also the eigenvalues. The phase factor is $`\varphi (L)=e^{iL(\mathrm{tr}_f(L))/3}=e^{iL(\mathrm{tr}_m(L))/3}`$.

In the case of a linear matter density of the form $`A(r)=A+Br`$, we obtain $`\overline{A}(r)=A+Br/2`$. Similarly, in the case of a step function like matter density, relevant to the matter distribution of the Earth, we have

$$A(r)=\{\begin{array}{cc}A_1,\hfill & 0rL_1,\hfill \\ A_2,\hfill & L_1rL_1+L_2,\hfill \\ A_1,\hfill & L_1+L_2r2L_1+L_2,\hfill \end{array},\text{where }2L_1+L_2L\text{,}$$

which leads to

$$\overline{A}(r)=\{\begin{array}{cc}A_1,\hfill & 0rL_1,\hfill \\ A_2\left(1\frac{L_1}{r}\right)+A_1\frac{L_1}{r},\hfill & L_1rL_1+L_2,\hfill \\ A_1\left(1\frac{L_2}{r}\right)+A_2\frac{L_2}{r},\hfill & L_1+L_2r2L_1+L_2,\hfill \end{array}.$$
(18)
Finally, in the case of an exponentially decreasing matter density $`A(r)=Ae^{r/r_0}`$, where $`A`$ and $`r_0`$ are parameters relevant to the matter distribution of the Sun, we obtain

$$\overline{A}(r)=A\frac{r_0}{r}\left(1e^{r/r_0}\right).$$
(19)
We can see here that in the semi-classical approximation the influence of the density $`\overline{A}`$ decays as $`1/L`$ with distance $`L`$ from the matter. The evolution should therefore be continued with the vacuum evolution operator as soon as the neutrinos leave the matter region.

## IV Digression on the semi-classical approximation

Let us consider the semi-classical (s.c.) evolution operator further. It can be written as

$$U(r)_m^{\mathrm{s}.\mathrm{c}.}=e^{i_m^{(1)}(r)},$$
(20)
where

$$_m^{(1)}(r)H_mr+A^{(1)}(r)K_m.$$

Here $`A^{(1)}(r)_0^rA(s)𝑑s`$. Now, the equation of motion for the full evolution operator $`U_m`$ is

$$i\frac{d}{dr}U_m(r)=_m(r)U_m(r).$$
(21)
This can be integrated to give the equation

$$U_m(r)=1i_0^r_m(s)U_m(s)𝑑s.$$
(22)
Upon differentiating the semi-classical evolution operator above, we see that, although

$$\frac{d}{dr}_m^{(1)}(r)=_m(r)=H_m+A(r)K_m,$$
(23)
we can equate $`i\frac{d}{dr}U_m^{\mathrm{s}.\mathrm{c}.}(r)`$ with $`_m(r)U_m^{\mathrm{s}.\mathrm{c}.}(r)`$ only when the commutator $`[_m^{(1)}(r),_m(r)]`$ can be neglected. This commutator can be calculated to be

$$[_m^{(1)}(r),_m(r)]=_0^rs\frac{dA}{ds}(s)𝑑s[H_m,K_m].$$
(24)
Thus, when $`_0^rs\frac{dA}{ds}(s)𝑑s=\frac{r^2}{2}\frac{dA}{ds}(\xi )`$ for $`0\xi r`$ is small, the semi-classical approximation to the evolution operator is a good approximation to the full evolution operator. For constant matter density this is of course true. For linear matter density $`A(r)=A+Br`$ the coefficient $`B`$ should be small, i.e., $`BrA`$ or at least $`Br<A`$.

Equation (22) can be solved in a systematic way by iteration, leading to

$$U_m(r)=1i_0^r_m(s)𝑑s+(i)^2_0^r_m(s)_0^s_m(s^{})𝑑s^{}𝑑s+\mathrm{}.$$
(25)
The Hamiltonian $`_m(r)`$ is the one given in Eq. (23). The result to second order is then given by

$$U_m(r)1ir\left(H_m+\overline{A}(r)K_m\right)+(i)^2\frac{r^2}{2}\left(H_m+\overline{A}(r)K_m\right)^2+(i)^2\frac{r^2}{2}\left(\overline{\overline{A}}(r)\overline{A}(r)\right)[H_m,K_m],$$
(26)
where

$$\overline{A}(r)A^{(1)}(r)/r,A^{(1)}(r)_0^rA(s)𝑑s,\overline{\overline{A}}(r)2A^{(2)}(r)/r^2,A^{(2)}(r)_0^rA^{(1)}(s)𝑑s.$$

By inspection we see that the expression in Eq. (26) deviates from an expansion of the semi-classical approximation by the terms proportional to the commutator $`[H_m,K_m]`$ and higher order terms in $`H_m`$ and $`K_m`$. In fact, the commutator between the Hamiltonian at different points $`s`$ and $`s^{}`$ is

$$[_m(s),_m(s^{})]=(A(s)A(s^{}))[H_m,K_m]=\frac{dA}{ds}|_{s=s^{}}(ss^{})[H_m,K_m]+\mathrm{},$$
(27)
which vanishes only for $`A(s)=A(s^{})`$. This is in general true only for constant matter densities, $`A(s)=A`$. When $`\frac{dA}{ds}`$ is large, the contribution of the commutator cannot be neglected.

We can therefore sum the semi-classical approximation terms and write the solution as

$$U_m(r)=U_m^{\mathrm{s}.\mathrm{c}.}(r)+a_1[H_m,K_m]+\mathrm{},$$
(28)
where

$$a_1=(i)^2\frac{r^2}{2}\left(\overline{\overline{A}}(r)\overline{A}(r)\right)$$

and the dots represent higher order terms that vanish when the commutator $`[H_m,K_m]`$ is neglected. When $`\left|\overline{\overline{A}}\overline{A}\right|`$ is small, the correction terms are small.

## V Probability amplitudes and transition probabilities

In the previous sections, we have calculated the evolution operator in the semi-classical approximation. Below we will study the corresponding probability amplitudes and transition probabilities.

The probability amplitude $`A_{\alpha \beta }`$ for $`\nu _\alpha \nu _\beta `$ transition is simply defined as the $`(\beta ,\alpha )`$-matrix element of the evolution operator in flavor basis, i.e.,

$$A_{\alpha \beta }\beta |U_f(L)|\alpha ,\alpha ,\beta =e,\mu ,\tau .$$
(29)
We now consider transition probabilities for neutrino oscillations in the semi-classical approximation given by Eq. (16). Inserting Eq. (16) into Eq. (29) gives

$$A_{\alpha \beta }=\varphi (L)\underset{a=1}{\overset{3}{}}e^{iL\lambda _a(L)}P_a(L)_{\beta \alpha }$$
(30)
where

$$P_a(L)_{\beta \alpha }=\frac{(\lambda _a^2+c_1)\delta _{\beta \alpha }+\lambda _a\stackrel{~}{T}_{\beta \alpha }+(\stackrel{~}{T}^2)_{\beta \alpha }}{3\lambda _a^2+c_1}$$
(31)
is the matrix element of the projector $`P_a(L)`$. Here $`\delta _{\alpha \beta }`$ is Kronecker’s delta. Note that $`\stackrel{~}{T}_{\alpha \beta }=\stackrel{~}{T}_{\beta \alpha }`$ and $`(\stackrel{~}{T}^2)_{\alpha \beta }=(\stackrel{~}{T}^2)_{\beta \alpha }`$. The transition probability $`P_{\alpha \beta }`$ for $`\nu _\alpha \nu _\beta `$ transition is defined as the absolute value squared of the probability amplitude $`A_{\alpha \beta }`$. Hence, the transition probabilities in matter are given by the formulas

$$P_{\alpha \beta }=|A_{\alpha \beta }|^2=\delta _{\alpha \beta }4\underset{a<b}{\underset{a=1}{\overset{3}{}}\underset{b=1}{\overset{3}{}}}P_a(L)_{\beta \alpha }P_b(L)_{\beta \alpha }\mathrm{sin}^2\stackrel{~}{x}_{ab},\alpha ,\beta =e,\mu ,\tau ,$$
(32)
where $`\stackrel{~}{x}_{ab}(\lambda _a(L)\lambda _b(L))L/2`$.

## VI Applications and Discussion

The main results of our analysis are given by the evolution operator for the neutrinos when passing through matter with varying matter density in the “semi-classical approximation”, Eq. (16), and the corresponding expressions for the transition amplitudes in Eq. (30) and the transition probabilities in Eq. (32), all expressed as finite sums of simple functions in the matrix elements of $`_f`$ (or $`_m`$) and integrals involving $`A(r)`$, the varying matter density.

As applications, we have calculated the transition probability $`P_{\mu e}`$ for neutrino oscillations for different matter density profiles. Our calculations compare two different cases:

1. An “exact” numerical evolution operator method based on the product of $`N`$ evolution operators using

$$U_f(L)=U_f(r_Nr_{N1})U_f(r_{N1}r_{N2})\mathrm{}U_f(r_2r_1)U_f(r_1r_0),$$
(33)
with the formula (8) used for each step with the appropriate matter density; and
2. The semi-classical approximation method based on one single evolution operator using

$$U_f(L)=\overline{\varphi }Ue^{iL\overline{T}}U^1.$$
(34)

In all our examples discussed below, we have used the Earth center crossing neutrino traveling path length, except for the last example in which we discuss the Sun.

Let $`R_{}6371`$ km be the radius of the Earth and $`r_{}3486`$ km be the radius of the core with this approximation and in a numerical simulation based on Eq. (9) with a step length of $`L/N=2R_{}/100=127.42`$ km. The thickness of the mantle is then $`R_{}r_{}2885`$ km with the matter density parameter $`A_1=A_{\mathrm{mantle}}1.7010^{13}`$ eV ($`\rho _1=\rho _{\mathrm{mantle}}4.5`$ $`\mathrm{g}/\mathrm{cm}^3`$), whereas the matter density parameter of the core is $`A_2=A_{\mathrm{core}}4.3510^{13}`$ eV ($`\rho _2=\rho _{\mathrm{core}}11.5`$ $`\mathrm{g}/\mathrm{cm}^3`$).

Neutrinos traversing the Earth towards a detector close to the surface of the Earth, pass through the matter of varying density densities $`A(r)`$ where the distances $`L_i`$, $`i=1,2`$, are functions of the nadir angle $`h`$, where $`h\pi \theta _z`$; $`\theta _z`$ being the zenith angle. As $`h`$ varies from $`0`$ to $`\pi /2`$, the cord $`L=L(h)`$ of the neutrino passage through Earth becomes shorter and shorter. At an angle larger than $`h_0=\mathrm{arcsin}(r_{}/R_{})33.17^{}`$, the distance $`L_2=0`$, and the neutrinos no longer traverse the core.

The mass squared differences ($`\mathrm{\Delta }M^2\mathrm{\Delta }m_{32}^2`$ and $`\mathrm{\Delta }m^2\mathrm{\Delta }m_{21}^2`$) and the vacuum mixing angles ($`\theta _1,\theta _2,\theta _3`$) used here are chosen to correspond to those obtained from analyses of various neutrino oscillation data. We have taken

$$\mathrm{\Delta }M^2=3.210^3\mathrm{eV}^2,\mathrm{\Delta }m^2=0,\theta _1=45^{},\theta _2=5^{},\theta _3=45^{}.$$

The values of $`\mathrm{\Delta }M^2`$ and $`\theta _1`$ are governed by atmospheric neutrino data  and the values of $`\mathrm{\Delta }m^2`$ and $`\theta _3`$ (LMA) by solar neutrino data , where LMA stands for large mixing angle (matter) solution. The value of $`\theta _2`$ is below the CHOOZ upper bound, which is $`\mathrm{sin}^22\theta _2=0.10`$ . These choices are the most optimistic ones for obtaining any effects in long baseline (LBL) experiments from the sub-leading $`\mathrm{\Delta }m^2`$ scale . We should mention though, that these data are taken from two neutrino flavor model analyses.

As a first example, we have investigated the Earth’s matter density profile using the published Stacey model for the Earth’s matter density profile . The resulting curves of the (exact) numerical evolution operator method for a mantle-core-mantle step function approximation of the Earth’s matter density profile, the semi-classical approximation method, and for reference the (exact) numerical evolution operator method of the Earth’s matter density profile are shown in Fig. 1. It has earlier been found by Freund and Ohlsson  that a mantle-core-mantle step function approximation of the Earth’s matter density profile is a good approximation (even in a three neutrino scenario).

The numerical evolution operator method results were carried out using $`N=100`$, i.e., they consist of a product of 100 evolutions with different constant matter densities for each evolution step, whereas the semi-classical approximation method result was obtained with just one single evolution with the average matter density of the Earth’s matter density profile $`\overline{A}_{}`$ ($`\overline{\rho }_{}7.8\mathrm{g}/\mathrm{cm}^3`$, which is also the reason why the semi-classical approximation curve only has got one resonance peak at $`E_\nu =\overline{E}_{\nu ,}5.410^9\mathrm{eV}`$. This peak of course lies inbetween the both resonance peaks (ideally at $`E_\nu =E_{\nu ,\mathrm{core}}3.710^9\mathrm{eV}`$ and $`E_\nu =E_{\nu ,\mathrm{mantle}}9.410^9\mathrm{eV}`$) of the exact numerical evolution calculation, since $`A_{\mathrm{mantle}}\overline{A}_{}A_{\mathrm{core}}`$.

In Figs. 2 and 3, we have used two different step function matter density profiles.

There appear no double peaks in these figures even though the step function matter density profiles consist of two different $`A_k`$’s. Furthermore, in Figs. 2 and 3, the differences between the values of the $`A_1`$’s and $`A_2`$’s are the same. The step function matter density profile in Fig. 2 could simulate the Earth’s matter density profile if the Earth has a core, which is much less dense than has been found by geophysics. Note that the absolute error between the two curves in Fig. 2 is larger than in Fig. 3, whereas the relative error of the curves in Fig. 3 is larger than in Fig. 2.

Next, in Figs. 4 and 5, we have studied constant matter density profiles.

The constant matter density used in Fig. 4 is the average density of the Earth. The curves in this figure could be compared with the dotted curve in Fig. 1. In Fig. 5, the constant matter density was chosen to be equal to the average density in the mantle of the Earth. In both these figures, the semi-classical approximation method gives an excellent agreement with the exact numerical evolution operator method as it should, since in the case of constant matter density $`\overline{A}=A=\mathrm{const}.`$, which means that the two methods are equivalent. Thus, the very small deviations seen in the figures are only due to numerics.

Then, in Figs. 6 and 7, we discuss linear matter density profiles.

The parameter $`B`$ used in Fig. 6 is larger than that in Fig. 7. The smaller $`B`$ is the better the agreement between the numerical evolution operator method and the semi-classical approximation method becomes. However, linear matter density profiles only have theoretical interest, since they are not to be found in Nature at least at large distance scales. They could, however, be used on shorter distances scales though, e.g., LBL experiments like K2K, MINOS, and CERN-LNGS , where the neutrinos traverse the Earth mantle with cord-like paths.

Finally, in Fig. 8, we investigated the Sun’s matter density profile, which is an exponentially decreasing matter density profile. The semi-classical approximation method does not work as well in this case as for step function, constant, and linear matter density profiles, since this matter density profile is varying too quickly.A remark about the neutrino energy interval shown in Fig. 8 is in place. Solar neutrinos have energies of the order of 0.1 MeV - 10 MeV, i.e., energies much smaller than what is shown in Fig. 8. However, the resonance energy due to the large mass squared difference $`\mathrm{\Delta }M^2=3.210^3\mathrm{eV}^2`$ is given in Fig. 8. Thus, the neutrino energy interval shown in Fig. 8 is only of theoretical interest and the figure shows how well the semi-classical approximation works in the region where the transition probability $`P_{\mu e}`$ changes most.

In conclusion, the semi-classical approximation will be a good approximation for some types of matter density profiles. In certain cases, of slowly varying matter density, it is even an excellent approximation. The major advantage of the semi-classical approximation method as compared to the exact numerical evolution operator method is that we only need to calculate one single evolution operator for one single average matter density, the average matter density parameter of the considered matter density profile $`\overline{A}(L)`$, i.e., we can make the replacement

$$U_f(L)=\underset{n\mathrm{operators}}{\underset{}{\underset{i=1}{\overset{n}{}}U_f(L_i,A_i)}}U_f(L)=\underset{\mathrm{one}\mathrm{operator}}{\underset{}{\overline{\varphi }e^{iL\overline{T}}}}.$$

###### Acknowledgements.

We would like to thank Martin Freund for useful comments. This work was supported by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT) \[T.O.\], the Wenner-Gren Foundations \[T.O.\], the “Sonderforschungsbereich 375 für Astro-Teilchenphysik der Deutschen Forschungsgemeinschaft” \[T.O.\], and the Swedish Natural Science Research Council (NFR), Contract No. F 650-19981428/2000 \[H.S.\].


format: text

text:
# Simulation of Einstein-Podolsky-Rosen experiments in a local hidden variables model with limited efficiency and coherence.

\[

## Abstract

We simulate correlation measurements of entangled photons numerically. The model employed is strictly local. The correlation is determined by its classical expression with two decisive difference: we sum up coincidences for each pair individually, and we include the effect of polarizer beam splitters. We analyze the effects of decoherence, detector efficiency and polarizer thresholds in detail. The Bell inequalities are violated in these simulations. The violation depends crucially on the threshold of the polarizer switches and can reach a value of 2.0 in the limiting case. Existing experiments can be fully accounted for by limited coherence and non-ideal detector switches. It seems thus safe to conclude that the Bell inequalities are no suitable criterium to decide on the nonlocality issue.

03.65.Ud,03.65.Ta

\]

The Einstein-Podolsky-Rosen (EPR) problem has long occupied a central place in the understanding of quantum mechanics . Bell’s inequalities in conjunction with correlation measurements seemed to prove that reality in microphysics is manifestly nonlocal . Furthermore, the experimental evidence seems to contradict even the notion of an independent reality . Both of these features, if true, are highly problematic. The former, because no field propagating with a velocity exceeding c has ever been observed. The latter, because without an independent reality there is no guarantee that theoretical models can at all be contradicted. And without the possibility of contradiction scientific progress follows no clear rules.

For these reasons the EPR problem is far more important than the experiments alone indicate. Consequently, a large amount of work has been devoted to this field. Two years ago, the standard reference on EPR - the book by Afriat and Selleri  \- more or less highlighted the dilemma. But in the same year Deutsch and Hayden  could show, by an analysis of the information flow in such an experiment, that there is in fact no nonlocal connection between the two measuring devices. All information about the two angles of polarization, $`\varphi _1`$ and $`\varphi _2`$, is stored locally. Even though this information cannot locally be accessed. It is probably due to this new field of research, quantum information theory, that the problem is even more important today than it was ten years ago. Consequently, a number of papers in the last two years have analyzed the paradox from different angles, and the analysis brought two features into focus: (i) The validity and significance of Bell’s inequalities ; and (ii) the relevance of a photon’s phase for the correlations . From the viewpoint of information flow a violation of Bell’s inequalities is no proof of nonlocality . From a formal point of view it seems that the standard inequalities cannot be derived without violating established notions about the measurement process. For a detailed discussion see Sica or Adenier . The notion of a phase seemed initially problematic because the phase e.g. of a wavefunction cannot be related to physical properties of a photon. But as shown later, the phase indicates the phase of a photon’s electromagnetic field . And it could be established that the existence of a phase connection between the two points of measurement, a connection which arises at the process of emission from a common source, is sufficient to explain correlations between two measurements in space-like separation. It was also emphasized that measurements cannot in general be factorized without loosing the linearity of the fields between the two polarizers.

In this Letter we pursue a different strategy. We perform numerical simulations of actual experiments. We sought to include the features of the experimental situation as far as possible. For this reason we shall give not only the results of ideal measurements, but also measurements with limited efficiency, limited coherence, and under the condition of dead angles of our polarizers. It will be seen that all these effects have a bearing on the actual data.

Setup. \- The setup of the experiment is shown in Fig. 1 (a). A source between two polarizers emits a pair of photons along the $`z`$ axis. Two polarizers, at the positions $`L_1`$ and $`L_2`$, respectively, measure the angle of polarization. The angle of polarizer one is varied by a half cycle, $`\pi `$, during the experiment. At every position of the polarizer a set of 1000 photon pairs is emitted and measured. The switch of the polarizer is shown in Fig. 1 (b). If $`\mathrm{cos}^2(\varphi _1\alpha )`$ is larger than $`0.5+\mathrm{\Delta }s`$, then the event is recorded as a transmission (+). If it is less than $`0.5\mathrm{\Delta }s`$, it is recorded as an adsorption (-). No action is taken for values between these two boundaries. The same switch is applied to both measurement devices. The threshold $`\mathrm{\Delta }s`$ formalizes three separate features of the experiments: (i) The transmission characteristics of the polarizer beam splitters; (ii) the threshold of the photodetectors; and (iii) the electronic evaluation of events, since double counts ($`+`$ and $``$) are excluded.

A single simulation run starts with the initialization of the random number generator . The generator is initialized only once, at the beginning of a full simulation cycle. The random number is mapped onto the initial phase from 0 to $`2\pi `$ of the photon pair. Simulations are generally made with a phase difference of $`\pi /2`$ between the angles of polarization of photon one and photon two. After covering the distance to polarizer one and two the photons are measured. We assumed, without lack of generality, that both distances are integer multiples of the wavelength. After a single pair has been measured, we record the coincidences (++,+-,-+,–). The procedure is repeated for all 1000 pairs, then the polarization angle of device one is changed by $`\pi /100`$. A run ends, when all 1000 pairs at the final position of polarizer one have been measured ($`\pi `$). In all figures we only plot the coincidences $`N_{++}`$.

We accounted for limited efficiency and decoherence in the following way. Limited efficiency means that not all pairs emitted are actually measured. In this case we simply did not evaluate all pairs, depending on the efficiency of the setup. 50 % efficiency, for example, means that only every second pair is actually recorded. To simulate decoherence we created an independent random input for a certain fraction of a half cycle of $`\pi `$. 100 % decoherence here means that half a wavelength of the photon’s optical path is random. This translates into a polarization angle random in the interval \[0,$`\pi `$\]. Both effects reduce the maximum of the output measured, but it will be seen presently that they have very different effects.

Ideal measurements. \- Initially we simulated an ideal measurement. The efficiency in this simulation is 100%, the fields of the two photons are fully coherent throughout the distance between the two polarizers, and the experimental devices are supposed to have ideal characteristics. The result of this simulation is shown in Fig. 2. We did two separate simulations, one with a polarization difference between the two photons of 0 (full squares), the other with $`\pi /2`$ (full circles). It can be seen that neither of the simulations comes close to the theoretical prediction of a $`sin^2(\beta \alpha )`$. Instead, the curves representing ideal measurements would be of angular shape. However, the maximum of the correlation ($`N_0/2`$, where $`N_0`$ is the number of pairs) and the minimum (zero) are exactly obtained in the extreme cases.

It should be noted that the results given in these plots reflect the ”classical” formulation of a coincidence, given by the equation :

$$P(\alpha ,\beta )=𝑑\lambda \mathrm{cos}^2(\lambda \alpha )\mathrm{cos}^2(\lambda \beta )$$
(1)
with two decisive differences: first, the summation is performed over single pairs, as in the actual experiments, rather than over the two polarizers separately. The latter procedure, formalized in the given integral, includes not only photons of one pair, but also sums up contributions of different pairs. And second, we accounted for the digital output of the polarization devices. The integral is no suitable representation of the digital results in current experiments. If, for example, $`\lambda =\varphi _1`$ and $`\varphi _2=\varphi _1+\pi /2`$, then for $`\alpha =\beta =0`$ the integral yields $`N_0/8`$. But the actual count, under the condition that $`\mathrm{cos}^2\varphi _1>0.5`$ and $`\mathrm{cos}^2(\varphi _1+\pi /2)>0.5`$ is zero. The digitalization, necessary to obtain formal agreement with spin measurements in quantum mechanics, is usually achieved by means of a polarizer beam splitter .

Dead angles. \- In our simulations we find that the curves obtained are not very sensitive to the threshold of the polarizer switches. We have performed simulations where $`\mathrm{\Delta }s`$ was varied from 0.00 to 0.20. Apart from a reduction of the absolute yield the increase of the threshold only affects the width of the minimum at the ultimate angles. This effect is equal to a retardation of the onset of the correlation function at its minimum position. The threshold therefore does not change the functional form of the correlations.

Efficiency. \- The detection of photons is one of the problems experimenters are still confronted with. The efficiency is in fact so low (less than 10% ), that the correlations found in Aspect’s measurements  were disputed on the grounds of an ”efficiency loophole”. In our simulations such a conclusion would not be justified. Even though the shape of the curve changes somewhat and the statistical spread is dramatically increased in the low efficiency range, the maxima and minima are still clearly distinguishable. The minimum, moreover, remains zero. We included a dead angle of detectors by a threshold of 0.05. The same threshold will be used in subsequent simulations. The coincidence rates due to detector efficiency varying from 50% to 10% are shown in Fig. 3. From this result we conclude that efficiency is not the decisive issue to estimate the relevance of an experiment.

Decoherence. \- The polarizers in current measurements are more than 400 m apart . Furthermore, there is no vibration damping or cooling to very low temperatures involved in such a measurement. This feature of the measurements is bound to cause random motion of system components. From surface science the range of motion without damping can be estimated, it should be for an isolated surface no less than a few nanometer or more than one percent of the photon’s wavelength. Considering that we deal with three coupled components and optical paths in between it seems safe to increase this estimate by one order of magnitude. In this case we have to include random motion of our system in the range of about 5-10% of the wavelength. This translates, in our simulation, into a rate of decoherence of 10-20% (100% means that half a wavelength of the photon’s optical path is random).

Simulation with a decoherence rate of 10%, 50% and 100% are shown in Fig. 4. The interesting feature of decoherence is that it renders the resulting distribution more sinusoidal than the correlations of an ideal measurement. The fully decoherent simulation proves that correlations are independent of the setting as required, but it also shows the noise due to a random distribution of the initial phase of the coupled system. In practice all effects analyzed will to a greater or lesser extent be present in any single measurement.

Bell violations. \- Finally, we demonstrate the influence of the polarizer threshold on the violation of Bell’s inequalities. To this end we simulate the counts at four selected angular positions of the polarizers $`\alpha `$ and $`\beta `$ (0, 45, 22.5, 67.5). These positions yield the maximum violation of Bell’s inequalities in the standard framework. We performed the simulations for varying threshold values from 0.0 (no threshold) to 0.2 (nearly half the photons remain undetected). For every setting we performed 10 separate runs, each with 10000 pairs of photons, the efficiency of the detectors was assumed to be 100%. Fig. 5 gives the result of our simulation. The violation (computed according to the version of Clauser et al. (CHSH) ) increases with increasing threshold. Furthermore, the limit of violation is close to 2.0 (CHSH value of 3.90) in the final setting. The actual threshold in the experiments can be estimated from the visibility of the correlation function. For decoherence rates of 10 - 20 %, as in the experiments, we can only obtain agreement with reported values (97% visibility ), if we increase the threshold $`\mathrm{\Delta }s`$ to more than 0.10. Our simulations indicate that the violation depends crucially on the threshold. The maximum violation can reach a value of as much as 2.0 (the limiting case). It seems thus safe to conclude that the Bell inequalities are no suitable criterium to decide on the nonlocality issue.

Summary.\- In summary we have presented a numerical simulation of EPR experiments under the assumption of strict locality and analyzed the effects of polarizer thresholds, limited efficiency, and decoherence in detail. We could show that Bell’s inequalities are violated in these simulation, and that the violation depends crucially on the threshold of the polarizer switches. The limit of violation in this model is about 2.0 (CHSH value of 4.0). We found that existing measurements can be fully accounted for by limited coherence and non-ideal polarizer switches.

Acknowledgements. \- Discussions about the EPR problem over the last year proved very helpful. I am especially indebted to A. Fisher, J. Gavartin, R. Stadler, and B. Wolkow. Computing facilities at the HiPerSPACE center were funded by the Higher Education Funding Council for England.


text:
# Classical programmability is enough for quantum circuits universality in approximate sense

## 1 Introduction

In the paper  was discussed conception of programmable quantum gate arrays, i.e., some quantum circuits are acting on a system in form $`|d;P|d|P`$ considered as data register $`|d`$ and program register (or simply program) $`|P`$. Similar with conception of usual classical computer, it was considered, that circuit acts as some fixed unitary transformation $`𝖴`$ on whole system and different transformations $`U_P`$ of data related only with content $`P`$ of program register, i.e:

$$𝖴\left(|d|P\right)=(U_P|d)|P^{}.$$
(1)
It should be emphasized, here $`U_P`$ is same for any state of data register $`|d`$ i.e., depends only on program $`|P`$, and states of these two registers are not entangled before and after application of $`𝖴`$.

In relation with such a definition in  was noticed, that if to write the Eq. (1) for two different programs $`|P`$ and $`|Q`$ and corresponding unitary data transformations $`|U_P`$ and $`|U_Q`$, it is simple to find by considering scalar products of both parts, that if $`U_P\varphi U_Q`$ for some complex number $`\varphi `$ then the states of program register must be orthogonal for different programs:

$$PQ=0.$$
(2)
The property Eq. (2) treated in  as demonstration that “… no universal quantum gate array (of finite extent) can be realized. More specifically, we show that every implementable unitary operation requires an extra Hilbert space dimension in the program register.” Due to such a problem in the article  was suggested idea of (exactly) universal stochastic programmable gate array developed further by other authors .

It was no-go theorem for “universality in exact sense” , but fortunately already in initial definition of universal quantum computer  was shown possibility of using some discrete everywhere dense subset in the whole continuous space of unitary operators, it is called sometime “universality in approximate sense” .

Here is shown, that in the approximate sense such universal programmable quantum gate arrays with possibility of approximation of any unitary transformation of data register with given precision are really exist and constructions are quite simple and may be described directly. But the property Eq. (2) is again very important, because if all different possible states of program register are orthogonal, it is possible without lost of generality to find implementation of same array with all different programs are elements of computational basis and superpositions of the states for the program register are never used and so any such array can be designed with possibility of using only classical programs. It is shown and discussed with more extent in relation with particular design used below.

## 2 Universal programmable quantum gate arrays

It is more convenient here to use notation with program register first, i.e., $`|P;d`$ instead of $`|d;P`$. The reason is more simple form of operations like CONTROL-$`U`$:

$$\left(\begin{array}{cccc}\hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill u_{00}& \hfill u_{01}\\ \hfill 0& \hfill 0& \hfill u_{10}& \hfill u_{11}\end{array}\right).$$
(3)
Such operator acts on second qubit as $`U`$ only if first qubit is $`|1`$, but $`|0`$ does not change anything. A straightforward generalization for arbitrary $`N\times N`$ matrix $`U`$ and one program qubit is $`2N\times 2N`$ matrix:

$$\left(\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill U\end{array}\right),$$
(4)
where 1 and 0 are $`N\times N`$ unit and zero matrices.

Let us now consider simple programmable gate array with dimension of program is $`M`$ and data register is $`N`$. For quantum computation with qubits $`M=2^m`$ and $`N=2^n`$. Let us introduce a matrix (cf. ) $`𝖲=𝖲_{\{U_1,\mathrm{},U_M\}}`$:

$$𝖲=\left(\begin{array}{cccc}\hfill U_1& & & \hfill 0\\ & \hfill U_2& & \\ & & \hfill \mathrm{}& \\ \hfill 0& & & \hfill U_M\end{array}\right).$$
(5)
where $`U_1,\mathrm{},U_M`$ are $`M`$ matrices $`N\times N`$. It is clear, that such $`NM\times NM`$ matrix $`𝖲=𝖲_{\{U_1,\mathrm{},U_M\}}`$ acts as:

$$𝖲\left(|P|d\right)=|P(U_P|d).$$
(6)
and it corresponds to Eq. (1) with $`P=P^{}`$ is simply number of matrix $`U_P`$.

Such matrix was suggested for description of conditional quantum dynamics in , but Eq. (2) introduced in  and discussed here shows also, that we do not need more general control with superpositions of basis states of program (control register).

Let us consider a case with $`\{U_1,\mathrm{},U_M\}`$ are universal set of gates (in approximate sense) for some quantum circuits with $`n`$ qubits together with unit matrix $`U_0\mathrm{𝟏}`$. Here $`N=2^n`$ and for design with matrix $`𝖲`$ size of program is small enough — number of qubits in the register is $`m=\mathrm{log}_2M`$. But if due to technical problem with implementation we should use set of CONTROL-$`U_i`$ gates Eq. (4) instead of $`𝖲`$ Eq. (5), then number of program qubits are $`m=M`$.

We may use small universal sets with only one and two-qubit gates . For example, if to choose a set with $`n+2`$ one-qubit gates and $`n1`$ two-gates acting on pairs of neighboring qubits from  together with unit matrix, then number of qubit in control register is $`m=\mathrm{log}_2(2n+2)`$ in best case and even in worst case discussed earlier it is $`2n+2`$.

Let us now extend program register up to $`km`$ qubits, i.e. state of arrays described now as $`|P_k;\mathrm{};P_2;P_1;d`$ and together with matrix $`𝖲`$ acting only on last $`m+n`$ qubits $`|P_1;d`$ as described above let us consider unitary matrix $`𝖱`$ acting on program register as right cyclic shift:

$$𝖱|P_k;\mathrm{};P_2;P_1;d=|P_1;P_k;\mathrm{};P_2;d$$
(7)
It is clear, that if $`𝖴=\mathrm{𝖱𝖲}`$, then operator $`𝖴^k`$ can perform with data arbitrary sequence of up to $`k`$ operators $`\{U_1,\mathrm{},U_M\}`$:

$$𝖴^k\left(|0;\mathrm{};0;P_l;\mathrm{};P_1|d\right)=|0;\mathrm{};0;P_l;\mathrm{};P_1(U_{P_l}\mathrm{}U_{P_1}|d)$$
(8)
here $`U_0`$ is unit operator used to fill out $`kl`$ positions for $`l<k`$ and $`U_{P_l}\mathrm{}U_{P_1}U`$ is arbitrary product of gates from the finite universal set $`\{U_1,\mathrm{},U_M\}`$ for approximation of any data gate $`U`$ with necessary precision. Usually number of terms in product $`lM`$. If minimal error of simulation is given, it is possible to find necessary number $`k`$ and then $`𝖴^k`$ is universal circuit with bounded chosen accuracy.

The idle $`kl`$ steps make more clear advantage of a design well known in classical theory of computations — instead of consideration of whole circuit $`𝖴^k`$ it is possible to introduce notion of one computational step $`𝖴`$. In such approach action of quantum circuit similar with usual CPU “timing” with repeating of same $`𝖴`$ till halting.

It is also useful, then due to algorithm of approximation we need to apply same loop many times, for example instead of some operator $`U`$ may be used “$`j`$-th root” $`U^{1/j}`$  and the circuit should be applied $`j`$ times and can be expressed as $`𝖴^{kj}`$.

Here were examples of some general architecture for programmable quantum gate arrays — we have $`n`$ data qubits and program with $`p+m`$ qubits of two different kinds. Here are $`m`$ control qubits of quantum controller together with $`n`$ data qubits form input for “quantum step operator” $`𝖲`$ described by Eq. (5).

Generalization of $`𝖱`$ operator is not necessary simple rotation, because such model may use too many space. The $`𝖱`$ is equivalent of reversible classical111“Qu-ERCC” — Quantum Equivalent of Reversible Classical Circuit, or simply “classical gate” in terminology used by  circuit acting on $`p+m`$ (qu)bit program register with purpose to generate necessary index of universal quantum gate for quantum controller during simulation.

It should be mentioned also yet another advantage and special property of “pseudo-classical”, i.e., orthogonal states of program register. Only for such kind of programs it is possible to make measurements of state without destruction and so use tools like halt bit. It is also further justification of notion of variable lenght algorithms already used above.

The consideration shows possible model of programmable “quantum chip”. It has three different kinds of “wiring”: quantum, intermediate and classical “buses” with $`d`$, $`m`$ and $`p`$ (qu)bits respectively and two different kinds of circuits: quantum controller acting on quantum bus and controlled by intermediate bus and reversible classical (or Qu-ERCC) circuit with input via classical bus for control of intermediate bus.

Due to result about orthogonality of programs  the design is general for deterministic circuits with pure states, and does not need some special quantum control, i.e, superpositions in program register, only data is necessary quantum.

On the other hand, for circuits with mixed states, stochastic programmable quantum arrays  control already is not limited by the (pseudo) classical case discussed here.


text:
# Magnetohydrodynamic turbulence in warped accretion discs

## I Introduction

Warped accretion discs have been a part of the astronomical vocabulary since the discovery of the 35-day cycle of the X-ray binary Her X-1 tananbaum ; katz ; roberts . Incidentally Her X-1 was the first occulting X-ray pulsar for which an optical counterpart was found as discussed by Neta Bahcall at the 6th Texas Symposium in 1972 bahcall . Later on warped accretion discs have been found in a multitude of systems. In later years one of the most interesting examples has been the maser source in the active galactic nucleus of NGC 4258 miyoshi .

While the warped accretion disc offered a simple interpretation of the observations, it was not easy to describe it theoretically at the hydrodynamic level. The problems have been both to explain the excitation mechanism of the warp and its coherence. Pringle pringle  showed that the radiation pressure from the central radiation source may produce a warp in the outer disc, and Schandl & Meyer SM  described a similar mechanism in which the irradiation produces a wind, which in its turn excites a warp.

A crucial condition for any of these mechanisms to work is that the tendency of the disc to straighten itself must be sufficiently weak. There are two different forces that strive to produce a flat disc. Firstly there is the usual viscous stress due to the local turbulence, which is also driving the accretion, but in general it is insignificant compared to the hydrodynamic stress due to the epicyclic shear flow which is driven by the warp itself PP ; PL . The amplitude of the epicyclic motion is inversely proportional to the ordinary turbulent viscosity. For that reason the hydrodynamic stress due to the epicyclic motion will also be inversely proportional to the turbulent viscosity, and the time scale for flattening the disc will be anomalously short compared to the ordinary viscous time scale.

The warping instability therefore requires a mechanism that can damp the epicyclic motion much more efficiently than it transports the angular momentum in the radial direction. That would for instance be the case if the turbulent viscosity was strongly anisotropic. The intention of this paper is to estimate how anisotropic the turbulent viscosity is and to check whether there are any other mechanisms that can limit the amplitude of the epicyclic motion. We describe the numerical model that we use for these estimates in Sect. 2. There are then two ways in which we have studied the interaction between the turbulence and the epicyclic motion. Firstly we have studied the free decay of an epicyclic motion (Sect. 3) and secondly we have studied the motion that results from a radial forcing (Sect. 4). We discuss and summarize our results in Sect. 5.

## II The mathematical model

We solve the magnetohydrodynamical equations in a Keplerian shearing box BNST ; HGB ; MS . Our units are chosen such that $`H=GM=\mu _0=1`$, where $`G`$ is the gravitational constant, $`M`$ the mass of the central object, and $`H`$ is the Gaussian scale height of the shearing box. The density distribution assuming isothermality can then be written as $`\rho =\rho _0e^{z^2/H^2}`$, where we put $`\rho _0=1`$. The physical size of the shearing box is $`L_x:L_y:L_z=1:2\pi :4`$, and the box is positioned such that $`x`$ and $`z`$ vary between $`\pm \frac{1}{2}L_x`$, and $`\pm \frac{1}{2}L_z`$, respectively, while $`y`$ goes from 0 to $`L_y`$, and the distance of the origin to the central object, $`R_0=10`$. This gives the orbital period $`T_0=199`$, and the mean internal energy $`e_0=\mathrm{7.4\hspace{0.17em}10}^4`$. To stop the box from heating up we add a cooling function

$$Q=\sigma _{\mathrm{cool}}\left(ee_0\right),$$
(1)
where $`\sigma _{\mathrm{cool}}`$ is the cooling rate, which typically corresponds to a time scale of 1.5 orbital periods. Our boundary conditions are (sliding) periodic in the ($`x`$-) $`y`$-direction. In the $`z`$-direction they are impenetrable and stress-free for the velocity, and acts as a perfect conductor with respect to the magnetic field.

We start the simulations from a snapshot from a previous simulation in which the magnetohydrodynamic turbulence is already fully developed. In the first set of simulations we then add a radial velocity of the form $`u_x=u_0\mathrm{sin}(\pi z/L_z)`$, which will have the time evolution of an epicyclic motion. In the second set of simulations, we do not modify the velocity field of the initial snapshot, but rather add a radial forcing term to the equation of motion. This forcing gives an acceleration with a harmonic time-dependence on the orbital time scale and the same $`z`$-dependence as the velocity above.

## III The free decay of an epicyclic motion

In the first set of simulations we add an epicyclic motion to the initial state of the simulation, and then follow the decay of the epicyclic motion. These simulations have previously been described in torkel . With a weak epicyclic motion, its maximum Mach number is initially 0.38, it is difficult to follow the evolution of the epicyclic motion as it is comparable in size to the turbulent velocities. When the amplitude is increased to a Mach number of 3.3, we can distinguish two stages in the damping. After a brief period of essentially no damping the epicyclic velocity quickly drops by a factor of 2 to 3. This damping is followed by an extended phase of exponential decay with a time scale of 25 orbital periods. The time scale of the exponential decay can be translated to a Shakura-Sunyaev SS  $`\alpha `$-parameter of 0.006, which is within a factor of two of the values usually derived from turbulence simulations, e.g. BNST ; stone . The preceding rapid damping is a new phenomenon though, which we interpret in terms of that the epicyclic motion is decaying to inertial waves via a parametric instability GGO . We also note that the toroidal magnetic field in the shearing box reverses its direction during the damping of the strong epicyclic motion (Fig. 1).

## IV Driven epicyclic motion

The dynamics in a real accretion disc is significantly different from the case we have studied above. In reality the epicyclic motion will be driven by the radial pressure gradient that is set up by the warp. To mimic this we have carried out a new set of simulations in which we drive the epicyclic motion by adding a time-periodic radial force to the equation of motion.

We plot the horizontally averaged velocities $`u_x`$ and $`u_y`$ in Fig. 2. The epicyclic motion is in this case comparable to that of Run 1b in torkel . The results are similar in the two cases, and the epicyclic motion is difficult to distinguish at $`z=1.17`$ due to the effect of the turbulent stresses, while it is easily distinguishable at $`z=0.29`$.

## V Conclusions

In this paper we have studied the dynamics of an epicyclic flow in a Keplerian shear flow. In our simulations we find two damping mechanisms for the epicyclic motion. The turbulent stresses can damp the motion in a way which can be described in terms of a turbulent viscosity comparable in strength to that driving the radial angular momentum transport, but a sufficiently fast epicyclic motion can lose significant amounts of energy by exciting inertial waves through a parametric instability.

## Acknowledgments

Computer resources from the National Supercomputer Centre at Linköping University are gratefully acknowledged. UT is supported by the Swedish Research Council (formerly the Natural Sciences Research Council, NFR), and RFS is supported by NASA grant NAG5-4031. This work was supported in part by the Danish National Research Foundation through its establishment of the Theoretical Astrophysics Center (ÅN).


text:
# Variability in the extreme helium star LSS 5121Based on observations made at the University of Texas McDonald Observatory, Fort Davis, Texas and the Jacobus Kapteyn Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

## 1 Introduction

Extreme helium stars (EHes) are luminous blue stars with spectra that indicate surfaces composed almost entirely of helium, with traces of carbon and nitrogen and very little (usually $`1`$ per cent) hydrogen. Their temperatures fall in the range $`8000\mathrm{K}\mathrm{T}\mathrm{32\hspace{0.17em}000}\mathrm{K}`$. The EHes are of interest because, among other things, some of them pulsate or are otherwise variable, they have very unusual chemical compositions, and they are quickly evolving, with the radii and temperature changes of some being measurable over a period of a decade or less (Jeffery et al. jshp01 (2001)). Although their evolutionary history is still not completely certain, it is evident that they are post AGB stars evolving to become white dwarfs. EHes evolution may be linked with other hydrogen deficient star classes such as RCrB stars and helium subdwarf B stars (He sdBs). It may be that a sizable fraction of stars go through an EHe phase.

Pulsations in EHes are excited through iron group (Z-bump) opacities (Saio s93 (1993)). The pulsations in some EHes allow accurate estimates to be made of their physical properties, including their mass and radii (Lynas-Gray et al. lsh84 (1984), lks87 (1987), Woolf & Jeffery wj00 (2000)). The study of pulsations in stars with atmospheres not dominated by hydrogen (e.g. EHes and He sdBs) provides additional constraints on pulsation theory. Fewer than 30 EHes are known (Jeffery et al. jhh96 (1996)) and few of those have been studied in any detail. The physical properties and the variability of additional EHes must be determined if a statistical study of the group is to have meaning.

LSS 5121 is a hot extreme helium star. The temperature of LSS 5121 has been estimated as 28 300 K (Heber et al. hjd86 (1986)) using UBV color measurements and 29 772 K (Jeffery et al. jshp01 (2001)) by comparing optical photometry and IUE ultraviolet spectra with the calculated theoretical flux distributions. Earlier spectroscopic observations by Lawson and Kilkenny (lk98 (1998)) demonstrated that its absorption lines vary in strength on a timescale of a few days. The line variations are similar to those observed in the ($`\mathrm{T}20\mathrm{kK}`$) EHe BD$-9°4395$, a non-radial pulsator (Jeffery & Heber jh92 (1992)). The spectrum of LSS 5121 is very similar to that of the ($`\mathrm{T}29\mathrm{kK}`$) EHe HD 160641 (Heber et al. hjd86 (1986)), which has been found to pulsate non-radially with several pulsation periods which range between 8.4 hours and 1.77 days (Lynas-Gray et al. lks87 (1987)). The spectral similarity indicates similar temperatures, gravities, and chemical compositions. The gravity of LSS 5121 has not been determined, but it is clear from its spectrum that its lies between those of the EHes HD 160641 and LS IV$+6°2$. As BD$`9^{}`$4395, LSS 5121 and HD 160641 all lie near or above the Z-bump instability boundary (Saio & Jeffery sj99 (1999)) it would be reasonable to suppose that LSS 5121 may also pulsate. Given the similarity of the physical properties of LSS 5121 and HD 160641, one might also expect LSS 5121 to pulsate non-radially. EHes with physical properties like those of LSS 5121 tend to vary on timescales of hours to days, rather than minutes or months. Our observations were designed to detect variability on the expected timescales.

LSS 5121 was too dim, $`m_\mathrm{V}=13.3`$, for Lawson and Kilkenny (lk98 (1998)) to test for photometric variations with the 0.5 m telescope used in their observations. Although preliminary observations of LSS 5121 at the Jacobus Kapteyn Telescope (JKT) in 1999 June showed probable variability, LSS 5121 was not the primary target for the observing run, so the observations did not cover enough time for the period search to be conclusive (Aznar Cuadrado et al. arj2000 (2000)). The observations we now report were carried out in 2000 May and June to provide better measurements of the star’s variability.

## 2 Observations and data reduction

### 2.1 JKT photometry

Photometric measurements of LSS 5121 were performed using the JAG-CCD camera on the 1-m JKT at the Roque de los Muchachos Observatory, La Palma on the nights of 2000 May 13 and 14. Observations were made using the Cousins V filter. Pixels on the $`2048\times 2048`$ chip covered 0.33 arcsec per pixel, giving a usable field of about 10 arcmin. Each night the observations covered about 4.2 hours, with exposures being taken every 3 or 4 minutes.

### 2.2 McDonald photometry

LSS 5121 was observed on the nights of 2000 June 10, 11, 13, 14, 15, and 21 using the Prime Focus Corrector of the McDonald Observatory 0.8-m telescope. Observations were made using the Bessel V filter. The camera gives 1.3 arcsec per pixel using the $`2048\times 2048`$ pixel CCD, giving a 46 arcminute field. As LSS 5121 is in a fairly crowded field, we were able to reduce the portion of the chip read out, thus reducing the readout time and disk space requirements, while still including enough comparison stars in the observation. The photometric data from both McDonald and JKT were obtained for use as differential photometry. The times covered by observations at JKT and McDonald are listed in Table 1.

### 2.3 McDonald spectroscopy

Spectra of LSS 5121 were obtained on the nights of 2000 June 13, 14, and 15 using the Sandiford Echelle spectrograph (McCarthy et al. mccarthy (1993)) on the 2.1-m Otto Struve Telescope. Repeated 15 to 20 minute exposures were made during the times listed in Table 1 to allow us to detect changes in the spectrum due to pulsation. Spectra typically had a signal to noise ratio between 7 and 11. The wavelength range of the echelle orders overlapped, so the spectra give continuous coverage between 4410 and 4960 Å. We used a 1.65 arc sec slit width, which gives a spectral resolution of $`\lambda /\mathrm{\Delta }\lambda \mathrm{40\hspace{0.17em}000}`$.

### 2.4 Data reduction

Photometric data were reduced using standard iraf packages to subtract bias and divide by flats. iraf’s phot package was used to do aperture photometry on the images. Stars in uncrowded parts of the field and falling on parts of the CCD without bad pixels or columns were chosen for use as comparison stars for the differential photometry by visually inspecting the images.

Spectra were reduced using standard iraf packages for bias and flat field correction, reducing echelle orders to one dimensional spectra, and applying the wavelength scale using thorium-argon arc spectra. Velocity corrections for Earth’s motion were found for each exposure using rvcorrect and were applied using dopcor.

## 3 Analysis

### 3.1 Differential photometry

Comparison stars were chosen for each night’s data. To be chosen, a star needed to be present on a good part of the chip for every LSS 5121 image for the night. The flux measured for LSS 5121 and the combined flux measured for the comparison stars were used to find $`\mathrm{\Delta }\mathrm{𝑚𝑎𝑔}=m_\mathrm{V}(\mathrm{LSS}5121)m_\mathrm{V}(\mathrm{comparison})`$. Corrections were then made to put all the data on the same scale, taking into account that different reference stars were used on separate nights. Differential photometry for the two observing runs are shown in Fig. 1.

### 3.2 Radial velocities

Velocity shifts were measured using the cross correlation package fxcor in iraf. The spectra were coadded to provide a template with a signal to noise ratio of about 35. Only the two echelle orders covering 4460 to 4500 Å and 4880 to 4930 Å had features, notably the 4471.5, 4920.8 and 4921.9 Å He i lines, strong enough for the cross correlation to yield a velocity. We gave more weight to the velocities found using the longer wavelength order, as the order containing the 4471.5 Å line had lower signal. Other orders contained lines that may be useful for velocity determination in studies using larger telescopes, but which did not give reliable cross correlation peaks with our spectra. The velocities thus found were applied and the resulting spectra were coadded to give a template with smaller velocity smearing for a second iteration. Velocities found for three nights are displayed in Fig. 2.

We estimate uncertainties in radial velocity determinations to be about $`\pm 10\mathrm{km}\mathrm{s}^1`$. This is a rough estimate based on how accurately we could determine the center of the cross correlation peak.

## 4 Results

The uncertainties for the radial velocity data are large compared to differences between the calculated velocities: we cannot say that radial velocity variations are obvious. The differential photometry shows evidence that LSS 5121 varies on a timescale of a few hours, however no period was immediately apparent. The photometry measured at JKT shows a magnitude change of about 0.04 mag over a period of about 2.5 hours. The photometry measured the next month at McDonald shows similar variation, though the data are noisier. Uncertainties for the photometric data varied with the weather, but in general we estimate uncertainties of $`\pm 0.003`$ mag for the JKT data and $`\pm 0.015`$ mag for the McD data. Noise was reduced by finding the mean for the photometric data in 0.02 day bins (Figs. 3 & 4). Tests with non-variable comparison stars with magnitudes similar to LSS 5121 gave differential photometry which did not vary more than about 0.01 mag.

Periodograms and window functions were found for the photometric and velocity data using dipso’s pdgram and pdgwin features (Scargle sc82 (1982)). The velocity periodogram has no major peaks, which is not surprising, given the small number of data points. The photometry periodogram has no major peaks that are not also in the data’s window function (Fig. 5). The major peaks in both are near multiples of the one day alias.

Data were phased to the periods indicated by minor peaks in the periodogram not present in the window function, but this provided no further evidence of periodic behavior. We cannot rule out that LSS 5121 varies with a period of about 1 day or that it varies with longer periods and/or smaller amplitudes than our observations could detect.

## 5 Discussion and Conclusions

Differential photometry shows that LSS 5121 varies in brightness, with changes of up to 0.04 mag over a span of 2.5 hours. LSS 5121 was previously known to have varying spectral lines. We have now shown that it varies in brightness as well. Analysis of data from observations on 8 nights in 2000 May and June did not yield evidence of any periodicity in the variations. Our search for variability in radial velocity was inconclusive. A successful radial velocity search will require a larger telescope to provide higher signal to noise spectra.

Saio (s95 (1995)) showed that the region of the HR diagram where pulsations should occur overlaps part of the region where EHes appear. Figure 6 shows the location of all known EHes with known $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$, and the Fe-group opacity instability boundary for 0.7 M helium stars in the $`\mathrm{log}g\mathrm{log}T_{\mathrm{eff}}`$ plane.

The EHes mentioned in the text are labeled in the figure. Helium stars above the boundary should be unstable against pulsation. The exact location of the instability boundary depends on stellar mass and chemical composition.

In general, the observed pulsation timescales in EHes decrease with higher temperature (Saio & Jeffery sj88 (1988)). Two EHes, LSS 3184 (BX Cir) and V652 Her pulsate in their fundamental radial modes with well determined periods of about 2.6 hours (Kilkenny et al. klr96 (1996), kkj99 (1999)). Pulsations in the coolest EHes (e.g. BD$+1°4381$) seem to be non-periodic or chaotic with variation timescales of 8 to 25 days (Lawson et al. lkw93 (1993)). Other EHes with known variability seem to pulsate non-radially, some with multiple periods. For example, Lynas-Gray et al. (lks87 (1987)) report finding periods of 0.35, 0.71, 1.12, and 1.77 days in the photometry of HD 160641, an EHe with a temperature similar to that of LSS 5121 and Jeffery et al. (jsh85 (1985)) report 3.5 day, 11.2 day, and possibly additional periods for BD$-9°4395$. It is possible that chaotic pulsations also affect some of the hotter EHes, so that some of the multiple periods found for the hot non-radial pulsators would more accurately be referred to as variability timescales or “quasi-periods,” and that different periods or no periods would be found if data from longer, more continuous observations were available.

Any pulsations of LSS 5121 are either non-periodic or have periods which our data cannot resolve. Our inability to find a period is consistent with LSS 5121 being a non-radial or a chaotic pulsator, though it does not constitute proof. Its pulsations are probably similar to those of the hot extreme helium stars BD$-9°4395$ and HD 160641, both of which have been reported to pulsate non-radially with multiple periods (Jeffery & Heber jh92 (1992), Lynas-Gray et al. lks87 (1987)).

###### Acknowledgements.

We acknowledge financial support from the UK PPARC (grant Refs PPA/G/S/1998/00019 and PPA/G/O/1999/00058). VW thanks Marcel Bergmann for instruction on use of the 0.8-m telescope and Nairn Baliber for a refresher on photometry reduction. GP received support from grant AST 9618414 from the NSF.