text:
# A Stochastic Model of Human Gait Dynamics

\[

## Abstract

We present a stochastic model of gait rhythm dynamics, based on transitions between different “neural centers”, that reproduces distinctive statistical properties of normal human walking. By tuning one model parameter, the hopping range, the model can describe alterations in gait dynamics from childhood to adulthood — including a decrease in the correlation and volatility exponents with maturation. The model also generates time series with multifractal spectra whose broadness depends only on this parameter.

\]

Many physical and physiological processes exhibit complex fluctuations characterized by scaling laws , some with monofractal  and others with multifractal behavior . For example, normal gait dynamics of adults are monofractal , while healthy heartbeat fluctuations are multifractal . The origin of these complex fluctuations and the factors that contribute to differences in their behavior are largely unknown. Here we develop a physiologically-motivated model that may be helpful in explaining some of the properties that contribute to complex dynamics. We focus on one class of signals — the time series of the inter-stride-intervals (ISI) between successive strides in human gait.

Walking is a voluntary process, but under normal circumstances, stride to stride regulation of gait is controlled by the nervous system in a largely automatic fashion . Gait is regulated in part by “neural centers” within the cortex and the spinal cord . As a result of time-varying inputs, neural activity at different times are thought to be dominated by different centers, leading to complex fluctuations in the ISI “output.” Furthermore, as the central nervous system matures from infancy to adulthood, the interaction between neural centers becomes richer .

To understand the underlying regulatory mechanisms of walking, deterministic and stochastic models have been proposed. For example, classic “central pattern generator” models are based on oscillatory neural activity, where the interaction between neural centers helps regulate gait dynamics . A stochastic version of a central pattern generator model reproduces certain fractal properties of the ISI series . However, existing models do not explain observed changes in scaling exponents , and volatility (magnitude) correlations  (a new finding reported in the present study) that occur during gait maturation from childhood to adulthood.

We propose a stochastic model consisting of a random walk (RW) on a chain, the elements of which represent excitable neural centers . A step of the RW between element $`i`$ and element $`j`$ represents the “hopping” of the excitation from center $`i`$ to center $`j`$. The increase of neural interconnectedness with maturation is modeled by increasing the range of “jump” sizes of the RW, since larger jump sizes will allow exploration of more neural centers. This property mimics one aspect of the increasing complexity of the adult nervous system.

Previous studies  have identified neural centers with pacemaker-like qualities that fire with frequency $`f_i`$, so we represent the network of neural centers by different frequency modes. One mode is activated at a given time ($`ISI1/f_i`$), and the $`f_i`$ are Gaussian distributed. The model is based on the following assumptions (Fig. 1) :

$``$ Assumption (i) is that the $`f_i`$ have finite-size correlations, $`f_if_{i+\delta }/f_i^2=e^{\delta /\delta _0}`$. We assume finite-range correlations among $`f_i`$ because neighboring neurons are likely to be influenced by similar factors . This assumption effectively creates “neuronal zones” composed of neural centers (modes) along the chain with a typical size $`\delta _0`$.

$``$ Assumption (ii) concerns the rule followed by the RW process. The active neural center is determined by the location of the RW. The “jump” sizes of the RW follow a Gaussian distribution of width $`C`$.

$``$ Assumption (iii) is that a small fraction of noise is added to the output of each mode to mimic biological noise not otherwise modeled. The output $`y`$ becomes $`y(1+A\eta )`$ where $`A`$ is the noise level and $`\eta `$ is Gaussian white noise with zero mean and unit variance .

The model has three parameters $`\delta _0`$, $`C`$, and $`A`$. We find that the best agreement with the data is achieved when $`A=0.02`$ and $`\delta _0=25`$. In order to simulate changes with maturation, we vary only the third parameter, $`C`$, as a function of age, $`C=`$ (age – 2) for ages 3 to 25 years. Increasing the hopping range with age is consistent with the fact that neural transmission is not fully developed until the late teens .

Examples of ISI time series are shown in Fig. 2a. The ISI series of the adult subject has smaller fluctuations compared with the more variable ISI of the child. Two examples of the model’s output are shown in Fig. 2b; using a large value of $`C`$ ($`C=25`$) simulates the ISI series of an adult, while using a small value ($`C=3`$) mimics the ISI series of a child. Surprisingly, we find that the normalized probability distributions of the child ISI increment series, $`\mathrm{\Delta }(ISI)_i=(ISI)_{i+1}(ISI)_i`$, and of the model with small $`C`$ are close to exponential (Fig. 3a). In contrast, the distribution converges to Gaussian in the adult and for large $`C`$ in the model (Fig. 3b).

To further test the model, we study time correlations in the ISI series. We calculate the function $`F(n)`$, which corresponds to the rms fluctuations of the integrated ISI series of the adult (Fig. 2a) and of the model for $`C=25`$ (Fig. 2b) . Estimating the scaling exponent from the slope of $`\mathrm{log}F(n)`$ vs $`\mathrm{log}n`$, we find, for both data and model, that the ISI time series have correlations with scaling exponents $`0.8`$ (Fig. 4a).

The scaling exponents of the original signal provide an indication of the linear properties (two-point correlations). Certain nonlinear aspects may be associated with the presence of long-range correlations in the magnitudes of ISI increments $`|\mathrm{\Delta }(ISI)_i|`$, an index of “volatility” . We find (Fig. 4b) that the adult magnitude series is uncorrelated (the scaling exponent is close to 0.5). The model shows similar behavior for large $`C`$ (Fig. 4b).

Next we inquire if changes in gait dynamics from childhood to adulthood might be reflected in changes in the scaling exponents of the ISI series. The short-range scaling exponents of the ISI series decrease as children mature . We compare the short range scaling exponents for a group of 50 children  with those of 10 adults ; this exponent decreases from $``$ 1.0 to $``$ 0.7 (Fig. 5a). By altering $`C`$, the model simulates these maturation-related changes in two-point correlation properties.

The magnitude series exponent of the ISI series also decreases with maturation (Fig. 5b). Our results for the magnitude series suggest that the gait pattern of children is more volatile (and thus more nonlinear) than the

usual walking pattern of adults. We note that: (i) unlike adults, young children have difficulty keeping their walking speed constant. A large “walking error” in one direction is likely to be followed by a large compensatory walking error in the opposite direction — i.e., a large (or small) increment of ISI is likely to be followed by a large (or small) decrement over a given range of scales. This instability may lead to increased magnitude correlations. (ii) Adults can voluntarily simulate the less stable gait of children and thereby increase the magnitude correlations of the ISI. Children, however, cannot maintain the less volatile dynamics of adults. Thus, in general, the more mature adult gait dynamics are likely to be richer. The model shows (Fig. 5b) a similar decrease of the magnitude exponents when increasing $`C`$ and is within the error bars of the data.

The decrease of the magnitude exponent with maturation is consistent with the possibility that the ISI time series (during goal-directed walking) becomes less multifractal as the individual matures. This prediction cannot be tested on the available data since the time series of the children are too short for multifractal analysis; but this prediction can be tested on the model. Fig. 6 summarizes a systematic multifractal analysis on ISI series generated by the model. We find narrower multifractal spectra with increasing $`C`$, suggesting the possibility of a decrease of multifractality with maturation during normal walking. The model’s prediction agrees with a recent finding regarding the monofractality of normal gait in adults . The model also serves as a generator of random series with different multifractal spectra — multifractal for $`C1`$ and monofractal for $`C1`$.

In summary, we find that a simple stochastic model captures multiple aspects of gait dynamics, and their changes with maturation, including: (i) the shape of the probability distribution of the ISI increments; (ii) correlation properties of the ISI; and (iii) correlations properties in the magnitudes of the ISI increments . Further, by varying only a single “hopping-range” parameter, $`C`$, a wide array of multifractal dynamics can be generated. The model can also be altered by “knocking out” certain frequency modes (akin to what may occur during very advanced age or in response to neurodegenerative disease). Simulation with drop-out of frequency modes predicts increased gait variability, with (i) increased magnitude exponents, and (ii) decrease of long-range correlations. Our preliminary analysis of the ISI series of older adults prone to falls is consistent with this prediction. Generalization of the model to two and three dimensional networks to describe other types of neurological activities is underway.

Partial support was provided by the NIH/National Center for Research Resources (P41 RR13622) and the NIA (AG14100). We thank J.J. Collins, S. Havlin, V. Schulte-Frohlinde, and C.-K Peng for helpful discussions.


format: text

text:
# Verification of the semiclassical method for an electron moving in a homogeneous magnetic field

## I Introduction

In the semiclassical approach, the charge motion is described in a classical way, and quantum transitions are calculated on the base of matrix elements . This method is known to be fairly in the ultrarelativistic case at ultrahigh electron energies. For example, at energies of $``$ 2.5 GeV, which are typical for current storage rings with a characteristic parameter $`H\rho 10^7`$ Oe cm the electron energy levels amount to $`n10^{17}`$. For radiative transitions to lower energy levels, $`\mathrm{\Delta }n/n10^7`$. In this case, according to the uncertainty principle, all quantum processes occur in a region as small as a few angstroms (see  for details). Clearly, this is a negligible small value for an orbital radius of about 10 m. Hence, the classical, or more precisely, semiclassical concept of the electron trajectory can be used here.

In the theory of synchrotron radiation, the semiclassical method allows all quantum corrections to the synchrotron radiation power, including the effects with spin-flip to be calculated . This method is straightforward to use. What is more, it provides a visual picture of physical processes occurring at ultrahigh electron energy (recoil effects, radiative spin self-polarization, spin magnetic-moment radiation, mixed charge and magneton radiation, radiation associated with the anomalous part of electron magnetic moment, etc.).

Despite of the impressive success of the semiclassical synchrotron-radiation theory rigorous mathematical substantiation of this method based on exact solutions to the Klein-Gordon equation or to the Dirac-Pauli equation is, to our knowledge, lacking heretofore. The basic statements of the theory were, in fact, postulated on the basis of the uncertainty principle (see ).

In this paper, solutions to the classical equation of charge motion and to the Bargmann-Michel-Telegdi (BMT) spin precession equation  are constructed in the framework of quantum mechanics, using nonstationary wave functions (see general statements in ). The later are, in their turn, exact solutions to the Klein-Gordon equation or to the Dirac-Pauli equation for a charged particle moving in a homogeneous magnetic field. A simplified version of this method was used in  to reveal a physical pattern underlying the behavior of longitudinal electron-spin polarization in a magnetic field. This approach is adopted in  to describe quantum neutron spin-flip transitions.

Our task is to explore the time evolution of the expectation values of the momentum and spin operators for a charged particle moving in a homogeneous magnetic field. To demonstrate the capability of the method discussed, let us examine the problem at hand first. In Sec. II, we will consider a scalar particle and compare the time evolution of the expectation value of the momentum operator for this particle with the solution to the classical equation of charge motion. These are found to differ by a factor. Sec. III will focus on a spin-$`\frac{1}{2}`$ particle. In addition to evolution of the expectation value of the momentum operator, we will discuss the behavior of the expectation value of the spin operator. The latter also differs from the solution to the BMT equation by factor only. Finally, in Sec. IV, we will introduce a procedure that will help to eliminate the difference between the quantum and classical approaches. To this end, the state of a particle will be represented by a wave packet so that the time evolution of expectation values of the corresponding operators with respect to this nonstationary wave function will coincide with that predicted classical mechanics.

## II Motion of a scalar particle

Let us assume that state of scalar particle is formed as a wave packet by superposing exact solutions to the Klein-Gordon equation corresponding to three neighboring energy levels $`n1`$, $`n`$, and $`n+1`$:

$$\mathrm{\Psi }(𝐫,t)=\underset{m=n1}{\overset{n+1}{}}A_m\psi _m(𝐫)\mathrm{exp}(\frac{i}{\mathrm{}}m_0c^2B_mt),$$

where $`\psi _m(𝐫)`$ are the stationary solutions to the Klein-Gordon equation in a homogeneous magnetic field (see Appendix A) and $`A_m`$ are the expansion coefficients.

To determine the coefficients $`A_m`$ in the semiclassical approximation ($`n1`$), the probability that a particle will be found at each of the levels is assumed to be the same. For $`n1`$ this is natural assumption. Normalizing the wave function $`\mathrm{\Psi }(𝐫,t)`$ to unit probability yields a relation of the form

$$\underset{m=n1}{\overset{n+1}{}}A_m^{}A_m=1.$$

In what follows an explicit form of $`A_m`$ is immaterial but different relations between these coefficients are of great importance. In particular, we need the following equation

$$𝔄_{\mathrm{KG}}=A_n^{}A_{n1}+A_{n+1}^{}A_n=\frac{2}{3}.$$

Now let us calculate the expectation value of the momentum operator $`\widehat{𝐏}`$ in the state $`\mathrm{\Psi }(𝐫,t)`$. In the first step we calculate matrix elements of $`\psi _m^{}|\widehat{𝐏}|\psi _m`$ type with respect to stationary solutions to the Klein-Gordon equation. Then, using the semiclassical condition ($`n1`$) for each matrix element \[see Eqs. (A)\] we write $`\mathrm{\Psi }|\widehat{𝐏}|\mathrm{\Psi }_t`$ in explicit form to give

$$\widehat{P}_x_t=\frac{i}{2}m_0cb_{}\left[𝔄_{\mathrm{KG}}\mathrm{exp}(i\omega t)𝔄_{\mathrm{KG}}^{}\mathrm{exp}(i\omega t)\right]=\frac{2}{3}m_0cb_{}\mathrm{sin}\omega t,$$

$$\widehat{P}_y_t=\frac{1}{2}m_0cb_{}\left[𝔄_{\mathrm{KG}}\mathrm{exp}(i\omega t)𝔄_{\mathrm{KG}}^{}\mathrm{exp}(i\omega t)\right]=\frac{2}{3}m_0cb_{}\mathrm{cos}\omega t,$$

$$\widehat{P}_z_t=m_0cb_z.$$

Here

$$\omega =\frac{m_0c^2}{\mathrm{}}(B_{n+1}B_n)=\frac{e_0H}{m_0c^2B_n}\frac{e_0H}{m_0c^2\gamma }$$

is the frequency which, in essence, coincides with the cyclotron frequency of rotation of a classical particle in a plane perpendicular to the magnetic field vector. In this approximation, $`B_n\gamma `$, where $`\gamma `$ is the Lorentz factor.

This leads us to conclude that the expectation values differ from the corresponding classical solutions to the equation of motion of a charged particle only by the factor ($`2/3`$) in the $`\widehat{P}_x_t`$ and $`\widehat{P}_y_t`$ components. Note that the same result was also obtained in  for an ordinary one-dimensional quantum harmonic oscillator. In the discussion below we will show that this factor can be made as large as unity.

## III Charge motion and spin precession of a Dirac-Pauli particle

When considering the Dirac-Pauli particle the spin polarization of the particle must be taken into account. A superposition of exact solutions to the Dirac-Pauli equation corresponding to three neighboring energy levels $`n1`$, $`n`$, and $`n+1`$ has the form

$$\mathrm{\Psi }(𝐫,t)=\underset{\zeta }{}\underset{m=n1}{\overset{n+1}{}}A_{\zeta m}\psi _{m\zeta }(𝐫)\mathrm{exp}(\frac{i}{\mathrm{}}m_0c^2B_{m\zeta }t),$$
(1)
where $`\psi _m(𝐫)`$ are the stationary solutions to the Dirac-Pauli equation (see Appendix B) and $`A_{\zeta m}`$ are expansion coefficients. Note that in this case, allowance is made for polarization states defined by a quantum number $`\zeta =\pm 1`$. Let us assume that a given longitudinal-polarization state at time $`t=0`$***It is known  that, with regard to the anomalous magnetic moment of an electron, the operator $`(𝝈\widehat{𝐏})`$ is not an integral of motion. for an electron located at level $`m`$ is of the form

$$(𝝈\widehat{𝐏})\underset{\zeta }{}A_{\zeta m}\psi _{m\zeta }=\lambda _mA_{\zeta m}\psi _{m\zeta }.$$
(2)
Having solved Eq. (2) we find the eigenvalue

$$\lambda _m=\epsilon \sqrt{b_{}^2+b_z^2},\epsilon =\pm 1.$$

In the semiclassical approximation, it can be supposed that for all values of $`m`$, the spin coefficients are the same, i.e., $`C_i(m,\zeta )C_i(n,\zeta )`$ and, moreover, $`\lambda _m\lambda _n`$ \[see Eqs. (B8) and (B9)\]. Then we have

$$A_{1m}=\kappa A_{1m},\kappa =\frac{b_z+\epsilon b\sqrt{b_{}^2+b_z^2}}{B_{n\zeta }b_{}}.$$

Let us also require that the wave function $`\mathrm{\Psi }(𝐫,t)`$ be normalized to unit probability. As a result, we get

$$\underset{\zeta }{}\underset{m=n1}{\overset{n+1}{}}A_{\zeta m}^{}A_{\zeta m}=1.$$
(3)
In what follows we need the relations

$$𝔄_1=\underset{m=n1}{\overset{n}{}}(A_{1m}^{}A_{1m+1}+A_{1m}^{}A_{1m+1})=\frac{2}{3},$$
(5)
$`𝔄_2`$ $`=`$ $`{\displaystyle \underset{m=n1}{\overset{n}{}}}A_{1m}^{}A_{1m+1}`$ (6)
$`=`$ $`{\displaystyle \underset{m=n1}{\overset{n}{}}}A_{1m}^{}A_{1m+1}={\displaystyle \frac{2}{3}}{\displaystyle \frac{\kappa }{\kappa +1}},`$ (7)
$$𝔄_3=\underset{m=n1}{\overset{n+1}{}}A_{1m}^{}A_{1m}=\frac{\kappa }{\kappa ^2+1},$$
(8)
$$𝔄_4=\underset{m=n1}{\overset{n+1}{}}(A_{1m}^{}A_{1m}A_{1m}^{}A_{1m})=\frac{\kappa ^21}{\kappa ^2+1}.$$
(9)
To simplify the calculations we suppose that, as an electron makes a transition from one level to another, the projection of the momentum of the electron onto the direction of the magnetic field is conserved. In addition, random deviations of the orbital center of the electron due to radiation are neglected. These restrictions imply that $`b_z^{}=b_z`$ and $`s^{}=s`$ (see Appendix B).

Now, by analogy to a spinless particle, we find the expectation value of the operator $`\widehat{𝐏}`$ in Dirac-Pauli particle state (1) (see the corresponding matrix elements in Appendix B). As a result of straightforward but cumbersome calculations, we derive, in view of these restrictions, the following equations:

$$\widehat{P}_x_t=\frac{i}{2}m_0cb_{}\left[𝔄_1\mathrm{exp}(i\omega t)𝔄_1^{}\mathrm{exp}(i\omega t)\right]=\frac{2}{3}m_0cb_{}\mathrm{sin}\omega t,$$
(11)
$$\widehat{P}_y_t=\frac{1}{2}m_0cb_{}\left[𝔄_1\mathrm{exp}(i\omega t)𝔄_1^{}\mathrm{exp}(i\omega t)\right]=\frac{2}{3}m_0cb_{}\mathrm{cos}\omega t,$$
(12)
$$P_z_t=m_0cb_z.$$
(13)
Here

$$\omega =\frac{m_0c^2}{\mathrm{}}(B_{n+1\zeta }B_{n\zeta })=\frac{e_0H}{m_0c^2B_n\zeta }\frac{e_0H}{m_0c^2\gamma }$$

is the frequency which in the BMT approximation (the charge motion is independent of the spin precession) coincides with the cyclotron frequency of rotation of an electron in the plane perpendicular to the magnetic-field vector. In this approximation, $`B_{n\zeta }\gamma `$.

A similar procedure applies to calculations of the expectation value of the spin operator $`\widehat{S}^\mu `$ (see ) to yield

$$\widehat{S}^\mu =(\frac{1}{m_0c}(𝝈\widehat{𝐏}),\rho _3𝝈+\frac{1}{m_0c}\rho _1\widehat{𝐏}+\frac{\mu _a}{m_0c^2}\rho _3𝐇).$$

Taking into account Eqs. (3), we obtain

$$\widehat{S}_x_t=\frac{2}{3}\zeta _{}(\mathrm{cos}\omega t\mathrm{sin}\mathrm{\Omega }_at+b\mathrm{sin}\omega t\mathrm{cos}\mathrm{\Omega }_at),$$
(15)
$$\widehat{S}_y_t=\frac{2}{3}\zeta _{}(\mathrm{sin}\omega t\mathrm{sin}\mathrm{\Omega }_atb\mathrm{cos}\omega t\mathrm{cos}\mathrm{\Omega }_at),$$
(16)
$$\widehat{S}_z_t=\frac{B_{n\zeta }}{b}\zeta _z+\frac{b_zb_{}}{b}\zeta _{}\mathrm{cos}\mathrm{\Omega }_at,$$
(17)
$$\widehat{S}^0_t=\frac{b_z}{b}\zeta _z+B_{n\zeta }\frac{b_{}}{b}\zeta _{}\mathrm{cos}\mathrm{\Omega }_at.$$
(18)
Here $`\zeta _z^2=1\zeta _{}^2`$, $`\zeta _{}=2\kappa (\kappa ^2+1)^1`$ are constants which, as we will be seen later, are analogous to some constants in the classical spin theory. Thus, as in the classical approach, the total precession is determined not only by the cyclotron frequency $`\omega `$ but also by the anomalous frequency

$`\mathrm{\Omega }_a={\displaystyle \frac{m_0c^2}{\mathrm{}}}(B_{n\zeta }B_{n\zeta ^{}})=`$ $`2{\displaystyle \frac{\mu _aH}{\mathrm{}}}{\displaystyle \frac{b}{\sqrt{b_z^2+b^2}}}`$
$`{\displaystyle \frac{\alpha }{2\pi }}{\displaystyle \frac{e_0H}{m_0c}}\sqrt{1\beta _z^2}.`$
Note that in the BMT approximation, the term comprising the anomalous magnetic moment in the spin operator $`\widehat{S}^\mu `$ is immaterial for result (III).

We see that by analogy with the spinless charged particle, the results obtained in Eqs. (11), (12), (15) and (16) are slightly different from the corresponding solutions to the classical equations of charge motion and spin precession by the factor ($`2/3`$).

## IV Semiclassical correspondence principle

To construct a more precise semiclassical theory of charge motion and spin precession, we will consider a wave packet involving exact solutions to Dirac-Pauli equation (B4) corresponding to closely spaced energy levelsIt can be shown that there is a simpler form of this construction for the motion of a spinless charged particle.. In this case, wave function (1) can be written in the following form:

$$\mathrm{\Psi }(𝐫,t)=\underset{\zeta }{}\underset{m\mathrm{}}{}A_{\zeta m}\psi _{m\zeta }(𝐫)\mathrm{exp}(\frac{i}{\mathrm{}}m_0c^2B_{m\zeta }t),$$
(19)
where $`\mathrm{}`$ is a set of values of the principle quantum number. Let us assume that each value of $`m`$ in Eq. (19) belonging to the set is much higher than unity, i.e., $`m1`$. At the same time, we have

$$nN=m_{\mathrm{max}}m_{\mathrm{min}}1,$$

where $`m_{\mathrm{max}}`$ and $`m_{\mathrm{min}}`$ are maximum and minimum values belonging to $`\mathrm{}`$, $`n`$ is the quantum number from $`\mathrm{}`$ defining some energy level to which a given classical trajectory corresponds. Let us further assume that the other energy levels from the set are symmetric about level $`n`$.

The next natural step in our method is to derive the coefficients $`A_{\zeta m}`$. The central idea of the derivation is the same as in Eq. (3), i.e.,

$$A_{1m}=\kappa A_{1m},\underset{\zeta }{}\underset{m\mathrm{}}{}A_{\zeta m}^{}A_{\zeta m}=1.$$

Then, relations (3) can be easily extended to the case under consideration, namely,

$$𝔄_1=\underset{m=m_{\mathrm{min}}}{\overset{m_{\mathrm{max}}1}{}}(A_{1m}^{}A_{1m+1}+A_{1m}^{}A_{1m+1})=\frac{N1}{N},$$
(21)
$`𝔄_2`$ $`=`$ $`{\displaystyle \underset{m=m_{\mathrm{min}}}{\overset{m_{\mathrm{max}}1}{}}}A_{1m}^{}A_{1m+1}`$ (22)
$`=`$ $`{\displaystyle \underset{m=m_{\mathrm{min}}}{\overset{m_{\mathrm{max}}1}{}}}A_{1m}^{}A_{1m+1}={\displaystyle \frac{N1}{N}}{\displaystyle \frac{\kappa }{\kappa +1}},`$ (23)
$$𝔄_3=\underset{m\mathrm{}}{}A_{1m}^{}A_{1m}=\frac{\kappa }{\kappa ^2+1},$$
(24)
$$𝔄_4=\underset{m\mathrm{}}{}(A_{1m}^{}A_{1m}A_{1m}^{}A_{1m})=\frac{\kappa ^21}{\kappa ^2+1}.$$
(25)
Following the above procedure of calculating the expectation values of the momentum and spin operators and using Eqs. (IV) we obtain

$$\widehat{P}_x_t=\frac{N1}{N}m_0cb_{}\mathrm{sin}\omega t,$$

$$\widehat{P}_y_t=\frac{N1}{N}m_0cb_{}\mathrm{cos}\omega t,$$

$$\widehat{P}_z_t=m_0cb_z,$$

$$\widehat{S}_x_t=\frac{N1}{N}\zeta _{}(\mathrm{cos}\omega t\mathrm{sin}\mathrm{\Omega }_at+b\mathrm{sin}\omega t\mathrm{cos}\mathrm{\Omega }_at),$$

$$\widehat{S}_y_t=\frac{N1}{N}\zeta _{}(\mathrm{sin}\omega t\mathrm{sin}\mathrm{\Omega }_atb\mathrm{cos}\omega t\mathrm{cos}\mathrm{\Omega }_at),$$

$$\widehat{S}_z_t=\frac{B_{n\zeta }}{b}\zeta _z+\frac{b_zb_{}}{b}\zeta _{}\mathrm{cos}\mathrm{\Omega }_at,$$

$$\widehat{S}^0_t=\frac{b_z}{b}\zeta _z+B_{n\zeta }\frac{b_{}}{b}\zeta _{}\mathrm{cos}\mathrm{\Omega }_at.$$

Thus, for $`N1`$, the difference of $`\widehat{P}_x_t`$, $`\widehat{P}_y_t`$, $`\widehat{S}_x_t`$, and $`\widehat{S}_y_t`$ from the results obtained by classical theory is eliminated, and the time evolution of the expectation values $`\widehat{𝐏}_t`$ and $`\widehat{S}^\mu _t`$ are made to coincide with the solutions to the corresponding classical equations .

It is easy to see that the anomalous magnetic moment of an electron has effect not only on the time behavior of spin projection onto the direction of motion (longitudinal polarization)  but also on the behavior of spin in the plane perpendicular to the magnetic-field vector, i.e., on the time evolution of $`\widehat{S}_x_t`$ and $`\widehat{S}_y_t`$.

## V Concluding remarks

It should be noted that $`\widehat{S}^\mu _t`$ derived in Sec. IV satisfy the relations given by the classical spin theory:

$$\widehat{S}^\mu _t\widehat{P}_\mu _t=0,\widehat{S}^\mu _t\widehat{S}_\mu _t=1.$$

Moreover, the procedure under review can be extended to the tensor operator

$$\widehat{\mathrm{\Pi }}^{\mu \nu }=(\widehat{𝚽},\widehat{𝚷}),$$

$$\widehat{𝚽}=\frac{1}{m_0c}\rho _3(𝝈\times \widehat{𝐏})+\frac{\mu _a}{m_0c^2}\rho _1(𝝈\times 𝐇),$$

$$\widehat{𝚷}=𝝈+\frac{1}{m_0c}\rho _2(𝝈\times \widehat{𝐏})+\frac{\mu _a}{m_0c^2}𝐇,$$

which, in the BMT approximation, is related to $`\widehat{S}^\mu `$ by the formula :

$$\widehat{\mathrm{\Pi }}^{\mu \nu }=\frac{1}{m_0c}\epsilon ^{\mu \nu \alpha \beta }\widehat{S}_\alpha \widehat{P}_\beta .$$

Thus, the expectation values $`\widehat{S}^\mu _t`$ and $`\widehat{𝐏}_t`$ as well as $`\widehat{\mathrm{\Pi }}^{\mu \nu }_t`$ calculated following the procedure discussed in Sec. IV obey the classical equation of charge motion and spin precession. The approach introduced in this work substantiates the use of the semiclassical method for gaining an insight into physical phenomena associated with high-energy charged particles. In particular, the semiclassical approach can be the basis for a classical model for certain purely quantum processes. For example, it can be used to establish a relationship between spin-flip transitions and classical spin precession .

###### Acknowledgements.

The authors would like to thank Professor V.G. Bagrov for support and valuable discussion.

## A

Motion of a scalar particle in a constant and homogeneous magnetic field $`𝐇=(0,0,H)`$ is given by the Klein-Gordon equation

$$(\frac{1}{c^2}\widehat{E}^2\widehat{𝐏}^2+m_0^2c^2)\psi (𝐫,t)=0,$$

where $`\widehat{E}=i\mathrm{}/t`$ is the energy operator, $`\widehat{𝐏}=\widehat{𝐩}(e/c)𝐀`$ is the kinetic momentum operator, $`𝐀=(Hy/2,Hx/2,0)`$, $`e=e_0<0`$ is the charge, and $`m_0`$ is the rest mass of the particle.

It is known that in this case, the wave function $`\psi (𝐫,t)`$ is to be the eigenfunction of the energy operator $`\widehat{E}`$ as well as $`z`$-components of the kinetic momentum operator $`\widehat{𝐏}`$ and orbital angular momentum operator $`\widehat{𝐋}=𝐫\times \widehat{𝐩}`$ onto the direction of magnetic field vector, i.e.,

$$\widehat{E}\psi =E_n\psi ,\widehat{P}_z\psi =p_z\psi ,\widehat{L}_z\psi =\mathrm{}l\psi .$$

Then, in terms of cylindrical coordinates ($`r`$, $`\phi `$, $`z`$), which are the most suitable ones in this case, the solution to the Klein-Gordon equation has the form

$$\psi (𝐫,t)=\frac{\sqrt{e_0H}}{\sqrt{L}\sqrt{2\pi \mathrm{}c}}e^{iE_nt/\mathrm{}}e^{ip_zz/\mathrm{}}e^{il\phi }I_{n,s}(\rho ).$$
(A1)
The function $`I_{n,s}(\rho )`$ in Eq. (A1) is defined by the Laguerre polynomials $`Q_s^l(\rho )`$ with help of the following relation

$$I_{n,s}(\rho )=\frac{1}{\sqrt{n!s!}}e^{\frac{1}{2}\rho }Q_s^{ns}\rho ^{\frac{ns}{2}},\rho =\frac{e_0H}{2\mathrm{}c}r^2.$$

Here $`n=l+s=0,1,2,\mathrm{}`$ is the principle quantum number, $`s=0,1,2,\mathrm{}`$ and $`l=0,\pm 1,\pm 2,\mathrm{}`$ are radial and azimuthal quantum numbers.

The energy of a particle is written as

$$E_n=m_0c^2B_n=m_0c^2\sqrt{1+b_z^2+b_{}^2},b_{}=2\sqrt{\frac{\mu _0H}{m_0c^2}(n+\frac{1}{2})},$$

where $`b_z=p_z/m_0c`$ is the projection of the momentum onto the direction of the field and $`\mu _0=e_0\mathrm{}/2m_0c`$ is the Bohr magneton.

Matrix elements of the momentum operator of scalar particle have the form

$$\psi _m^{}|\widehat{P}_x|\psi _m=\frac{i}{2}m_0cb_{}(\delta _{m^{}m,+1}\delta _{m^{}m,1}),$$
(A3)
$$\psi _m^{}|\widehat{P}_y|\psi _m=\frac{1}{2}m_0cb_{}(\delta _{m^{}m,+1}+\delta _{m^{}m,1}),$$
(A4)
$$\psi _m^{}|\widehat{P}_z|\psi _m=m_0cb_z\delta _{m^{}m}.$$
(A5)
## B

In the paper we use the solutions to the Dirac-Pauli equation

$$i\mathrm{}\frac{\psi }{t}=\widehat{H}\psi ,\widehat{H}=c(𝜶\widehat{𝐏})+\rho _3m_0c^2+\mu _a\rho _3(𝝈𝐇),$$

where $`\widehat{𝐏}=\widehat{𝐩}(e/c)𝐀`$ is the kinetic momentum operator, $`𝐀=(Hy/2,Hx/2,0)`$, $`𝐇`$ is the external magnetic field, $`\mu _a=(\alpha /2\pi )\mu _0`$ is the anomalous part of electron magnetic moment, $`e=e_0<0`$ and $`m_0`$ are the charge and rest mass of electron.

It is known  that in this case complete set of commuting operators characterizing the quantum state of the particle consists of the Hamiltonian $`\widehat{H}`$, the projections of momentum operator $`\widehat{𝐏}`$ and total angular momentum operator $`\widehat{𝐉}=\widehat{𝐋}+\frac{1}{2}\mathrm{}𝝈`$ onto the direction of magnetic field vector which, for definiteness, is oriented along the $`Z`$-axis, and polarization operator. According to , in a homogeneous magnetic field the operator

$$\widehat{\mathrm{\Pi }}_z=\sigma _z+\frac{1}{m_0c}\rho _2(𝝈\times \widehat{𝐏})_z+\frac{\mu _aH}{c}$$

is taken as spin integral of motion characterizing the polarization of electron spin relative to the magnetic field vector. Then we have

$$\widehat{H}\psi =E_{n\zeta }\psi ,\widehat{P}_z\psi =p_z\psi ,\widehat{J}_z\psi =\mathrm{}(l\frac{1}{2})\psi ,$$
(B2)
$$\widehat{\mathrm{\Pi }}_z\psi =\zeta b\psi ,\zeta =\pm 1.$$
(B3)
In terms of cylindrical coordinates ($`r`$, $`\phi `$, $`z`$), the wave function being a solution to the Dirac-Pauli equation and satisfying Eqs. (B) has the form

$$\psi (𝐫,t)=\frac{\sqrt{e_0H}}{\sqrt{L}\sqrt{2\pi \mathrm{}c}}e^{iE_{n\zeta }t/\mathrm{}}e^{ip_zz/\mathrm{}}e^{i(l1)\phi }f(\rho ,\phi ),$$
(B4)
where

$$f(\rho ,\phi )=\left(\begin{array}{c}C_1I_{n1,s}(\rho )\hfill \\ iC_2I_{n,s}(\rho )e^{i\phi }\hfill \\ C_3I_{n1,s}(\rho )\hfill \\ iC_4I_{n,s}(\rho )e^{i\phi }\hfill \end{array}\right).$$

Here $`n=l+s=0,1,2,\mathrm{}`$ is the principle quantum number, $`s=0,1,2,\mathrm{}`$ is the radial quantum number, and $`l=0,\pm 1,\pm 2,\mathrm{}`$ is the azimuthal quantum number. In dimensionless form, the spin coefficients $`C_i`$ are determined as

$`C_1`$ $`=`$ $`{\displaystyle \frac{\zeta }{2}}\sqrt{{\displaystyle \frac{1}{2}}(1+\zeta {\displaystyle \frac{1}{b}})}\left[\sqrt{1+{\displaystyle \frac{b_z}{B_{n\zeta }}}}+\zeta \sqrt{1{\displaystyle \frac{b_z}{B_{n\zeta }}}}\right],`$ (B5)
$`C_2`$ $`=`$ $`{\displaystyle \frac{\zeta }{2}}\sqrt{{\displaystyle \frac{1}{2}}(1\zeta {\displaystyle \frac{1}{b}})}\left[\sqrt{1{\displaystyle \frac{b_z}{B_{n\zeta }}}}\zeta \sqrt{1+{\displaystyle \frac{b_z}{B_{n\zeta }}}}\right],`$ (B6)
$`C_3`$ $`=`$ $`{\displaystyle \frac{\zeta }{2}}\sqrt{{\displaystyle \frac{1}{2}}(1+\zeta {\displaystyle \frac{1}{b}})}\left[\sqrt{1+{\displaystyle \frac{b_z}{B_{n\zeta }}}}\zeta \sqrt{1{\displaystyle \frac{b_z}{B_{n\zeta }}}}\right],`$ (B7)
$`C_4`$ $`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{1}{2}}(1\zeta {\displaystyle \frac{1}{b}})}\left[\sqrt{1+{\displaystyle \frac{b_z}{B_{n\zeta }}}}+\zeta \sqrt{1{\displaystyle \frac{b_z}{B_{n\zeta }}}}\right],`$ (B8)
where the energy of electron is given by the expression

$$B_{n\zeta }=\frac{E_{n\zeta }}{m_0c^2}=\sqrt{b_z^2+\left(\sqrt{1+b_{}^2}+\zeta \frac{\mu _aH}{m_0c^2}\right)^2},$$
(B9)
$$b_z=\frac{p_z}{m_0c},b_{}=2\sqrt{\frac{\mu _0H}{m_0c^2}n},b=\sqrt{1+b_{}^2}.$$

Let us consider exact forms of all necessary matrix elements of the momentum and spin operators.

$$\psi _{m^{}\zeta ^{}}|\widehat{P}_x|\psi _{m\zeta }=\frac{i}{2}m_0cb_{}(\delta _{m^{}m,+1}\delta _{m^{}m,1})\delta _{\zeta ^{}\zeta },$$

$$\psi _{m^{}\zeta ^{}}|\widehat{P}_y|\psi _{m\zeta }=\frac{1}{2}m_0cb_{}(\delta _{m^{}m,+1}+\delta _{m^{}m,1})\delta _{\zeta ^{}\zeta },$$

$$\psi _{m^{}\zeta ^{}}|\widehat{P}_z|\psi _{m\zeta }=m_0cb_z\delta _{m^{}m}\delta _{\zeta ^{}\zeta },$$

$$\psi _{m^{}\zeta ^{}}|\widehat{S}_x|\psi _{m\zeta }=\frac{i}{2}(b\zeta (\delta _{m^{}m,+1}\delta _{m^{}m,1}))(\delta _{m^{}m,+1}\delta _{m^{}m,1})\delta _{\zeta ^{}\zeta },$$

$$\psi _{m^{}\zeta ^{}}|\widehat{S}_y|\psi _{m\zeta }=\frac{1}{2}(b\zeta (\delta _{m^{}m,+1}\delta _{m^{}m,1}))(\delta _{m^{}m,+1}+\delta _{m^{}m,1})\delta _{\zeta ^{}\zeta },$$

$$\psi _{m^{}\zeta ^{}}|\widehat{S}_z|\psi _{m\zeta }=(\zeta \frac{B_{n\zeta }}{b}\delta _{\zeta ^{}\zeta }+\frac{b_{}b_z}{b}\delta _{\zeta ^{}\zeta })\delta _{m^{}m},$$

$$\psi _{m^{}\zeta ^{}}|\widehat{S}^0|\psi _{m\zeta }=(\zeta \frac{b_z}{b}\delta _{\zeta ^{}\zeta }+\frac{B_{n\zeta }b_{}}{b}\delta _{\zeta ^{}\zeta })\delta _{m^{}m}.$$


text:
# 1 Energy level diagram

Robust quantum memory via quantum control

Andrew D. Greentree, S. G. Schirmer and A. I. Solomon

Quantum Processes Group, Departments of Physics and Astronomy and Applied Maths, The Open University, Milton Keynes, MK7 6AA, United Kingdom

Tel: +44-1908-652326, Fax: +44-1908-652140, S.G.Schirmer@open.ac.uk

## Abstract

A general scheme for building a quantum memory by transferring quantum information to an essentially decoherence-free memory transition using quantum control is presented and illustrated by computer simulations.

© 2024 Optical Society of America
OCIS codes: 270.0270, 270.1670

Quantum computation  has been a fruitful area of research lately. Some of the most promising schemes involve encoding qubits into ions and neutral atoms in high Q optical cavities. One of the greatest limitations to such schemes is the decoherence which occurs at the optical transitions. This decoherence is the limiting factor in determining the temporal length of a sequence of pulses to perform a given computation. Many schemes have been suggested to overcome this decoherence, especially quantum error correction , which uses redundant information to compensate for the losses, and decoherence-free substates , which use combinations of states which are robust against decay because of quantum interference. In this paper we employ Lie group decompositions  to derive a promising scheme for a quantum memory, with the idea being to transfer quantum information from the important channel for quantum information processing to a “memory transition” which holds the quantum information.

Our scheme can be illustrated with a simple example. Consider a four-level atom as depicted in figure 2 with two degenerate ground states $`|g_1`$ and $`|g_2`$ and two non-degenerate excited states $`|e_1`$ and $`|e_2`$. The ground states might be for example Zeeman sublevels of an atom such as Rubidium, although our scheme is far more general than this. We assume that the $`|g_1|e_1`$ transition is an optical transition that forms part of a quantum information processing scheme. After performing some quantum logical operations we wish to protect the quantum information stored in this transition in a long-lived system by applying a series of Gaussian control pulses (derived from optical fields) to transfer the quantum information onto the decoherence-free $`|g_1|g_2`$ transition. Although the transition pairs $`|g_1|e_1`$, $`|g_2|e_1`$ and $`|g_1|e_2`$, $`|g_2|e_2`$ are at equal frequencies, they can be individually addressed through the choice of appropriate field polarisations, allowing complete controllability of the system.

Formally, the problem can be stated as follows. We wish to map the density matrix $`\rho `$ representing the (initial) state of the system, whose elements are $`\rho _{g_1g_1}`$, $`\rho _{e_1e_1}`$, $`\rho _{g_1e_1}`$ and $`\rho _{g_2g_2}=\rho _{e_2e_2}=\rho _{g_1g_2}=\rho _{g_1e_2}=\rho _{g_2e_1}=\rho _{g_2e_2}=0`$, onto a density matrix $`\rho ^{}`$ such that $`\rho _{g_1g_1}^{}=\rho _{g_1g_1}`$, $`\rho _{g_2g_2}^{}=\rho _{e_1e_1}`$, $`\rho _{g_1g_2}^{}=\rho _{g_1e_1}`$ and $`\rho _{e_1e_1}=\rho _{e_2e_2}=\rho _{g_1e_1}=\rho _{g_1e_2}=\rho _{g_2e_1}=\rho _{g_2e_2}=0`$ by applying a sequence of simple control pulses. Because the populations and coherence have been mapped to degenerate energy levels with no allowed transitions between them, these states will be extremely long lived. The quantum information can then be returned to the optical transition simply by using the reverse of the quantum control scheme. In order to realize the mapping of $`\rho `$ onto $`\rho ^{}`$ we find a unitary operator $`U`$ such that $`\rho ^{}=U\rho U^{}`$ and decompose the operator $`U`$ into a product of simple unitary operators, each realized dynamically by applying a control field 1 (2) (3) that drives the $`|g_2|e_2`$ ($`|g_1|e_2`$) ($`|g_1|e_1`$) transition, respectively. Concretely, note that we have

$$\underset{\rho ^{}}{\underset{}{\left(\begin{array}{cccc}\rho _{g_1g_1}& \rho _{g_1e_1}& 0& 0\\ \rho _{e_1g_1}^{}& \rho _{g_2g_2}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)}}=\underset{U}{\underset{}{\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right)}}\underset{\rho }{\underset{}{\left(\begin{array}{cccc}\rho _{g_1g_1}& 0& \rho _{g_1e_1}& 0\\ 0& 0& 0& 0\\ \rho _{e_1g_1}^{}& 0& \rho _{e_1e_1}& 0\\ 0& 0& 0& 0\end{array}\right)}}\underset{U^{}}{\underset{}{\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right)}}.$$
(1)
Furthermore, we can express $`U`$ as a product $`U=V_1V_2V_3V_2V_1`$ where

$$V_1=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\end{array}\right),V_2=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right),V_3=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right).$$
(2)
Observing that

$$V_1=\mathrm{exp}\left[\frac{\pi }{2}\left(|g_2e_2||e_2g_2|\right)\right]$$
(3)
we see that the operator $`V_1`$ can be dynamically realized by applying a pulse $`f_1(t)=A_1(t)e^{i\omega _{g_2e_2}+\pi /2}`$ with total pulse area $`|A_1(t)|𝑑t=\pi /(2d_{g_2e_2})`$ and appropriate polarisation, where $`\omega _{g_2e_2}`$ is the frequency and $`d_{g_2e_2}`$ the absorption oscillator strength of the $`|g_2|e_2`$ transition. Similarly,

$$V_2=\mathrm{exp}\left[\frac{\pi }{2}\left(|g_1e_2||e_2g_1|\right)\right],V_3=\mathrm{exp}\left[\frac{\pi }{2}\left(|g_1e_1||e_1g_1|\right)\right]$$
(4)
shows that $`V_2`$ \[$`V_3`$\] can be dynamically realized by applying a pulse $`f_2(t)=A_2(t)e^{i\omega _{g_1e_2}+\pi /2}`$ \[$`f_3(t)=A_3(t)e^{i\omega _{g_1e_1}+\pi /2}`$\] with appropriate polarisation and total pulse area $`|A_2(t)|𝑑t=\pi /(2d_{g_1e_2})`$ \[$`|A_3(t)|𝑑t=\pi /(2d_{g_1e_1})`$\] where $`\omega _{g_1e_2}`$ \[$`\omega _{g_1e_1}`$\] is the frequency and $`d_{g_1e_2}`$ \[$`d_{g_1e_1}`$\] the absorption oscillator strength of the $`|g_1|e_2`$ \[$`|g_1|e_1`$\] transition.

Hence, generation of $`U`$ requires a sequence of five pulses, where the first pulse drives the transition $`|g_2|e_2`$ and has pulse area $`\pi /(2d_{g_2e_2})`$; the second pulse drives the transition $`|g_1|e_2`$ and has pulse area $`\pi /(2d_{g_1e_2})`$; the third pulse drives the transition $`|g_1|e_1`$ and has pulse area $`\pi /(2d_{g_1e_1})`$; the fourth pulses drives again the transition $`|g_1|e_2`$ and has pulse area $`\pi /(2d_{g_1e_2})`$ and the last pulse drives again the transition $`|g_2|e_2`$ and has pulse area $`\pi /(2d_{g_2e_2})`$. Observe that only the total area of each pulse matters and hence the precise pulse envelopes can be adjusted to suit experimental constraints.

Figure 2 shows the pulse sequence and evolution of a system initially in the superposition state $`(|g_1+2|e_1)/\sqrt{5}`$. Note that initially $`\rho _{g_1g_1}=0.2`$, $`\rho _{e_1e_1}=0.8`$, $`\rho _{g_1e_1}=0.4`$ and all other matrix elements of $`\rho `$ are zero. Observe that at the final time we have indeed $`\rho _{g_1g_1}^{}=0.2`$, $`\rho _{g_2g_2}^{}=0.8`$, $`\rho _{g_1g_2}=0.4`$ and all other elements of $`\rho ^{}`$ are zero. The time unit in all plots is $`1/\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ is the Rabi frequency of the transition $`|g_1|e_1`$.

One of the authors (ADG) would like to acknowledge the financial support of the EPSRC.


text:
# Solid-State Quantum Computer Based on Scanning Tunneling Microscopy

## I Introduction

Recently, we suggested a solid-state nuclear spin quantum computer based on magnetic resonance force microscopy (MRFM) . In this proposal, a ferromagnetic particle, which is attached to the tip of the cantilever, moves along a chain of paramagnetic ions (atoms) located below the surface of silicon. The ferromagnetic particle targets a particular ion, providing selective excitation of its electron or nuclear spin, a measurement of the state of the nuclear spin, and initial polarization and one-qubit rotations of nuclear spins. A two-qubit Control-Not (CN) gate for two neighboring nuclear spins utilizes the dipole-dipole interaction between the nuclear spin of one ion and the electron of the neighboring ion.

There are three basic disadvantages in the suggested MRFM proposal: 1) The motion of the ferromagnetic particle along the chain of paramagnetic ions can induce additional dephasing in the spin system. This dephasing must be taken into account. 2) Fluctuations of the position of a ferromagnetic particle cause additional decoherence in a spin system. 3) The accurate positioning of the ferromagnetic particle relative to an ion is a complicated experimental problem.

In this paper, we suggest a quantum computer architecture which seems to be free of all these disadvantages. At the same time, it retains the main advantage of the MRFM proposal – it does not require an array of single-atom gates, which is an essential part of well-known solid-state proposals based on silicon technology . We give estimates for characteristic parameters required for realization of a quantum computer based on the STM technique.

## II STM detection of an electron spin resonance and a nuclear spin state

Our proposal utilizes recent results in STM experiments which seem to allow detection of the precessing electron spin from a single ion (atom). In these experiments , the STM technique was used to measure the frequency of the Larmor precession of an electron spin of an individual iron atoms in silicon in the presence of a small, static applied magnetic field. This measurement utilizes the interaction between a localized electron spin in a surface of a sample and the current produced by tunneling electrons. Due to this interaction the tunneling probability is affected by the Larmor precession of a localized spin, and the tunneling current is modulated at the Larmor frequency.

The mechanism causing the tunnel current modulation is, however, not well understood. While this requires additional consideration, in this paper we assume that the STM technique is able to detect the frequency of electron spin precession for an individual ion (atom) in a silicon surface. This results in the possibility for the STM technique to detect the state of a nuclear spin in the same ion. Indeed, suppose that the hyperfine interaction between electron and nuclear spins of the ion is greater than the Zeeman interaction between a high external magnetic field, $`\stackrel{}{B}_0`$, and the nuclear spin. Then, the frequency of an electron spin precession depends on the state of the nuclear spin. As an example, following , we consider a $`{}_{}{}^{125}Te^+`$ ion (A-center) in the surface of silicon. The Hamiltonian of the electron-nuclear spin system of the ion has the form,

$$=g_e\mu _BB_0S_z+g_n\mu _nB_0I_zA\stackrel{}{S}\stackrel{}{I},$$
$`(1)`$
where $`\stackrel{}{S}`$ and $`\stackrel{}{I}`$ are electron and nuclear spin operators ($`S=I=1/2`$); $`B_0`$ is a high magnetic field pointed in the positive $`z`$-direction; $`g_e2`$, $`g_n0.882`$; and $`A/2\pi \mathrm{}3.5`$GHz is the constant of the hyperfine interaction. The magnetic moments of the $`{}_{}{}^{125}Te`$ nucleus and the electron are both negative. Fig. 1, demonstrates the energy levels of the electron and nuclear spins of the $`{}_{}{}^{125}Te^+`$ ion in a high magnetic field.

If the nuclear magnetic moment points “up”, the modulation frequency is,

$$f=f_{e0}=g_e\mu _BB_0/2\pi \mathrm{}+A/4\pi \mathrm{}.$$
$`(2)`$
In the opposite case,

$$f=f_{e1}=g_e\mu _BB_0/2\pi \mathrm{}A/4\pi \mathrm{}.$$
$`(3)`$
The hyperfine splitting is $`A/2\pi \mathrm{}`$. If the tunneling current is observed to oscillate with the frequency $`f_{e0}`$, the nuclear spin of the ion is in its ground state. Otherwise the nuclear spin is in the excited state.

## III Selective excitation of
a nuclear spin

To provide a selective excitation of the nuclear spin (qubit), we assume that a permanent magnetic field is slightly non-uniform in the $`x`$-direction (see Fig. 2). Assuming that the distance between the neighboring ions is $`a`$, we have, for the difference in the precession frequencies for these ions,

$$\mathrm{\Delta }f_e=(g_e\mu _B/2\pi \mathrm{})a\frac{B_0}{x}.$$
$`(4)`$
The corresponding difference for nuclear magnetic resonance (NMR) frequencies is,

$$\mathrm{\Delta }f_n=(g_n\mu _n/2\pi \mathrm{})a\frac{B_0}{x}.$$
$`(5)`$
(In Eq. (4) we assume that nuclear spins of the neighboring ions are in the same state. In Eq. (5) we assume that the electron spins of the neighboring ions are in the same state.) As an example, taking $`a=5`$nm and $`B_0/x=10^5`$T/m, we obtain for 125Te: $`\mathrm{\Delta }f_n6.75`$kHz.

For selective excitation of a nuclear spin, one can apply a $`\pi `$\- or $`\pi /2`$-pulse of a rotating magnetic field, $`\stackrel{}{B}_{}`$. The nuclear Rabi frequency, $`f_{nR}`$, is given by,

$$f_{nR}=(g_n\mu _n/2\pi \mathrm{})B_{}.$$
$`(6)`$
For selective excitation, this quantity must be less than $`\mathrm{\Delta }f_n`$. It imposes a restriction on the duration, $`\tau `$, of a $`\pi `$-pulse,

$$\tau =\pi /2\pi f_{nR}>1/2\mathrm{\Delta }f_n74\mu s.$$
$`(7)`$
Using selective nuclear $`\pi `$-pulses one can drive a nuclear spin chain into the ground state. For this purpose one prepares the system in a high magnetic field at a low temperature, so that all electron spins are in their ground state. (As an example, $`B_010`$T, $`T1`$K.) Then, one measures the modulation frequency of the tunneling current at any position, $`x_k`$. If the modulation frequency fits the value, $`f_{e1}(x_k)`$, one applies a selective nuclear $`\pi `$-pulse with the frequency, $`f_n(x_k)`$, which drives the nuclear spin into the ground state.

## IV A Control-Not gate

The ability to implement a CN gate is important for quantum computation as any quantum algorithm can be decomposed into a sequence of one-qubit rotations and CN gates . We propose to use a scheme close to that outlined in our MRFM proposal . It consists of three steps:

1) One applies an electronic $`\pi `$-pulse with the frequency $`f=f_{e1}(x_k)`$. This pulse drives the electron spin of the ion “k” into the excited state only if the nuclear spin of the ion “k” (a control qubit) is in the excited state. The hyperfine splitting of the ion must be smaller than the difference, $`\mathrm{\Delta }f_e`$, of precession frequencies of neighboring ions. For $`B_0/x=10^5`$T/m we have: $`\mathrm{\Delta }f_e=14`$MHz. So, this condition is satisfied, and the frequency, $`f_{e1}(x_k)`$, is unique in the spin chain. (The dipole-dipole contribution to the oscillation is small in comparison to the hyperfine splitting.)

2) One tunes the frequency of the nuclear $`\pi `$-pulse for the ion “k+1” (or “k-1”) taking into consideration the dipole field produced by an electron spin on the nuclear spin of the ion “k+1” (a target qubit). If the electron spin of the ion “k” did not change its direction, the NMR frequency for the ion “k+1” (taking into account a dipole field of electron spins) is,

$$f=f_n(x_{k+1})f_{nd}f_{nd}^{},$$
$`(8)`$
where $`f_{nd}`$ is the dipole shift due to two neighboring electron spins, and $`f_{nd}^{}`$ is the dipole shift due to all other electron spins. (For $`a=5`$nm, $`f_{nd}200`$Hz, $`f_{nd}^{}<40`$Hz, the nuclear-nuclear dipole interaction is negligible.) If the electron spin of the ion “k” has changed its direction, then the NMR frequency for the ion “k+1” is,

$$f=f_n(x_{k+1})f_{nd}^{}.$$
$`(9)`$
The difference, $`f_{nd}`$, between the frequencies (8) and (9) must be smaller than $`\mathrm{\Delta }f_n`$. In our case this condition is satisfied, so both frequencies (8) and (9) are unique in a spin chain. If the frequency of a nuclear $`\pi `$-pulse is tuned to (9), and the Rabi frequency, $`f_{nR}`$, of the nuclear $`\pi `$-pulse is less than $`f_{nd}`$, then the nuclear spin of the ion “k+1” changes its direction only when the control qubit is in the excited state. This provides an implementation of a CN gate. This step imposes a most severe restriction on the duration of a quantum gate: $`\tau =\pi /2\pi f_{nR}>2.5`$ms. The same step imposes a restriction on the electron relaxation time which must be greater than $`\tau `$. (Available data  show that the electron relaxation time for paramagnetic impurities can be much greater than $`2.5`$ms.)

3) Finally, one repeats the first step to drive an electron spin of the ion “k” back to the ground state, if the control qubit is in the exited state.

## V Conclusion

We proposed a nuclear spin solid-state quantum computer based on scanning tunneling microscopy (STM). Our proposal does not require sophisticated single-atom electrostatic gates. A single-spin measurement is provided by measuring the modulation of the tunneling current. Selective excitation of a nuclear spin can be achieved using an inhomogeneous permanent magnetic field. A Control-Not gate is provided by utilizing the dipole-dipole interaction between the electron and nuclear spins. Our envisioned STM quantum computer would operate at high magnetic field ($`10`$T) and at low temperature $`1`$K.
Acknowledgments
This work was supported by the Department of Energy (DOE) under contract W-7405-ENG-36. The work of GPB and VIT was partly supported by the National Security Agency (NSA) and by the Advanced Research and Development Activity (ARDA).


text:
# Calculation of atmospheric neutrino flux

## 1 Introduction

The interpretation of the data on the atmospheric neutrinos obtained by different experiments  strongly relies on the calculations of neutrino fluxes . Starting from a common assumption that the primary cosmic rays arrive in the vicinity of the Earth isotropically the calculations differ in the level of sophistication used by different authors to describe the production of the secondaries in the hadronic interactions of primary cosmic particles (mostly protons) with the atmosphere and in the way the secondaries are propagated in the Earth’s magnetic field. Different calculations producing considerably different results, the validity of neutrino flux calculations can only be checked by comparing the predictions of these calculations with the experimental data. Although the difference in the measured and calculated atmospheric fluxes of upgoing and downgoing neutrinos is considered as a clear evidence for neutrino oscillations  it is nevertheless recognized , that the same calculations  seem to fail to properly describe the measured east-west asymmetry of atmospheric neutrino fluxes  related to a well known  observation that due to the influence of the Earth’s (nearly dipole) magnetic field the fluxes of the primary predominantly positive (p,$`\alpha `$-particles) hadronically interacting cosmic rays reaching the atmosphere contain more particles coming from the west than from the east. Thus, more accurate and detailed calculations describing propagation of primary cosmic rays and their interactions with the atmosphere is needed.

The recent data from the Alpha Magnetic Spectrometer (AMS)  provide a vast amount of information on the fluxes of both primary and secondary particles detected in the near Earth orbit covering most of the Earth’s surface. The AMS data display not only the common features of the magnetic cutoff dependence on the distance from the magnetic equator but also reveal the presence of particles well under the cutoff . These particles originate from the production and eventual decay of the secondaries produced in the inelastic interactions of the primary cosmic particles with the atmosphere and consequently carry the information about neutrinos produced in decays of the secondaries. A comparison of the results obtained in the calculations with the fluxes of the secondaries measured with AMS represents a means to prove or disprove the validity of the results of the calculations concerning the predictions of atmospheric neutrino fluxes.

In the following the result of calculation of the near Earth fluxes of charged particles and atmospheric neutrinos is presented. First, the model used and its parameters is described. Section 3 describes the origin of the particles near the Earth. Section 4 is devoted to results of calculations of fluxes of charged particles at the altitude of 400 km from the Earth’s surface corresponding to the AMS flight altitude. The simulated fluxes of charged particles are compared with those measured with AMS in section 5. Section 6 is devoted to the results of calculations of atmospheric neutrino fluxes at the ground level in general. The neutrino fluxes at different experimental sites are discussed in section 7. The last section gives a summary of the results and draws some conclusions.

## 2 Description of the model

The calculation uses the Monte Carlo technique and is based on the GEANT3/GHEISHA package . This package is extensively used for the simulation of processes in various detectors of high energy physics experiments. To use this package for the simulation on a scale of propagation of cosmic rays in the Earth’s magnetic field, their interactions with the Earth’s atmosphere and the Earth itself, followed by production, propagation and decay of secondary particles, the parameters of the corresponding media are adapted to scale them from a thousands of kilometer scale of cosmic ray travel and interaction to a centimeter scale characteristic for GEANT based simulations.

The Earth is modeled to be a sphere of 6378.14 km radius of a uniform density of 5.22 g/$`\mathrm{cm}^3`$ derived from the total mass of the Earth. The atomic and nuclear properties of the Earth are taken to be those of Ge - the closest in density (5.32 g/$`\mathrm{cm}^3`$) element. The Earth’s atmosphere is modeled by 1 km thick variable density layers of air extending up to 71 km from the Earth’s surface. The density change with altitude is calculated using a Standard Atmosphere Calculator .

The Earth’s magnetic field for the year 1998 is calculated according to World Magnetic Field Model (WMM-2000)  with 6 degrees of spherical harmonics. The WMM-2000 is a predictive model based on the magnetic field measurements on a global scale.

The flux ($`\mathrm{\Phi }`$) of high energy primary protons in the Solar system is parametrized on the basis of the existing experimental data  by a power low in rigidity (R), $`\mathrm{\Phi }(R)R^\gamma `$ with $`\gamma `$=2.79. At the energies below several GeV the spectrum is corrected for the solar activity to be similar to that of the year 1976  \- the year of a similar solar activity as in the year 1998 according to the 11 year solar cycle. Primary protons are isotropically emitted from a sphere of 10 Earth’s radii. The long distance from the Earth of the emission point is sufficient for a particle to endure sufficient influence of the Earth’s magnetic field to form an adequate topology of the flux in the Earth’s vicinity.

Secondary particles produced in the interactions with the atmosphere and the Earth are traced in the magnetic field following a normal GEANT interaction/decay procedure.

## 3 Origin of particles near Earth

The shape of magnetic field determines the behaviour of charged particles at their approach to the Earth and the eventual motion of primaries and secondaries in the Earth’s vicinity. The channeling of the particles in the Earth’s magnetic field leads to a well known effect that the flux of particles coming from space is at its maximum at the magnetic poles (see top Fig.1). This effect is associated with the magnetic cutoff of low momenta particles in the equator region. A peculiar feature of particle channeling from the outer space is that most of the particles detected in the pole region in the Earth’s vicinity have their origin in the equator region of the isotropically shining sphere at 10 Earth’s radii used in the simulation (bottom Fig.1). The efficiency of channeling of the particles depends on the emission parameters such as particle momentum and emission direction.

Figure 2 shows the emission parameters of the particles reaching the Earth. One can see in figure 2 that small momentum particles even emitted in the direction opposite to the Earth can reach the Earth’s vicinity whereas high momentum particles can do that only when the initial direction is close to a straight line between the emission point and the Earth. The dashed lines in figure 2 show the cuts on the emission parameters used in the simulation. An additional cut is set on the momentum of primary protons. It is set to be of 0.5 GeV/c for two reasons. First, due to the detector material thickness the minimum kinetic energy required for a proton to trigger AMS is $``$ 0.13 GeV corresponding to a proton momentum of 0.5 GeV/c. Secondly, the kinetic energy threshold of charged pion production in p-nucleus collisions is $``$ 0.16 GeV  corresponding to a proton momentum of 0.54 GeV/c. Protons with energies below the threshold can not produce charged pions which represent the source of detectable secondaries and eventually of atmospheric neutrinos. The cut on the proton momentum eliminates the few events in the region forbidden by the cut of the bottom Fig.2.

Of all protons with the momentum spectrum described above uniformly emitted from the sphere of 10 Earth’s radii only a $`7.410^4`$ part reaches the Earth’s atmosphere. In the subsequent simulation only particles with initial kinematics compatible with those reaching the Earth’s atmosphere are traced the whole way until they interact with the atmosphere or the Earth itself or/and continue the journey beyond a sphere of 15 Earth’s radii. The charged secondaries from the interactions are traced until they go below a limit of 0.125 GeV/c.

Figure 3 shows an example of the trajectory of a particle reaching the Earth. (Note that the AMS flight altitude of 400 km corresponds to 1.06 Earth’s radii). Due to influence of the Earth’s magnetic field the direction of a particle at its entrance to the atmosphere is by far different from the initial direction.

Primary cosmic particles may endure an elastic or inelastic scattering in the atmosphere thus changing their direction and/or momentum. Figure 4 shows the trajectory of such a particle which eventually reenter the atmosphere but at a place far away from the point of initial interaction. As it will be shown later, there is a considerable amount of such particles in the total flux of incoming cosmic rays.

Secondary particles produced in the inelastic interactions of primary cosmic rays with the atmosphere follow a usual pattern of hadronic or electromagnetic cascade and most of them quickly end up in the atmosphere itself or in the Earth. However, under the influence of the Earth’s magnetic field some (predominantly stable) particles from the atmospheric showers escape upward from the atmosphere and travel long distances before being absorbed in the atmosphere. The trajectory of such particle depends on the kinematics of the particle production and the geomagnetic coordinates of the origin. Unless they interact with the atmosphere the particles with the same rigidity (p/z) follow the same trajectory. Figure 5 shows the trajectory of a positive particle produced in the atmosphere in the northern hemisphere, reflected at a mirror point in the southern hemisphere before being absorbed at an approximatively the same latitude but a different longitude.

The spiraling of the particles, their motion from and to the Earth’s surface as well as their bouncing across and drifting along the magnetic equator is determined by the shape of the Earth’s magnetic field. Depending on the production place and direction with respect to the Earth’s surface secondary particles can make a return trip as illustrated in Fig.5, a one way trip across the magnetic equator, a multiple bounce (Fig.6) or even a round the globe trip (Fig.7). At the limit a secondary particle produced in the atmosphere over a certain limited region of the Earth (called South Atlantic Anomaly) can have a closed trajectory (Fig.8) and could live forever if it would not interact with the atmosphere at the place where it was created nor lose its energy on synchrotron radiation.

As one can see in figures 5d, 6d, 7d and 8d the particles in their motion pass the level of the AMS flight (1.06 Earth’s radii) going up and down several times. The geographic coordinates of the place over which a secondary particle can be detected at a certain altitude is connected with the place of origin of the particle.

One of the consequences of the fact that in its motion the same particle can cross a detector several times is that the measured flux of such upgoing and downgoing particles is the same. A build up of the measured rate, especially in the region close to the magnetic equator, connected with the multiple crossing of the detector plane rate, does not produce secondary particles unless the spiraling particle passes through and interacts with the atmosphere. The motion of the particle being determined by the Earth’s magnetic field, negative particles travel along the same trajectories as the positive ones but in the opposite direction.

The production of secondary pions and protons resulting from the interactions of cosmic (mostly He) nuclei with the atmosphere is treated using the superposition approximation .

The calculation presented below takes into account all effects connected with production, propagation, interaction and decay of all particles of the atmospheric shower.

## 4 Particle flux

The results of Monte Carlo simulation presented in the following correspond to a statistics of $`1.510^{10}`$ protons isotropically emitted from a sphere of 10 Earth’s radii. This is equivalent to a 0.62 ps exposure.

The dependence on the geographical coordinates of the rate of protons crossing a sphere surrounding the Earth at the altitude of 400 km is shown in Figure 9 for downward (top picture) and upward (bottom picture) protons. The rate includes all i.e. primary and secondary, protons and is averaged over all directions in the corresponding hemispheres.

A typical picture with high rate near the magnetic poles and the depleted by the magnetic cutoff regions near the magnetic equator is clearly seen in the downward proton rate behaviour. The geographical pattern of upgoing protons is similar in shape but at much lower particle rate.

The flux of protons going toward the Earth (Figs.10a and b) consists of original protons from space i.e. those with $`p/p_{init}>`$0.9999996, where $`p_{init}`$ is the initial cosmic proton momentum and $`p`$ is the momentum of the proton crossing the detection sphere (see Fig.3), protons scattered in the atmosphere (0.995$`<p/p_{init}<`$0.9999996. Fig.4) and secondary protons ($`p/p_{init}<`$0.995. Figs.5,6,7,8). The flux of upgoing protons (Figs.10c and d) contains reflected protons most of which have lost some energy in the atmosphere as well as scattered protons on their way from the Earth and the secondary protons.

The important feature of the primary flux relevant to neutrino production is that up to momenta as high as $``$40 GeV it becomes anisotropic for directions close to parallel to the Earth’s surface. (see Fig.11). This decrease is not in contradiction with the AMS measurement because the orientation of the AMS detector during the flight did not allow to measure the horizontal flux.

Figure 12 illustrates the behaviour of the flux of primary protons. Figure 12a shows the dependence of the cutoff on the magnetic latitude. The cutoff momentum is about 10 GeV/c at the equator and goes down to zero at the poles. One can see, however, that there are some under cutoff particles even in the flux of primary protons. These under cutoff particles are present at all magnetic latitudes (Figs.12 a and d). The other feature of these protons is that they are moving predominantly, and in the magnetic equator region exlusivly (Figs.12b and c), toward the Earth.

The scattered protons (Fig.13) are mostly detected in the regions of the flux maxima close to the magnetic poles (Fig.13a). These are mostly albedo protons (Figs.13b and d) although there is a considerable amount of those (Figs.13b and c) reentering the atmosphere (see fig.4).

Secondary protons (Fig.14) produced at the magnetic poles add up to the albedo flux (Fig.14b) whereas those produced in the regions closer to the equator follow the trajectories illustrated in Figs.5,6,7, 8 producing equal fluxes of upgoing and downgoing undercutoff particles (Fig.14b). The impact points of particles with multiple crossing of the detection sphere are clearly seen in Fig.14a. The rate of under cutoff protons in the magnetic equator region is dominated by the long traveling protons (see Figs.6,7).

Positrons and electrons resulting from the decays of secondary pions display the behaviour similar to that of secondary protons allowing for much faster energy losses in the atmosphere due to the difference in radiation and nuclear interaction lengths.

## 5 Particles detected by zenith facing AMS

To compare the results of calculations with the AMS measurements  the overall fluxes produced in the simulation were restricted to the AMS acceptance during the zenith facing flight period. . In accordance with the AMS construction no particle with kinetic energy below $`120`$ MeV can trigger AMS. The AMS orbit spans within the geographic latitude of $`\pm 51.7^o`$. The corresponding cuts are applied on the MC simulated data. The result of Monte Carlo is an absolute one. No renormalization is applied to make a comparison with the data. As in AMS publications the spectra are presented as a function of kinetic energy of the particles. Figure 15 shows the spectra of protons obtained in the simulation for different positions of AMS with respect to the magnetic equator. The result of simulation correctly reproduces the characteristic features present in the AMS data (black dots in the figure) such as:

* the value of the cutoff at different magnetic latitudes;
* the presence of protons with under cutoff momenta at all magnetic latitudes;
* the equality of fluxes of under cutoff protons going to and from the Earth;
* the energy spectra of under cutoff protons and their absolute flux.

Figures 16 and 17 show the spectra of secondary positrons and electrons correspondingly. Although the positrons and electrons can not contribute to atmospheric neutrino production they are (together with cosmic muons) the ”spectators” of neutrino production in the decays of particles in the atmospheric showers and consequently the correct prediction of the fluxes of these particles is a proof of validity of the whole simulation procedure including the neutrino flux predictions. The simulated positron and electron spectra are similar to those obtained in the AMS experiment  for the secondary leptons in terms of equality of fluxes of upgoing and downgoing particles and energy dependence of the flux at different magnetic latitudes.

The difference in the measured and simulated secondary leptonic (and to some extent secondary proton) fluxes at low magnetic latitudes (Figs.16a, 17a and 15a) is due to insufficient statistics of the present Monte Carlo study. As mentioned above, the measured rate in this region is largerly connected with the relatively rare long traveling particles multiply crossing the AMS orbit plane. Actually, the peak in the spectrum of figure 16a is due to just 3 spiraling positrons each crossing the AMS detector plane about 500 times in both upward and downward direction. Nevertheless the measured relative rate of secondary leptons  and their dependence on the magnetic latitude is well reproduced by the simulation (Fig.18). It is worth to emphasize that the impact of long traveling secondaries on the production of atmospheric neutrinos is negligible.

Apart from the data connected with those measured by AMS the present simulation includes a prediction of the fluxes of muons at the sea level. The histogram in Fig.19 shows a result obtained in this study for the vertical muon flux. The result is in good agreement with the measured fluxes reported in .

## 6 Atmospheric neutrino fluxes

As for the fluxes of charged particles discussed in the previous section all results that follow are obtained directly without any renormalization. The dependence of the neutrino flux averaged over the entire Earth’s surface on the neutrino direction with respect to nadir is presented in Fig.20a. The overall equality of neutrino fluxes to and from the Earth is clearly seen. The dashed line in the figure is a mirror image of the flux from the top. A small increase of the flux of almost horizontal neutrinos coming from inside the Earth is due to neutrinos produced by the cosmic particle’s showers in the Earth. The relative flux of different kinds of atmospheric neutrinos ($`\nu _\mu ,\stackrel{~}{\nu }_\mu ,\nu _e,\stackrel{~}{\nu }_e`$) is constant for all incidence angles (Fig.20b.). The averaged over all directions energy spectra of neutrino of different kind are presented in Fig.20c. The dependence of the flux on the azimuth (with respect to east direction) angle (Fig.20d) reflects some east-west asymmetry of average neutrino fluxes.

The actual picture of neutrino flux distribution is very different from the average one. Fig.21 shows the geographic maps of neutrino fluxes measured in KHz/$`\mathrm{m}^2\mathrm{sr}`$ for downgoing (top) and upgoing (bottom) atmospheric neutrinos. The downward flux is very nonuniform and, not surprisingly, strongly correlated with the flux of cosmic protons (see fig.9). The downward neutrino flux changes from more than 15 KHz/$`\mathrm{m}^2\mathrm{sr}`$ at the magnetic poles down to less than 3 KHz/$`\mathrm{m}^2\mathrm{sr}`$ at the magnetic equator whereas the upward flux is almost uniform at the level of 6-7 KHz/$`\mathrm{m}^2\mathrm{sr}`$ over the entire Earth’s surface. For different regions of the Earth not only the total flux of neutrinos is different but also the kinematical parameters of the neutrinos.

To allow the comparison with previous calculations of Ref. we compute the neutrino event rate in three equal area regions according to the geomagnetic latitude, although in view of neutrino flux rapidly changing with magnetic latitude (see Fig.21) such division appears to be rather crude.

As in  to compute the event rates the $`\nu `$ fluxes have been convoluted with the neutrino cross section ,,. The angle in the distributions (Figs.22,23,24,25,26,27) is the neutrino one with no experimental smearing or inefficiency included. Apparently, the results of the present study differ from previous calculations both in the energy spectra and the angular distributions of atmospheric neutrinos. The flux of sub-GeV neutrino predicted by this calculation depends much stronger on the geographical region following the intensity of the primary cosmic particle flux (Fig.9). The down/up asymmetry of the flux is much stronger. There is an important increase in the horizontal flux emerging from the Earth. This flux is connected with neutrinos produced under the Earth’s surface at a depth of several meters. A part of these neutrinos can be detected in the underground detectors. The actual amount of ”underground” neutrinos detected in a given detector depends on the detector position depth.

## 7 Atmospheric neutrino fluxes at major experimental sites

For the relatively small area corresponding to different experimental sites the statistics of the present study is high enough to make reliable conclusions for neutrinos with energies below $``$1.5 GeV. The parameters describing the neutrino flux at the Kamioka site are calculated for an area within the geographical coordinates $`35^o\pm 12^oN,142^o\pm 12^oE`$ for an underground detector at a depth of 520 m (2700 mwe). The averaged over all directions differential energy spectra of different neutrinos and antineutrinos are given in Fig.28 together with the results from other authors. The present study gives lower neutrino rates at higher energies. The fluxes of all kinds of neutrino are strongly down/up asymmetric with a considerable contribution of ”underground” neutrinos in the horizontal direction (Fig.29b). The ratio of muon($`\nu _\mu +\stackrel{~}{\nu }_\mu `$) to electron($`\nu _e+\stackrel{~}{\nu }_e`$) fluxes (Fig.29c) suggests a tendency to decrease for horizontal directions.

This systematic effect for neutrinos with the highest flux (see Fig.29b) if not taken into account may lead to a conclusion about $`\nu _\mu ,\stackrel{~}{\nu }_\mu `$ disappearance. The east-west asymmetry of electron and muon neutrino is visible in Fig.29d.

The energy dependence of the angular distribution of the flux of electron and muon neutrino is illustrated by Figs.30a and 30b correspondingly. Note the absence of the increase in the muon neutrino flux in the horizontal direction. Figure 30d shows the energy dependence of the downward/upward flux ratio for muon (solid line) and electron (dashed line) neutrinos. The ratio is the same for two different kinds of neutrino in the energy range below 0.31 and above 1.5 GeV. However, for the region between 0.31 and 1.5 GeV where the bulk of the experimental data is collected the down/up ratio for muon neutrinos is systematically higher than that for electron neutrinos. This systematic effect may fake muon neutrino disappearance.

Neutrino fluxes at the Soudan and Gran Sasso sites are calculated for the areas with geographical coordinates of $`48^o\pm 12^oN,98^o\pm 12^oW`$ and $`42^o\pm 12^oN,13^o\pm 12^oE`$ and the depth underground of 700 m and 750 m correspondingly. As for the Kamioka site the averaged over all directions differential energy spectra of different neutrinos and antineutrinos at the Soudan (Fig.31) and Gran Sasso (Fig.34) sites obtained in the present calculation give lower neutrino rates at higher energies than those obtained by other authors. The angular distributions of neutrino flux is strongly up/down asymmetric with predominance of the downward flux at the Soudan site (Fig.32b) whereas it is almost up/down symmetric for Gran Sasso (Fig.35b). The value and the energy dependence of the downward/upward flux ratio (Figs.33d,36d) for muon (solid line) and electron (dashed line) neutrinos is different for the two sites. At both experimental sites this ratio is almost the same for neutrinos of different kind with energies less than 1.5 GeV. There is no strong indication that the ratio of muon to electron fluxes at the Soudan site (Fig.32c,33d) depends on the neutrino direction. On the contrary the tendency for the ratio of muon to electron fluxes to decrease for directions close to a horizontal plane is clearly seen for neutrinos in the energy range 0.31 - 1.0 GeV in Gran Sasso (solid line in Fig.36c).

## 8 Summary and conclusions

The proton flux near the Earth consists of primary protons from space, primary protons scattered in the atmosphere and secondary protons produced in atmospheric showers. Each component of the proton flux has different angular distribution. The results of the above Monte Carlo study of propagation in space and interaction with the atmosphere of cosmic particles show that the primary proton flux in the vicinity of the Earth is essentially anisotropic at the directions approaching the horizontal plane. This result is in contradiction with the assumption of most atmospheric neutrino flux calculations. The anisotropy has an important impact on the production of higher energy neutrino.

The momentum spectra of primary protons from space display at different latitudes the expected  behaviour connected with the Earth’s magnetic field (cutoff). There exist, however, a non negligible amount of primary protons with momenta well below the magnetic cutoff.

The momenta of the considerable number of scattered primary protons present at all magnetic latitudes are close to the cutoff value. When backtraced, the scattered protons have ”forbidden” trajectories and consequently their contribution to atmospheric neutrino production was not taken into account in previous calculations of neutrino fluxes.

Predominantly secondary under cutoff particles (protons,$`e^+,e^{}`$) are present in equal amount and the same momenta in the fluxes of particles going to and from the Earth. The simulation correctly describes the fluxes of protons and leptons ($`e^+,e^{}`$) measured by AMS .

The atmospheric neutrinos produced in the simulation have the energy spectra, relative fluxes and angular distributions different from those obtained in the previous atmospheric neutrino flux calculations . The difference is a direct consequence of the anisotropy of the cosmic ray flux at the vicinity of the Earth. As it follows from the present study, the flux of primary cosmic particles at their entrance into the atmosphere is diminishing when approaching the horizontal direction (see Fig.11). Because of that, the fraction of high energy secondaries capable to decay before they lose most of their energy on the ionization of the atmosphere decreases. This results in smaller neutrino fluxes at higher energies.

The up/down asymmetry of the neutrino fluxes resulting from this study is essentially different from the previous calculations. Similar to previous 3-D calculations the flux of atmospheric neutrinos has a maximum in the angular region around the horizontal plane. However, as it is demonstrated in this study there is an important contribution to the horizontal flux from neutrinos produced just under the Earth’s surface. A part of the flux of the ”underground” neutrinos can be detected in the underground experiments. The horizontal region is the most important for determination of the up/down asymmetry of the atmospheric neutrino fluxes and consequently the result of the present study may affect the eventual experimental measurement of the neutrino oscillation parameter $`\mathrm{\Delta }\mathrm{m}^2`$. Calculated for different experimental sites the angular distributions of neutrino fluxes of different kind and their dependence on the neutrino energy suggest the presence of effects faking disappearance of upward going muon neutrinos. These effects are different at different experimental sites (Figures 30c,d; 33c,d,36c,d). This could probably be the reason why different experiments do not agree in their estimate of oscillation parameters.

Acknowlegements. I’m grateful to Prof.Yu.Galaktionov for very useful and stimulating discussions. I have to express my gratitude to Dr.V.Shoutko for reading the manuscript and giving suggestions for important improvements. Finaly I would like to thank Prof.S.C.C.Ting for his interest in this work and encouragement.


text:
# Comment on: necessary and sufficient condition of separability of any system

## A The content of the paper

Deciding whether a given quantum state is entangled or not is a central problem in quantum information theory. For the smallest non-trivial systems (with $`2\times 2`$ resp. $`2\times 3`$ dimensional Hilbert spaces) the positivity of the partial transpose is known to be an efficient necessary and sufficient criterion . Beyond these special cases however, no such calculable criterion is known. For this reason a paper titled necessary and sufficient condition of separability of any system should be of interest for the whole quantum information community (this is one reason, why we post this comment on the quant-ph server).

For the benefit of other readers of the archive, we recapitulate the main result in  and show that it is nothing but a reformulation of the definition of separability — which is naturally a necessary and sufficient criterion for itself.

Before we give a slightly reformulated but textual unchanged version of the main result in  let us introduce the following two-dimensional projectors:

$$P_{ij}=|ii|+|jj|ij.$$
(1)
Such a projection can be understood, as a projection from a larger Hilbert space $`𝐂^n`$ to the qubit space $`𝐂^2`$. Although the results in  are not expressed in terms of these projectors, we will make use of them in order to avoid lengthy and cumbersome notations. Any proposition in  has then a simple and rather intuitive translation. For example, take a vector $`|\mathrm{\Psi }=_{ij}A_{ij}|ij`$. To say that $`\{A_{11},A_{12}\}`$ is parallel to $`\{A_{21},A_{22}\}`$ is now equivalent to the fact that the respective two qubit projection $`P_{12}P_{12}|\mathrm{\Psi }`$ is separable . The main theorem in  then reads as follows:

###### Theorem 1

Let $`\rho `$ be a mixed state on $`𝐂^n𝐂^m`$ with eigenvectors $`\{x_k\}`$ (unnormalized). Let

$$\rho ^{ijkl}=P_{ki}P_{lj}\rho P_{ki}P_{lj}$$

denote the projection of $`\rho `$ to a two qubit space. Then $`\rho `$ is separable iff the following two conditions hold:

1. For all two-qubit projections $`\rho ^{ijkl}`$ the concurrence (denoted as $`a^r`$) as introduced by Wootters  is negative (or zero, if we take the original definition $`c=\mathrm{max}\{0,a^r\}`$) .
2. There exists an isometry $`U`$ corresponding to a decomposition $`\rho =_k|z_kz_k|`$ with $`|z_k=_jU_{kj}|x_j`$, such that every two qubit projection $`P_{ki}P_{lj}|z_k`$ is a product vector .

## B What it is about

Now let us have a closer look at the two conditions in Theorem 1. The first condition utilizes Wootters formula for the concurrence as separability criterion for the mixed two qubit states (even the main part of Wootters proof is repeated). However, this is rather confusing since their final result is completely independent of the choice of the separability criterion at this point. One could equivalently take the partial transposition criterion.

Condition 2 is equivalent to saying that all the vectors $`|z_k`$ have to be separable (see page 2 in ). But then the search for the transformation $`U`$ is nothing but a search for a decomposition into product vectors, i.e., just a rephrasing of the initial problem. Notwithstanding the fact that condition 1 is therefore redundant it is not even a good necessary criterion for separability, since the usual PPT-criterion is both, easier to calculate and stronger.

The necessary and sufficient separability criterion, which is “practical in judging” thus turns out to be a tautology. The fact that it works for an example of a bound entangled state is not surprising, since all the known examples of such states are just constructed in a way that it is easily seen that they are entangled (i.e., they do not have any product vector in their range).


text:
# Numerical computation of resonance poles in scattering theory

## Abstract

We present a possible way of computing resonance poles and modes in scattering theory. Numerical examples are given for the scattering of electromagnetic waves by finite-size photonic crystals.

Resonance poles are the main quantity of interest in scattering theory . They are poles of the scattering matrix, and may be considered as generalized eigenvalues to which generalized eigenmodes are associated. These poles are complex ones, i.e. they correspond to complex values of the energy in the scattering theory of Schrödinger equation and to complex frequencies for the scattering theory of Maxwell system.

When dealing for instance with finite size photonic crystals, i.e. periodic arrangements of scatterers that are finite in at least one direction of space , one cannot use Bloch waves theory to compute the electromagnetic properties of the structure and one has to retreat to the computation of the scattering matrix, that is the operator $`S\left(k\right)`$ such that $`U=S\left(k\right)U^i`$ where $`U^i`$ is the harmonic incident field, with wavevector $`k`$, and $`U`$ is the total field. The scattering matrix writes $`S\left(k\right)=I_d+T\left(k\right)`$ where $`T\left(k\right)`$ is the so-called scattering amplitude. Let us assume that there exists some pole $`k_p`$ of $`T`$ in some neighborhood $`𝒱`$ of the complex plane, then locally the scattering amplitude writes $`T\left(k\right)=\frac{P_p}{kk_p}+T_0\left(k\right)`$ where $`P_p`$ is a residu operator and $`T_0`$ is holomorphic in $`𝒱`$. Operator $`P_p`$ is a finite rank operator and its range is precisely the kernel of $`T^1\left(k_p\right)`$. It is define in an abstract way as the Cauchy integral:

$$P_p=\frac{1}{2i\pi }T\left(z\right)𝑑z$$
(1)
where integration takes place on a loop oriented in the direct sense enclosing the only pole $`k_p`$. Another way of defining the projection operator $`P_p`$ is to define it as the following limit

$$P_p=s\underset{kk_p}{lim}\left(kk_p\right)T(k)$$
(2)
The point of this note is to show that the first abstract definition (1) can be turned into a useful numerical tool for both the computation of the value of the pole and of the residu operator, and hence the generalized eigenmode, whereas the second is useless. From a numerical point of view, we of course only deal with finite rank operators and the scattering amplitude admits a representation as an operator on $`\mathrm{}^2\left(\right)`$, that is as a matrix, in the usual meaning, acting on double complex sequences. Once this representation is given, the residu operator can be computed provided that a region of the complex plane containing only one pole can be precised. This means that it suffices to know the value of the pole with a very poor precision to be able to compute the residu operator, which is not the case when the second definition (3) is used: in that last case numerical instabilities necessarily occur as it uses the product of a singular matrix by a term tending to zero, which is a very bad numerical situation. From a practical point of view, one has to define a path $`\gamma :t[0,1]\gamma \left(t\right)`$ whose graph is a loop enclosing $`k_p`$ and to compute numerically the integral $`_0^1T\left(\gamma \left(t\right)\right)\gamma ^{}\left(t\right)𝑑t`$ for which any reasonable numerical method works. However a precise computation of the pole is useful when one wishes to compute a map of the electromagnetic field of the pole, for in that case a particular basis such as Hankel-Fourier series are used, i.e. the field is expanded on the basis $`\left(H_n^{\left(1\right)}\left(k_pr\right)\mathrm{exp}\left(in\theta \right)\right)_n`$ . A possible way is to use a Müller like algorithm  and to compute a zero of the determinant of $`T^1\left(k\right)`$. However this matrix is generally badly conditionned and a better idea is to compute the smallest eigenvalue of $`T^1\left(k\right)`$. This works well in case of a finite size crystal, but this is not always the case: for instance, when modelizing photonic crystals by stacks of gratings  and introducing periodic defects, convergence problems may occur when using Müller algorithm . We suggest then to compute the following Cauchy integral:

$$\frac{1}{2i\pi }zT\left(z\right)𝑑z=k_pP_p.$$
(3)
Recalling that $`P_p`$ has finite rank and hence has only a finite number of eigenvalues, a simple comparison of this last integral with $`P_p`$ gives the value of $`k_p`$ with a very good accuracy. Of course formula (3) only holds when $`k_p`$ is a pole with multiplicity $`1`$: this case is a very common one. One should not mistaken the range of $`P_p`$ and the multiplicity of $`k_p`$: it is possible that the multiplicity of $`k_p`$ is $`1`$ while the rank of $`P_p`$ is greater than $`1`$ .

Let us now turn to some numerical applications. We deal with the structure depicted in figure 1. It is a collection of $`7\times 7`$ homogeneous fibers with relative permittivity $`\epsilon =9`$, the radius of the rods $`R=1/2`$ and the spacing is $`d=1`$ (these values are given in arbitary units). We use a rigorous modal theory of diffraction to compute the scattering matrix of this system . All the numerical results have been obtained using a standard PC computer. Removing a rod at the center of the crystal, a defect mode appears in the gap  (see fig. 2 for the transmission spectrum). To this peak in the transmission spectrum corresponds a pole $`k_p`$. The reference value that we use for convergence comparison is $`k_p=2.329197035861340.00378267987614i.`$ This value has been computed using Müller algorithm by minimising the smallest eigenvalue of $`T^1`$. For this given value of $`k_p`$ the smallest eigenvalue has modulus inferior to $`10^{14}`$. We then compute both Cauchy integrals (1,3). We use an integration path that is a triangle whose vertices have affixes: $`(2.3,2.40.1i,2.4)`$. We use the integration algorithm described in  and we denote by $`k_{N,p}`$ the numerical value obtained by using $`N`$ points of integration, which we compare with the above value $`k_p`$ which is the best numerical value that we can obtain . A very good precision is rapidly obtained (see fig. 3): for instance, using a discretisation of $`15`$ points we obtain $`6`$ exact figures though with such a rough discretization we only get operator $`P_p`$ with a low precision. In fact, it seems that the proportional coefficient between both integrals is not much affected by the precision with which $`P_p`$ is computed. A finer computation of integral (1) gives the defect mode. The convergence can be checked by looking at the non-zero eigenvalue of $`P_p`$ (fig. 4), here the reference eigenvalue (i.e. the best numerical value for a precision error of $`10^{15}`$) is obtained with $`N=150`$. A much finer discretization than in the case of the pole is required to get a good representation of the defect mode, though the computing time is perfectly accessible with a very basic PC.

In conclusion, we have shown that it was possible to turn a rather abstract mathematical object into a useful numerical tool. This technique applies as well for any situation in which a meromorphic operator with non essential poles is involved, which is the usual case.

Figures captions

Figure 1: Sketch of the 2D photonic crystal. The transmission ratio is computed as the flux of the Poynting vector through the segment indicated below the crystal.

Figure 2: Transmission ratio versus the wavenumber for an incident plane wave.

Figure 3: Convergence of the value of the pole versus the number of integration points.

Figure 4: Convergence of the eigenvalue of the projection operator versus the number of integration points.


text:
# Finite market size as a source of extreme wealth inequality and market instability

## Acknowledgements

We thank Dietrich Stauffer for very helpful discussions and comments. Z.F.H. acknowledges the financial support of SFB 341.