format: text

text:
# Andreev Bound States and Self-Consistent Gap Functions for SNS and SNSNS Systems

## I Introduction

In 1991, Tanaka and Tsukada  investigated the energy spectrum of the quasiparticle states in a superconducting superlattice based on a Kronig-Penney model potential. The calculations were done using the Andreev approximation  and the system considered was infinite in both the transverse and the longitudinal dimensions. In this paper we concentrate on three additional aspects: i) we do the exact calculation and investigate the reliability of the Andreev approximation, ii) we do the calculation for a limited number of layers, which relates better to a possible experimental situation, and iii) we investigate the transverse size dependence of the relevant properties such as the local density of states. The systems we consider cover the entire range from very narrow transverse dimensions to wider ones, thereby simulating 1, 2 and 3-dimensional systems.

The authors are aware of two studies focusing on the transverse size dependence of Andreev bound states and the quasiparticle local density of states. Šipr and Györffy  show the splitting of the degenerate bound states found in the Andreev approximation by studying the exact solution of the Bogoliubov-de Gennes equations. Their study is limited to a superconductor-normal metal-superconductor (SNS) junction and no self-consistency of the gap is considered. A second study, by Blaauboer et al. , is a rather global one, but in that study the multiple-scattering Green’s function formalism developed by Koperdraad  was applied for the first time. The local density of states of normal metal-superconductor (NS) and SNS junctions were calculated, but the breakdown of the Andreev approximation was not noticed.

We want to extend these studies by considering self-consistently obtained gap functions and by studying a system with more than two interfaces, looking for the influence of an additional superconducting layer to the Andreev bound state spectrum. The multiple-scattering Green’s function formalism to be used  is an extension of the formalism introduced by Tanaka and Tsukada  in that no Andreev approximation is made. In addition it is set up in a quite different way, such that the multiple scattering of the (quasi)particles by the different interfaces is exhibited explicitly.

The paper is organized as follows. First we give a concise account of the main features of the theory. In Sec. III the local density of states for the two systems considered is studied. In Sec. IV the gap function is calculated self-consistently for a bar-shaped superconductor and for the NS, SNS and SNSNS junctions. Concluding remarks are given in Sec. V.

Throughout the paper, Rydberg atomic units are used, in which the energy is in Rydberg, the distance is in Bohr ($`1`$ Bohr $``$ 0.5 Å), $`\mathrm{}=1`$, and the electronic mass is $`\frac{1}{2}`$.

## II Theory

Although in this paper a Green’s function formalism is used, we first give the Bogoliubov-de Gennes equations

$`\left[\begin{array}{cc}^2\mu & \mathrm{\Delta }(𝐫)\\ \mathrm{\Delta }^{}(𝐫)& ^2+\mu \end{array}\right]\mathrm{\Psi }(𝐫)=E\mathrm{\Psi }(𝐫)E\left(\begin{array}{c}u(𝐫)\\ v(𝐫)\end{array}\right)`$ (5)
used by Šipr and Györffy  to investigate the bound states of a SNS junction shown in the upper panels of FIG. 2. The main reason is that our formalism will be expressed in terms of the solutions of these equations. The inhomogeneity of the system is expressed by the space dependence of the gap $`\mathrm{\Delta }(𝐫)`$. In the superconducting parts the gap is a complex constant whereas in the normal metal it is zero. The spinor wavefunction describes quasiparticle excitations and the energy $`E`$ is measured with respect to the Fermi energy $`\mu `$.

Application of the Bogoliubov-de Gennes equations to the bar-shaped superconductor shown in FIG. 1 yields

$$\mathrm{\Psi }(𝐫)=\left(\begin{array}{c}u_S^\sigma e^{i\varphi /2}\\ u_S^\sigma e^{i\varphi /2}\end{array}\right)e^{i\sigma \nu k_S^\sigma x}\mathrm{sin}\left(\frac{n_y\pi }{L_y}y\right)\mathrm{sin}\left(\frac{n_z\pi }{L_z}z\right),$$
(6)
where $`\varphi `$ is the phase of the complex constant $`\mathrm{\Delta }`$, and

$`k_y`$ $`=`$ $`{\displaystyle \frac{n_y\pi }{L_y}};k_z={\displaystyle \frac{n_z\pi }{L_z}}\mathrm{with}n_y,n_z=1,2,3,\mathrm{}`$ (7)
$`u_S^\sigma `$ $`=`$ $`\sqrt{E+\sigma \sqrt{E^2|\mathrm{\Delta }|^2}}`$ (8)
$`k_S^\sigma `$ $`=`$ $`\sqrt{k_{F_x}^2+\sigma \sqrt{E^2|\mathrm{\Delta }|^2}}`$ (9)
$`k_{F_x}^2`$ $`=`$ $`\mu k_y^2k_z^2`$ (10)
The discretized nature of the transverse wavevectors is due to the vanishing of the wavefunctions on the transverse boundaries. The four standard solutions are labeled with the sign indices $`\sigma `$ and $`\nu `$ that can both be equal to $`\pm 1`$. By this convention, the index $`\sigma `$ refers to the type of particle (electronlike for $`\sigma =+1`$ and holelike for $`\sigma =1`$) and the index $`\nu `$ indicates the direction of propagation ($`\nu =+1`$ to the right and $`\nu =1`$ to the left). The complete solution of Eq. (5) is a linear superposition of Eq. (6) for different possible combinations of $`(\sigma ,\nu )`$. The above equations also hold for a normal metal by letting $`\mathrm{\Delta }=0`$, in which case the subscript $`S`$ is to be replaced by $`N`$.

The Green’s function formalism is outlined extensively by Koperdraad et al. . It is an extension of the microscopic theory used by Tanaka and Tsukada, in that the electron-hole scattering properties are treated exactly, and it traces back to the microscopic description of superconductivity by Gor’kov  and Ishii . For the sake of completeness and clarity we summarize its main features and add relevant new elaborations. The matrix Green’s function satisfies the equation

$$\left[\begin{array}{cc}i\omega _n+^2+\mu & \mathrm{\Delta }(𝐫)\\ \mathrm{\Delta }^{}(𝐫)& i\omega _n^2\mu \end{array}\right]G(𝐫,𝐫^{},i\omega _n)=\delta (𝐫𝐫^{})\mathrm{𝟏}$$
(11)
in which the differential operator is closely related to the operator in the Bogoliubov-de Gennes equations (5) apart from the replacement of $`E`$ by $`i\omega _n`$ where

$$\omega _n=\pi nk_BT\mathrm{with}n\mathrm{a}\pm \mathrm{odd}\mathrm{integer}.$$
(12)
The quantity $`\omega _n`$ is called the Matsubara frequency. Possible inhomogeneities of the system are fully represented by the $`𝐫`$ dependence of the gap function. As far as the transverse directions are concerned the general solution of Eq. (11) can be expressed as a linear combination of solutions (6) over all allowed values of $`k_y,k_y^{},k_z,`$ and $`k_z^{}`$. Thus

$$G(𝐫,𝐫^{},i\omega _n)=\frac{4}{L_yL_z}\underset{k_y,k_y^{},k_z,k_z^{}}{}G(x,x^{},k_y,k_y^{},k_z,k_z^{},i\omega _n)\mathrm{sin}k_yy\mathrm{sin}k_y^{}y^{}\mathrm{sin}k_zz\mathrm{sin}k_z^{}z^{}.$$
(13)
The transverse standing waves are normalized to unity, which implies the orthogonality condition

$$_0^{L_y}\mathrm{sin}k_yy\mathrm{sin}\overline{k}_yydy=\frac{L_y}{2}\delta _{k_y,\overline{k}_y},$$
(14)
for the y coordinate and a similar condition for the z coordinate. The Fourier coefficient $`G(x,x^{},k_y,k_y^{},k_z,k_z^{},i\omega _n)`$ can be obtained from Eq. (13) using Eq. (14). It takes the form,

$`G(x,x^{},k_y,k_y^{},k_z,k_z^{},i\omega _n)`$ $`=`$ $`{\displaystyle \frac{4}{L_yL_z}}{\displaystyle _0^{L_y}}\mathrm{sin}k_yy{\displaystyle _0^{L_y}}\mathrm{sin}k_y^{}y^{}{\displaystyle _0^{L_z}}\mathrm{sin}k_zz{\displaystyle _0^{L_z}}\mathrm{sin}k_z^{}z^{}`$ (16)
$`\times G(𝐫,𝐫^{},i\omega _n)dz^{}dzdy^{}dy.`$
Using Eq. (16), Eq. (11) can be written as

$$\left[\begin{array}{cc}i\omega _n+\frac{d^2}{dx^2}+k_{F_x}^2& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& i\omega _n\frac{d^2}{dx^2}k_{F_x}^2\end{array}\right]G(x,x^{},k_y,k_y^{},k_z,k_z^{},i\omega _n)=\delta (xx^{})\delta _{k_y,k_y^{}}\delta _{k_z,k_z^{}}\mathrm{𝟏}$$
(17)
Before we proceed we want to make the terminology clear. $`G(𝐫,𝐫^{},i\omega _n)`$ is the actual Green’s function. The Fourier coefficient $`G(x,x^{},k_y,k_y^{},k_z,k_z^{},i\omega _n)`$ is a Green’s function in the sense that it is the solution of Eq. (17), but this equation demonstrates that it is diagonal in $`k_y`$ and $`k_z`$. That facilitates both the notation, because the variables $`k_y^{}`$ and $`k_z^{}`$ can be omitted, and the calculation of $`G(𝐫,𝐫^{},i\omega _n)`$ according to Eq. (13). In calculating the local density of states and the self-consistent gap function, we will need the Green’s function for diagonal spatial coordinates. As long as we keep $`xx^{}`$, we can already take $`y=y^{}`$ and $`z=z^{}`$ in Eq. (13). Since we are only interested in the variations over the longitudinal direction $`x`$, we can average over the transverse dimensions. By that Eq. (13) simplifies to the series

$$G(x,x^{},i\omega _n)=\frac{1}{L_yL_z}\underset{k_y,k_z}{}G(x,x^{},k_y,k_z,i\omega _n).$$
(18)
First we give the Green’s function of a homogeneous bar-shaped superconductor

$$G_S^0(x,x^{},k_y,k_z,i\omega _n)=\underset{\sigma }{}d_S^\sigma \psi _S^\sigma (x_<)\stackrel{~}{\psi }_S^{\sigma +}(x_>)=\underset{\sigma }{}d_S^\sigma \psi _S^{\sigma +}(x_>)\stackrel{~}{\psi }_S^\sigma (x_<),\mathrm{with}d_S^\sigma =\frac{1}{4\mathrm{\Omega }_Sk_S^\sigma },$$
(19)
in which $`i\mathrm{\Omega }_S=\sqrt{(i\omega _n)^2|\mathrm{\Delta }|^2}`$. The wavefunction $`\psi _S^{\sigma \nu }\left(x\right)`$ is given by Eq. (6) after omitting the transverse solutions, and $`u_S^\sigma `$ and $`k_S^\sigma `$ are given by Eqs. (8) and (9) with the replacement $`Ei\omega _n`$. In addition the conjugate wavefunction is defined by

$$\stackrel{~}{\psi }_S^{\sigma \nu }(x)=\left(\begin{array}{cc}u_S^\sigma e^{i\varphi /2}& u_S^\sigma e^{i\varphi /2}\end{array}\right)e^{i\sigma \nu k_S^\sigma x}$$
(20)
The Green’s function (19) describes the propagation of excitations (holelike or electronlike) from the starting point $`x^{}`$ to the final point $`x`$ with the weighting factor $`d_S^\sigma `$. In this case of a superconducting bar, there is no scattering.

Now we give the Green’s function for a system with one interface. To account for the scattering at the interface, its appropriate form appears to be

$$G_{\nu j\nu ^{}j}(x,x^{},k_y,k_z,i\omega _n)=G_{\nu j}^0(x,x^{},k_y,k_z,i\omega _n)\delta _{\nu \nu ^{}}+\underset{\sigma \sigma ^{}}{}d_{\nu j}^\sigma d_{\nu ^{}j}^\sigma ^{}\psi _{\nu j}^{\sigma \nu }(x)t_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}\stackrel{~}{\psi }_{\nu ^{}j}^{\sigma ^{}\nu ^{}}(x^{})$$
(21)
where

$$G_{\nu j}^0(x,x^{},k_y,k_z,i\omega _n)=\underset{\sigma }{}d_{\nu j}^\sigma \psi _{\nu j}^{\sigma \mu }(x)\stackrel{~}{\psi }_{\nu ^{}j}^{\sigma ^{},\mu }(x^{})\mathrm{with}\mu =\mathrm{sgn}(xx^{}).$$
(22)
The label $`S`$ in Eq. (19) has been replaced with the more flexible label $`\nu j`$ indicating the position $`x_j`$ of the interface. The subscript $`+j`$ means $`x`$ is in the right-hand side of the interface $`x_j`$ and $`j`$ means that it is in the left-hand side. The first term is nonzero only if the starting and final points, $`x^{}`$ and $`x`$ respectively, are at the same side of the interface. The second term takes into account the scattering of the particle at the interface. This scattering is aptly described by the parameter we call scattering $`t`$-matrix, $`t_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}`$. For a superconducting bar the $`t`$-matrix is zero, because there is pure propagation and no scattering. The scattering $`t`$-matrix can be obtained by imposing the continuity of the Green’s function and its derivative at the interface. We then obtain

$$\underset{\sigma \nu }{}\nu d_{\nu j}^\sigma \psi _{\nu j}^{\sigma \nu }(x_j)t_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}=\nu ^{}\psi _{\nu ^{}j}^{\sigma ^{}\nu ^{}}(x_j).$$
(23)
In applying the above interface matching condition, it has come out to be convenient to use an extended definition

$$\psi _S^{\sigma \nu }(x)=\left(\begin{array}{c}u_S^\sigma e^{i\varphi /2}\\ u_S^\sigma e^{i\varphi /2}\\ i\sigma \nu k_S^\sigma u_S^\sigma e^{i\varphi /2}\\ i\sigma \nu k_S^\sigma u_S^\sigma e^{i\varphi /2}\end{array}\right)e^{i\sigma \nu k_S^\sigma x}$$
(24)
of the wavefunction $`\psi _S^{\sigma \nu }\left(x\right)`$ which includes the derivative with respect to $`x`$.

If there are more interfaces as in the systems depicted in FIG. 2, multiple scattering occurs and the Green’s function is given by

$`G_{\nu j\nu {}_{}{}^{}j_{}^{}}(x,x^{},k_y,k_z,i\omega _n)`$ $`=`$ $`G_{\nu j\nu {}_{}{}^{}j_{}^{}}^0(x,x^{},k_y,k_z,i\omega _n)(\delta _{\nu \nu ^{}}\delta _{jj^{}}+\delta _{\nu \nu ^{}}\delta _{j+\nu ,j^{}})`$ (26)
$`+{\displaystyle \underset{\sigma \mu \sigma {}_{}{}^{}\mu ^{}}{}}d_{\nu j}^\sigma d_{\nu {}_{}{}^{}j_{}^{}}^\sigma ^{}\psi _{\nu j}^{\sigma \mu }(x)T_{\nu j\nu {}_{}{}^{}j_{}^{}}^{\sigma \sigma {}_{}{}^{}\mu \mu ^{}}\stackrel{~}{\psi }_{\nu {}_{}{}^{}j_{}^{}}^{\sigma {}_{}{}^{}\mu ^{}}(x^{}).`$
This is the generalization of Eq. (21). There are features in Eq. (26) that do not appear in that equation. One of these features is the strange combination of Kronecker deltas in the first term. The first set of Kronecker deltas, $`\delta _{\nu \nu ^{}}\delta _{jj^{}}`$, insures that the first term is nonzero only if the starting and final positions are in the same layer. The other set, $`\delta _{\nu \nu ^{}}\delta _{j+\nu ,j^{}}`$, serves the same purpose but at the same time it takes care of the redundancy in the indexing of the layers. The structure of the two sets of Kronecker deltas assures us that there is no overlapping of their functions. If one set gives the value unity the other set must be zero and vice versa. Another feature of Eq. (26) is the presence of the quantity $`T_{\nu j\nu ^{}j^{}}^{\sigma \sigma ^{}\mu \mu ^{}}`$. Whereas the $`t`$-matrix describes the scattering of the particle at a particular interface, this quantity describes the multiple scattering of the particle at the interfaces along its path. We call this the multiple scattering $`T`$-matrix. The $`T`$-matrix is in fact a function of the $`t`$-matrix and is given by the multiple scattering equation

$$T_{\nu j\nu ^{}j^{}}^{\sigma \sigma ^{}\nu \mu ^{}}=t_{\nu j\mu {}_{}{}^{}j_{}^{}}^{\sigma \sigma {}_{}{}^{}\nu \mu ^{}}(\delta _{\mu {}_{}{}^{}\nu ^{}}\delta _{jj^{}}+\delta _{\mu ^{}\nu ^{}}\delta _{j+\mu ^{},j^{}})+\underset{\sigma {}_{}{}^{\prime \prime }\nu ^{\prime \prime }}{}t_{\nu j\nu {}_{}{}^{\prime \prime }j}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}d_{\nu {}_{}{}^{\prime \prime }j}^{\sigma ^{\prime \prime }}T_{\nu {}_{}{}^{\prime \prime }j\nu {}_{}{}^{}j_{}^{}}^{\sigma {}_{}{}^{\prime \prime }\sigma {}_{}{}^{},\nu {}_{}{}^{\prime \prime }\mu ^{}}.$$
(27)
To solve this equation it is necessary to first solve for the $`t`$-matrix at each interface using Eq. (23).

Until now no approximations have been made in using the solutions (6) of the Bogoliubov-de Gennes equations. In applying the formalism given above to the calculation of the density of states, $`i\omega _n`$ has to be replaced by $`E+i\delta `$. After that one can make the frequently used Andreev approximation, which amounts to the replacement

$$k_{\nu j}^\sigma k_{F_x}+\sigma \frac{\sqrt{E^2|\mathrm{\Delta }|^2}}{2k_{F_x}}$$
(28)
if $`k_{\nu j}^\sigma `$ occurs in the exponential and to $`k_{\nu j}^\sigma k_{F_x}`$ if $`k_{\nu j}^\sigma `$ occurs as a factor as shown explicitly in the third and fourth components of the wavefunction (24). This approximation is valid when $`E,|\mathrm{\Delta }|<<k_{F_x}^2`$. In the present paper we will investigate its limitation by looking at configurations in which $`E,\mathrm{\Delta }k_{F_x}^2`$.

## III The Local Density of States

Sec. II provides the basic machinery to determine the matrix Green’s function, which enables us to calculate the local density of states (LDOS) at a position $`x`$ using the equation

$$\mathrm{LDOS}(x,E)=\frac{1}{\pi L_yL_z}\underset{\delta 0}{lim}\underset{k_y,k_z}{}\mathrm{Im}G_{11}(x,x;k_y,k_z,E+i\delta )$$
(29)
in which $`G_{11}`$ is the upper left matrix element of the multiple scattering Green’s function (26). Although we only study Andreev bound states, which imply infinite peaks at the bound-state energies, the Green’s function formalism makes it possible to broaden the peaks by using a small value of $`\delta `$. The peaks acquire a finite height and must correspond to the bound-state energies.

### A The SNS Junction

We first apply the formalism discussed in Sec. II to a superconductor-normal metal-superconductor (SNS) junction. A schematic diagram of the system studied is given in the upper and middle panels of FIG. 2. The dotted lines in the upper panel serve as an indication of the position of the inner NS and SN interfaces in the SNSNS system to be treated in Sec. III B. In this first application we will show more explicitly how the different labels are used. The interface index $`j`$ has only two values $`1`$ and $`2`$. The $`T`$-matrix equation (27) can be recast into the form

$$T_{\nu j\nu {}_{}{}^{}j_{}^{}}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}=t_{\nu j\nu {}_{}{}^{}j_{}^{}}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}\delta _{jj^{}}+\underset{\sigma {}_{}{}^{\prime \prime }\nu ^{\prime \prime }}{}t_{\nu j\nu {}_{}{}^{\prime \prime }j_{}^{}}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}d_{\nu {}_{}{}^{\prime \prime }j}^{\sigma ^{\prime \prime }}T_{\nu {}_{}{}^{\prime \prime },j+\nu {}_{}{}^{\prime \prime },\nu {}_{}{}^{}j_{}^{}}^{\sigma {}_{}{}^{\prime \prime }\sigma {}_{}{}^{},\nu {}_{}{}^{\prime \prime }\nu ^{}}.$$
(30)
To implement this in matrix form, we must see to it that the elements of the $`T`$-matrix in both sides of the equation are arranged in the same manner. Thus, its elements in the right-hand side of the equation must be similarly displayed as in the left-hand side. This can be done by writing Eq. (30) in the form

$`\left[\begin{array}{cc}T_{\nu 1\nu {}_{}{}^{}1}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}& T_{\nu 1\nu {}_{}{}^{}2}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}\\ T_{\nu 2\nu {}_{}{}^{}1}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}& T_{\nu 2\nu {}_{}{}^{}2}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}\end{array}\right]=\left[\begin{array}{cc}t_{\nu 1\nu {}_{}{}^{}1}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}& 0\\ 0& t_{\nu 2\nu {}_{}{}^{}2}^{\sigma \sigma {}_{}{}^{}\nu \nu ^{}}\end{array}\right]+`$

$$\left[\begin{array}{cc}0& t_{\nu 1\nu {}_{}{}^{\prime \prime }1}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}d_{\nu {}_{}{}^{\prime \prime }1}^{\sigma ^{\prime \prime }}\delta _{\nu {}_{}{}^{\prime \prime }+}\\ t_{\nu 2\nu {}_{}{}^{\prime \prime }2}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}d_{\nu {}_{}{}^{\prime \prime }2}^{\sigma ^{\prime \prime }}\delta _{\nu {}_{}{}^{\prime \prime }}& 0\end{array}\right]\left[\begin{array}{cc}T_{\nu 1\nu {}_{}{}^{\prime \prime }1}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}& T_{\nu 1\nu {}_{}{}^{\prime \prime }2}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}\\ T_{\nu 2\nu {}_{}{}^{\prime \prime }1}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}& T_{\nu 2\nu {}_{}{}^{\prime \prime }2}^{\sigma \sigma {}_{}{}^{\prime \prime }\nu \nu ^{\prime \prime }}\end{array}\right].$$
(31)
Each element of the above matrices is itself a $`4\times 4`$ matrix with $`\sigma \nu `$ $`(=++,+,+,)`$ as row and column indices. So actually, the matrices appearing in Eq. (31) have dimensions $`8\times 8`$. In this form the $`T`$-matrix can be obtained by appropriate matrix inversion. The Green’s function is obtained using Eq. (26). The local density of states is determined by using Eq. (29).

To investigate the validity of the Andreev approximation in our SNS junction, we will focus on the choice of the transverse width $`L_y=L_z=L_t`$ of the junctions. The transverse wave components of the wavefunction (6) are standing waves proportional to $`\mathrm{sin}(k_yy)\mathrm{sin}(k_zz)`$ where $`k_y`$ and $`k_z`$ are given by Eq. (7). The different combinations of $`(k_y,k_z)`$ or $`(n_y,n_z)`$ are called modes whose allowed values are determined by

$$k_{F_x}^2=\mu \left(\frac{\pi }{L_t}\right)^2(n_y^2+n_z^2)>0.$$
(32)
When the transverse dimension is small, the second term in the right becomes large, as a result, only a few modes will be allowed. If this term exceeds the chemical potential $`\mu `$, $`k_{F_x}`$ becomes imaginary, the wavefunction is damped, and consequently, such mode cannot exist. For larger transverse dimensions, the second term is smaller whereupon more modes are allowed. Most of our calculations will be done for small $`L_t`$ so that only few modes will exist. We will tune $`L_t`$ such that $`k_{F_x}^2`$ is of the same order of magnitude as the gap energy $`\mathrm{\Delta }`$, in which regime the Andreev approximation (28) is not valid, and call such a $`L_t`$ value a critical width.

FIGs. 3 and 4 show the results for a configuration in which $`(n_y,n_z)=(2,2)`$ is the highest allowed mode. The chemical potentials in the superconductor and in the normal metal, $`\mu _S`$ and $`\mu _N`$ respectively, are assumed equal with magnitude 0.5. The longitudinal dimension $`L`$ of the normal-metal part is $`4000`$ Bohr and the gap $`\mathrm{\Delta }`$ is treated as real with magnitude $`0.0001`$ Ry. The LDOS in the normal-metal part at $`x=1000`$ Bohr is plotted against $`E/\mathrm{\Delta }`$ . The peaks represent discrete energy states . We make the width of the curves, determined by the the parameter $`\delta `$ in $`E+i\delta `$, wide enough so that the fundamental features can be seen. The numbers in parentheses denote the mode to which the energy belongs. In FIG. 4 the transverse width is determined by the condition that $`k_{F_x}^2=\mathrm{\Delta }`$ for the mode $`(2,2)`$, in which one finds that $`L_t=12.5676`$ Bohr and in FIG. 3 the transverse width is $`L_t=13`$ Bohr, which is slightly larger than the critical width but has the same allowed modes. In both figures the dashed curve is the local density of states calculated using Andreev approximation (AA) and the solid curves are the results from the exact calculation. In FIG. 3, the exact and the AA results coincide and just three states are found, one for each mode. For the critical width, the states for the first two modes are at almost the same position, but for the $`(2,2)`$ mode many states are found. We notice that the AA peaks have about the same height whereas the amplitudes of the exact peaks vary widely with the energy. The amplitude variation of the exact peaks is due to the specific position chosen for the LDOS. FIG. 5 depicts the exact LDOS at $`x=1000`$ Bohr and $`x=1500`$ Bohr. We observe that some of the peaks in the solid curve are suppressed while the corresponding peaks in the dashed curve are prominent and vice versa. The suppressed peaks are pulled down by the small magnitude of the weighting wavefunctions at those values of $`x`$. The degenerate traveling waves corresponding to the AA are split in the exact treatment into odd and even (sine and cosine) functions having different heights at different positions .

Finally, we want to comment on our choice of the transverse width up to now. In Sec. IV it will come out that superconductivity is suppressed for transverse widths in the order of $`20`$ Bohr or less. This means that so far our choice of transverse widths may seem not too appropriate. We chose those transverse widths to illustrate with clarity the fundamental features of the bound states. In order to see the influence of a wider transverse dimension, we show in FIG. 6 the LDOS of a SNS system with $`L_t=100`$ Bohr (solid curve) and $`L_t=99.8514`$ Bohr (dotted curve). The critical transverse width of $`99.8514`$ Bohr is obtained from the condition $`k_{F_x}^2=\mathrm{\Delta }`$ for the mode $`(19,12)`$. To illustrate its main features clearly, we only show the results for the Andreev approximation. It appears that higher modes are allowed for a wider transverse dimension. The peak at about $`E=0.95\mathrm{\Delta }`$ comes from the many lower modes, each supporting just one bound state. Only the highest mode gives rise to the distinct set of peaks starting at about $`E=0.074\mathrm{\Delta }`$. So effectively, the character of the results shown above is essentially unchanged.

### B The SNSNS System

A schematic representation of the system is shown in the lower panel of FIG. 2. The height of the middle superconductor is chosen to take the values $`h\mathrm{\Delta }`$ where the range of $`h`$ is $`0h1`$. The LDOS is calculated using Eq. (29). The process of determining the $`t`$-matrix, the $`T`$-matrix, and the Green’s function for the SNSNS junction is the same as in the SNS junction. The number of interfaces has increased by $`2`$, and the dimension of the matrices has become $`16`$ instead of $`8`$. We again choose a square transverse cross section whose dimension is $`L_t=12.5676`$ Bohr. This means that the set of modes is the same as in the SNS system. The lengths of the normal metal parts are each $`1000`$ Bohr and that of the middle superconductor is $`2000`$ Bohr.

FIG. 7 shows the development of the bound states as the gap of the middle superconductor is increased from $`0.25\mathrm{\Delta }`$ to $`\mathrm{\Delta }`$. All peaks belong to the highest $`(2,2)`$ mode, apart from a $`(1,1)`$ mode and a $`(1,2)`$ mode peaks just below $`E/\mathrm{\Delta }=1`$. The mode labels are shown only in FIG. 7d to avoid cluttering the figures. Both the exact solution and the states according to the Andreev approximation (AA) are shown. Again the lifting of degeneracy is observed in the exact results. The distribution of the peak heights shown corresponds to a calculation of the LDOS at $`x=1500`$ Bohr. Results for other positions just give other distributions of peak heights. From now on, we concentrate on the position of the AA peaks to facilitate the comparison of the different figures. The SNSNS system with $`h=0`$ is equivalent to the SNS system, so it is interesting to compare FIG. 4 with FIGs. 7a-d. One can see a general shift of the states in the latter figure to the right relative to the states in FIG. 4. For the system with $`h=0.25`$ the $`h=0`$ state at $`E=0.075\mathrm{\Delta }`$ suffers a large displacement. This can be understood as follows. States below $`E=0.25\mathrm{\Delta }`$ see a longitudinal length $`L_N`$ of $`1000`$ Bohr whereas the states above $`E=0.25\mathrm{\Delta }`$ see a length $`L_N`$ of $`4000`$ Bohr, $`L_N`$ being the length of the N metal. According to a semiclassical result , bound-state energies scale as $`\frac{v_F}{L_N}`$, so they are inversely proportional to the longitudinal length, which would imply a shift upwards to $`E=0.3\mathrm{\Delta }`$. However, the separation of $`2000`$ Bohr between the N part is smaller than the BCS coherence length $`\frac{v_F}{\pi \mathrm{\Delta }}4500`$ Bohr, which suggests that the N parts of the system are not completely decoupled yet. This is the reason why the lowest state is found below $`E=0.25\mathrm{\Delta }`$, namely at $`E=0.19\mathrm{\Delta }`$. The positions of the other peaks are hardly changed.

Looking at FIGs. 7b-d one sees that the lowest state stabilizes at a position of about $`E=0.26\mathrm{\Delta }`$. All other states are shifted more and more upwards for increasing $`h`$ values. For $`h=0.5`$, FIG. 7b, still a set of peaks is found above $`0.5\mathrm{\Delta }`$, in FIG. 7c the set starts at about $`E=0.75\mathrm{\Delta }`$. For $`h=1`$, FIG. 7d, only three states are seen just below $`E=\mathrm{\Delta }`$, apart from a state at $`E=0.78\mathrm{\Delta }`$, lying three times as high as the lowest state, in agreement with the semiclassical picture. A test calculation for a SNS system with $`L_N=1000`$ Bohr differs from FIG. 7d only as far as the position of the highest peak, the $`(1,1)`$ peak, is concerned.

The idea of decoupling can also be illustrated by FIG. 8. For the sake of clarity we just give the results for the Andreev approximation. FIG. 8(a) shows the LDOS for a SNSNS system with $`L=6000`$ Bohr. The lengths of the N parts are kept at $`1000`$ Bohr but now the length of the middle superconductor is $`4000`$ Bohr, slightly below the BCS coherence length of approximately $`4500`$ Bohr. In FIG. 8(b), the length $`L`$ of the SNSNS system is $`8000`$ Bohr and the lengths of the N parts are again kept at $`1000`$ Bohr. This makes the middle superconductor $`6000`$ Bohr in length, which is longer than the BCS coherence length. In both figures the gap of the middle superconductor is $`0.25\mathrm{\Delta }`$. Thus, there is effectively complete decoupling in FIG 8(b). In addition, one sees that the $`(2,2)`$ peaks are lying closer to one another for the longer system, which is in line with the $`\frac{V_F}{L_N}`$ behavior of the energies.

## IV Self-Consistent Calculation of the Gap

In the preceding discussions, the gap function $`\mathrm{\Delta }`$ used in the calculation is steplike, that is, it has a finite constant value in the superconducting part and is zero in the normal part. At the interfaces, it has a step discontinuity. This configuration is shown in FIG. 2. So in the calculation of the bound states in Sec. III, the proximity effect is not taken into account. However we want to demonstrate in this section that our formalism makes it possible to show that the actual gap function is not steplike, but decreases smoothly towards the interface and abruptly goes down to zero in the normal metallic layer. This proximity effect is studied in the present section. In addition, the actual gap function will be calculated self-consistently.

The self-consistency condition is given by

$$\mathrm{\Delta }(x)=VF(𝐫,𝐫,0^+)$$
(33)
in which $`F(𝐫,𝐫,0^+)`$ is the well-known anomalous Green’s function  given by the upper right element of the original $`(𝐫,\tau )`$-dependent matrix Green’s function. Taking into account the expansion of the matrix Green’s function over the Matsubara frequencies and the transverse wave vectors, we obtain

$$\mathrm{\Delta }(x)=\frac{V}{\beta L_yL_z}\underset{\omega _n,k_y,k_z}{}F(x,x,k_y,k_z,i\omega _n)$$
(34)
where $`F(x,x,k_y,k_z,i\omega _n)`$ is the upper right element of the matrix Green’s function in Eq. (26). The summation over the Matsubara frequencies is restricted by the Debye temperature $`\theta _D`$ according to the formula,

$$k_B\theta _D=\omega _{n_{max}}=n_{max}\pi k_BT.$$
(35)
The summation over the transverse wave vector is over all positive values of $`k_y`$ and $`k_z`$ according to Eq. (7). Evaluation of the summation takes much computer time, so we resort to some approximations. Since the cross section is a square, that is $`L_y=L_z`$, an excellent approximation to reduce the number of terms is to partition the transverse $`k_{}`$-plane by concentric circles with radius $`k_{}=\sqrt{k_y^2+k_z^2}`$ and a $`\delta k_{}=\frac{\pi }{L_t}`$. The number of allowed values of $`k_y`$ and $`k_z`$ in each ring is given by its surface divided by the density of the transverse $`k`$ states $`\left(\frac{\pi }{L_t}\right)^2`$. Another approximation can be made by noting that the terms involving $`k_{}>>\sqrt{\mu }`$ do not have significant contributions to the sum. This allows us to evaluate a finite number of terms instead of evaluating an infinite series. In our calculation we take $`k_{_{max}}=3\sqrt{\mu }`$. In determining the coupling constant $`V`$ according to the BCS-relation

$$T_c=1.13\omega _De^{\frac{1}{N(\mu )V}},\mathrm{with}N(\mu )=\frac{\mu }{4\pi ^2}$$
(36)
we use $`T_c=1.2`$ K and $`\omega _D=375`$ K for aluminum. We find $`V=9.516`$ Ry.

### A The Bar-Shaped Superconductor

The matrix Green’s function for a bar-shaped superconductor $`G_S^0(x,x,k_y,k_z,i\omega _n)`$ is given by Eq. (19). By substituting the upper right element of this Green’s function in Eq. (34) one straightforwardly obtains the following selfconsistency condition for the homogeneous bar-shaped superconductor,

$$\mathrm{\Delta }(x)=\frac{V|\mathrm{\Delta }|}{\beta L_yL_z}\underset{\omega _n,k_y,k_z}{}\frac{1}{4\mathrm{\Omega }_S}\left(\frac{1}{\sqrt{k_{F_x}^2+i\mathrm{\Omega }_S}}+\frac{1}{\sqrt{k_{F_x}^2i\mathrm{\Omega }_S}}\right).$$
(37)
In calculating the gap self-consistently, we first assume an initial value of the gap. By substituting it on the right-hand side of Eq. (37) we obtain a new value of the gap. Then this new value is again substituted on the right-hand side to get another new value. This procedure can be repeated until the difference between successive iterations is negligibly small. As shown in FIG. 9, the difference between the first ten iterations is still significant. After about 80 iterations, the value of the gap stabilizes to $`1.106\times 10^5`$ Ry. The transverse length is $`1000`$ Bohr and the temperature is $`0.6`$ K. The initial value of the gap is $`2.0\times 10^5`$ Ry.

FIG. 10 shows the plot of the gap against the transverse length $`L_t`$ for different temperatures. The number of iterations is $`100`$. It can be seen that there are oscillations of the gap. The amplitudes of the oscillations decrease as the transverse dimension increases. These oscillations can be attributed to the discreteness of the transverse wave vector. As the transverse width increases, the transverse wave vector approaches the continuous regime which can be gauged from the gap becoming closer to its bulk value obtained by integrating instead of summing over the transverse wave vectors. In the figure, we show the bulk value at $`0.2`$ K. Another interesting thing which can be seen in the figure is the suppression of superconductivity for narrower transverse dimensions. We notice that as the temperature increases the onset of the suppression of superconductivity occurs at higher values of $`L_t`$.

In FIG. 11 we show the temperature variation of the gap for different transverse widths $`L_t`$. A residual value of the gap for the bulk superconductor can still be observed beyond the critical temperature ($`1.2`$ K for aluminum). We also notice small-amplitude oscillations of the gap at higher temperatures, which only show up in the curves for larger transverse dimensions. For smaller transverse dimensions, the oscillations are suppressed. These oscillations are due to the cut-off in the summation over the Matsubara frequencies. For lower temperatures, the cut-off value $`n_{max}`$ in Eq. (35) becomes very large and the results are no longer sensitive to it.

### B The NS System

The matrix Green’s function appropriate for a normal metal-superconductor system is given in Eqs. (21) and (22). The first term on the right-hand side of Eq. (21) is the matrix Green’s function for the bar superconductor. The second term contains the elements of the $`t`$-matrix, which take into account the scattering of the quasiparticles at the interface. The latter acts as a perturbing term to the former and is therefore responsible for the spatial variations of the gap and the pair amplitude at the vicinity of the interface. To calculate the spatial variation of the gap according to Eq. (34), we substitute the upper right element of the matrix Green’s function, Eq. (21), using the value of the self-consistent gap for the bar-shaped superconductor.

Before presenting our results, we want to mention a computational problem which has to be solved. Due to the explicit presence of the position $`x_j`$ in Eq. (23), the matrices involved have alternating columns of very small and very large values which the computer cannot handle anymore. However this problem is not intrinsic to the formalism and can be solved by appropriate rescaling. By defining the matrix $`\widehat{t}`$

$$\widehat{t}_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}e^{i\sigma \nu k_{\nu j}^\sigma x_j}t_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}e^{i\sigma ^{}\nu ^{}k_{\nu ^{}j}^\sigma ^{}x_j}$$
(38)
the Green’s function (21), combined with Eq. (22), obtains the form

$`G_{\nu j\nu ^{}j}(x,x^{},k_y,k_z,i\omega _n)`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}d_{\nu j}^\sigma \psi _{\nu j}^{\sigma \mu }(x)\stackrel{~}{\psi }_{\nu j}^{\sigma ,\mu }(x^{})\delta _{\nu \nu ^{}}`$ (40)
$`+{\displaystyle \underset{\sigma \sigma ^{}}{}}d_{\nu j}^\sigma d_{\nu ^{}j^{}}^\sigma ^{}\psi _{\nu j}^{\sigma \nu }(xx_j)\widehat{t}_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}\stackrel{~}{\psi }_{\nu ^{}j}^{\sigma ^{}\nu ^{}}(x^{}x_j)`$
in which only position differences occur. The rescaled matrix $`\widehat{t}`$ is determined by the equation

$$\underset{\sigma \nu }{}\nu d_{\nu j}^\sigma \psi _{\nu j}^{\sigma \nu }(0)\widehat{t}_{\nu j\nu ^{}j}^{\sigma \sigma ^{}\nu \nu ^{}}=\nu ^{}\psi _{\nu ^{}j}^{\sigma ^{}\nu ^{}}(0),$$
(41)
found straightforwardly from the original equation (23).

FIG. 12(a) shows the spatial variation of the gap near the interface of a normal metal-superconductor system at $`T=0.6`$ K for different transverse widths. At a distance of 30000 Bohr from the interface, which is about six times the coherence length ($`4500`$ Bohr), the gap is slightly smaller than the one obtained for the bar-shaped superconductor. The differences are about $`7.79\%`$, $`6.57\%`$, and $`2.93\%`$ of their bar value for $`L_t=99.8514`$, $`L_t=100`$, and $`1000`$ Bohr, respectively. At $`4000`$ Bohr, which is of the order of the coherence length from the interface, the differences are about $`21.7\%`$, $`20.44\%`$, and $`14.71\%`$, respectively. These figures lead to the inevitable conclusion that for larger transverse dimensions, the proximity effect is less pronounced than for smaller ones. This may be due to the fact that for small transverse dimensions the superconductivity tends to be suppressed. We find that the difference seen for the two smaller widths can be attributed to the oscillations in the self-consistent gap shown in FIG. 10. We have not seen a special influence of the fact that the smaller width is a critical one. FIG. 12(b) shows the corresponding pair amplitude defined by $`\frac{\mathrm{\Delta }(x)}{V}`$ (see Eq. (33)) which has a finite value in the normal metallic region near the NS interface but it decays in the inner region of the normal metal.

### C The SNS System

In the superconductor-normal metal-superconductor system, there are two interfaces which we designate as $`x_1`$ and $`x_2`$. The $`t`$-matrices must be determined at these interfaces so that we can evaluate the $`T`$-matrix. The unpleasant singularities due to the explicit presence of the position of the interfaces in Eq. (23) can be removed by using Eq. (38) resulting to a corresponding transformation of the $`T`$-matrix given by

$$\widehat{T}_{\nu j\nu ^{}j^{}}^{\sigma \sigma ^{}\nu \nu ^{}}=e^{i\sigma \nu k_{\nu j}^\sigma x_j}T_{\nu j\nu ^{}j^{}}^{\sigma \sigma ^{}\nu \nu ^{}}e^{i\sigma ^{}\nu ^{}k_{\nu ^{}j^{}}^\sigma ^{}x_j^{}}.$$
(42)
By implementing this, the multiple-scattering Green’s function (26) and Eq. (27), which determines the $`T`$-matrix in terms of the $`t`$-matrices, can be modified straightforwardly.

The steps in calculating the gap and the pair amplitude self-consistently are the same as in Sec. IV B. FIG. 13 shows the self-consistent gap function and the pair amplitude for different transverse widths. The center of the system is at $`x=0`$. The spatial variation of the gap near the interfaces is clearly shown. The proximity effect is stronger than for the NS system shown in FIG. 12. Whereas in FIG. 12a for $`L_t=1000`$ Bohr the gap at a distance of $`6000`$ Bohr from the interface has a value of $`9.753\times 10^6`$ Ry, in FIG. 13a it has already increased to $`1.074\times 10^5`$ Ry, which is much closer to the bulk value of $`1.106\times 10^5`$ Ry. In FIG. 13b the pair amplitude in the N region does not decrease below a value of $`5\times 10^7`$. In FIG. 12b it decreases to zero and the value at $`2000`$ Bohr from the NS interface has already decreased to $`2.784\times 10^7`$. In FIG. 14 we show the gap functions for different values of $`L`$. It can be seen that for larger $`L`$ the gap function is lesser in magnitude. This is another manifestation of the proximity effect.

### D The SNSNS System

The extension of the steps outlined in Sec. IV C leading to the $`T`$-matrix and the matrix Green’s function for a SNSNS system is straightforward. In this case we have to consider four interfaces but the procedure is basically the same. In the outer superconductors, we again use the self-consistent value of the gap for the superconducting bar we have calculated in Sec. IV A. In calculating the gap in the middle superconductor self-consistently, we apply the recipe introduced by Tanaka and Tsukada . We start the iteration procedure by taking half of the gap in the outer superconductors for the gap in the middle. By substituting these values in the righthand side of Eq. (34), the middle $`x`$-dependent gap function is obtained and its spatial average is determined. This average value is used as the new value of the gap in the middle superconductor. The process is repeated until the difference of these average values between successive iterations is negligibly small. To determine the gap’s spatial variation, the self-consistent average value of the gap in the middle superconductor and the bar value in the outer ones are substituted in the right-hand side of Eq. (34).

FIG. 15 shows the gap and the pair amplitude of two SNSNS systems, one with a length $`L_N=3000`$ Bohr for the normal metallic parts and the other one with a larger length $`L_N=5000`$ Bohr. First note that the gap of the middle superconductor is smaller than the gap of the outer superconductors although the metallic parameters for these superconductors are taken equal. This supports our choice of a $`h`$ value smaller than $`1`$ in Sec. III, see FIG. 2. Further note that the gap and the pair amplitude are higher for the system with narrower normal-metal part, which is just another manifestation of the proximity effect. Although by definition the gap is zero in the normal-metal parts, the pair amplitude is different from zero. This is also a manifestation of the proximity effect. If the width of the normal-metal layers had been much larger than the sum of the coherence lengths of the outer and middle superconductors, the pair amplitude would have been zero there except in those regions near the interfaces.

## V Conclusions and Future Prospects

The multiple-scattering Green’s function formalism developed by Koperdraad  has been applied to determine the Andreev bound states in SNS and SNSNS junctions through the local density of states and to calculate the superconducting gap function self-consistently. We have shown that for transverse junction widths tuned such that the motion in the longitudinal direction is very slow, the Andreev approximation breaks down. For these transverse dimensions, the highest mode is supported by many bound states whose degeneracy is lifted when no approximation is applied. The thickness used for the normal metallic layers is chosen such that the lower modes are each supported by only one nondegenerate state. Results for the self-consistent gap functions exhibit various features of the proximity effect. Furthermore, our results show that for small transverse dimensions, superconductivity is suppressed.

The formalism is applicable to systems of an arbitrary number of layers. In addition, it allows for the calculation of the supercurrents through such junctions, to which self-consistent gap functions are necessary . This will be investigated in the near future.

## Caption for figures

FIG. 1. The geometry of the bar-shaped superconductor under consideration. It has finite transverse dimension and is infinite longitudinally. A rectangular cross-section is shown whose dimensions are $`L_y`$ and $`L_z`$.

FIG. 2. A two-dimensional picture of the systems studied. In the upper panel the vertical direction stands for a transverse direction. In the lower panels the potential is shown, which is zero in the normal part(s), and proportional to $`\mathrm{\Delta }`$ in the superconducting parts.

FIG. 3. Plot of the LDOS against $`E/\mathrm{\Delta }`$ for a SNS system in which $`L_t=13`$ Bohr.

FIG. 4. Plot of the LDOS against $`E/\mathrm{\Delta }`$ for $`x=1000`$ Bohr, $`L_t=12.5676`$ Bohr, L=4000 Bohr, and $`\mathrm{\Delta }=0.0001`$ Ry.

FIG. 5. Plot of the LDOS for a SNS system against $`E/\mathrm{\Delta }`$ at $`x=1000`$ Bohr (solid curve) and $`x=1500`$ Bohr (dashed curve). The transverse dimension is $`L_t=12.5676`$ Bohr and the length of the normal metal is $`L=4000`$ Bohr.

FIG. 6. Plot of the LDOS against $`E/\mathrm{\Delta }`$ for a SNS system in which $`L_t=100`$ Bohr (solid curve) and $`L_t=99.8514`$ Bohr (solid curve). The length of the normal-metal part is $`4000`$ Bohr.

FIG. 7. Plot of the LDOS at $`x=1500`$ Bohr against $`E/\mathrm{\Delta }`$ for a SNSNS system in which $`L_t=12.5676`$ Bohr for (a) $`h=0.25`$ (b) $`h=0.5`$ (c) $`h=0.75`$ and (d) $`h=1`$.

FIG. 8. The LDOS against $`E/\mathrm{\Delta }`$ for a SNSNS system in which $`L_t=12.5676`$ Bohr. In (a) $`L=6000`$ Bohr, the length of the middle superconductor is $`4000`$ Bohr and the LDOS is calculated at $`x=2500`$ Bohr. In (b) $`L=8000`$ Bohr, the length of the middle superconductor is $`6000`$ Bohr and the LDOS is calculated at $`x=3500`$ Bohr. The gap of the middle superconductor is $`0.25\mathrm{\Delta }`$. All unlabeled LDOS peaks belong to the mode $`(2,2)`$. The calculation is done in the Andreev approximation.

FIG. 9. The gap against the number of iterations for a bar-shaped superconductor. Note that the value of the gap stabilizes as the number of iterations increases. For the system considered, $`L_t=1000`$ Bohr and $`T=0.6`$ K.

FIG. 10. Plot of the self-consistent gap function against the transverse length $`L_t`$ for a bar-shaped superconductor at different temperatures. The number of iterations is 100.

FIG. 11. The temperature variation of the self-consistent gap for different transverse widths.

FIG. 12. (a) The gap and (b) the pair amplitude against the distance from the interface of a NS system at $`T=0.6`$ K. The interface is chosen at $`x=0`$.

FIG. 13. (a) The gap function and (b) the pair amplitude of a SNS system against the distance from the middle of the system chosen at $`x=0`$. The interfaces are located at $`x=\pm 2000`$.

FIG. 14. The gap function of a SNS system against the distance from the middle of the system chosen at $`x=0`$ for different values of $`L`$. The interfaces are located at $`\pm 2000`$ for $`L=4000`$ Bohr and at $`\pm 4000`$ for $`L=8000`$ Bohr.

FIG. 15. (a) The gap function and (b) the pair amplitude against the distance from the middle of a SNSNS system, which is chosen at $`x=0`$. The interfaces are chosen at $`x=\pm 2000,\pm 5000`$ for $`L_N=3000`$ Bohr (solid curve) and $`x=\pm 2000,\pm 7000`$ for $`L_N=5000`$ Bohr (dotted curve). The transverse width is $`L_t=100`$ Bohr.


text:
# The rare decay 𝑡→𝑐𝛾 in the General 2HDM type III

## I Introduction

The search for Flavor Changing Neutral (FCN) processes is one of the most interesting possibilities to test the Standard Model (SM), with the potential for either discovering or putting stringent bounds on new physics. In the SM, there are not FCN processes at the tree-level and at one-loop level they are GIM suppressed gim . For this reason, the branching ratios predicted within the SM for these kind of processes are in general, very small. Special interest has been devoted to the so-called rare top quark decays, e.g. $`tc(\gamma ,Z,g)`$ and $`tc(h^0,H^0,A^0)`$ sola . The reason for the interest in these decays is again that the observation of a single event of this kind in near future colliders would imply evidence of new physics which could be related with the Higgs sector of the theory.

In the framework of the SM the decay $`tbW`$ is by far the dominant one and the rates for the rare decays of the top quark into gauge bosons $`(\gamma ,Z,g)`$ plus c-quark are too much small. With the current experimental value of the top quark mass, the ratio for $`tc\gamma `$ is of the order of $`5\times 10^{13}`$, for the $`Z^0`$ boson $`1\times 10^{13}`$ and for the gluon $`4\times 10^{13}`$ sm . These numbers are so tiny because in the SM the FCN decays of the top quark are present at the one loop-level and they are controlled by down-type quarks and are beside strongly suppresed by the GIM mechanism gim .

But some of this undesirable features can change with physics beyond the SM. For instance, there is not any fundamental reason for having only one Higgs doublet. In fact, for two or more scalar doublets, FCN couplings could be generated at tree level unless and $`adhoc`$ discrete symmetry is imposed gw . Then it is possible to enhance the branching ratios for rare top quark decays in different kinds of models which include new physics.

Models like Minimal Supersymmetric Standard Model (MSSM) with R-parity, that relax the universality of soft breaking terms with a large flavor mixing only between second and third family can reach values such as $`B(tc\gamma (Z))10^8`$ and $`B(tcg)10^6`$ mssm . On the other hand, it is possible to introduce new couplings with R-parity broken and in this context, some of these channels could become observable. For example, very enhanced ratios are predicted in this framework such as $`B(tc\gamma )10^5`$, $`B(tcZ)10^4`$ and $`B(tcg)10^3`$ rpar  which could be within the observable threshold of near future experiments. MSSM needs a new Higgs doublet and therefore its Higgs sector is phenomenologically richer. Then, one possibility for new physics is that the electroweak symmetry breaking involves more that one Higgs doublet and without supersymmetry. In particular, if there is only one extra Higgs doublet, we are in the context of the so-called Two Higgs Doublet Model (2HDM). In such models, a new free parameter $`\mathrm{tan}\beta `$ arise from the breaking symmetry sector and as we shall see later, it could provide a significant enhancement for the FCN processes. There are three classes of these 2HDM which have been used to examine the ratios for rare top decays, called models I, II and III. The first two are characterized by a discrete symmetry which avoid tree-level FCN processes gw . In the model type I, one Higgs Doublet provides masses to the up-type and down-type quarks, simultaneously. In the model type II, one Higgs doublet gives masses to the up-type quarks and the other one to the down-type quarks. In models type I and II the corresponding branching fractions for $`tc(\gamma ,Z,g)`$ cannot approach in anyway the detectability threshold reaching only 3-4 orders of magnitude above the SM values sm .

But the discrete symmetry gw  is not compulsory and both doublets may generate the masses of the quarks of up-type and down-type simultaneously, in such case we are in the model type III modb . It has been used to search for FCNC at the tree level ARS , sher . In next sections we will show that in the context of the general 2HDM type III, the $`\mathrm{tan}\beta `$ enhancement might be enough to achieve the dectectability threshold. Recently, the 2HDM type III has been discussed and classified us , depending on how the basis for the vacuum expectation values (VEV) are chosen and according to the way in which the flavor mixing matrices are rotated. In brief, the reference us  considered two cases according to the VEV structure; in the case (a), both Higgs Doublets acquire a VEV, while in the case (b) only one of them get it. Further, there are two types of rotations which generate four different Lagrangians in the quark sector and two different ones in the leptonic sector. The well known 2HDM types I and II, could be generated from them in the limit in which the mixing vertices at the tree level vanish. In that paper, has been pointed out that the phenomenology of the 2HDM type III is highly sensitive to the rotation used for the mixing matrices. Also it was shown that case (b) is a particular ocurrence of case (a). The model which includes all these features mentioned above, will be called from now on the General Two Higgs Doublet Model (G2HDM) type III. The case (b) is the most familiar in the literature and calculations has been done ARS  for the branching fractions of top quark decays. In this particular version of model type III values up to $`B(tc\gamma )10^7,B(tcZ)10^6`$ and $`B(tcg)10^5`$ are gotten ARS .

On the other hand, we must take into account the sensitivity of future experiments to these kind of processes and we must prepare the analysis of the phenomenology in different models in order to distinguish among them in case of a single event detection of this kind. Up to now, CDF has set limits on the FCN decays $`tc(g,Z,\gamma )`$. In seeking $`tZq`$ they looked at $`p\overline{p}t\overline{t}qZbW,qZbZll+4jets`$ with high jet $`E_T`$, taking into account the possible background from the SM reactions, they are getting the upper limit $`B(tcZ)<0.33`$ at 95% C.Lcdf1 . To study $`tq\gamma `$, CDF has examined the reaction $`p\overline{p}t\overline{t}Wb\gamma q`$. If the W decayed leptonically, the signature would be $`\gamma +l+E_{\mathit{}}+(2)`$ jets, but if it were hadronically, the signature would be $`\gamma +(4)`$ jets with one jet b-tagged. CDF set the 95 % C.L. upper bound $`B(t\gamma q)<0.032`$ cdf1 . Run II will provide much greater sensitivity to these decays, getting to explore ratios up to $`B(tq\gamma )8.4\times 10^5`$ and $`B(tZq)6.3\times 10^4`$ at 100 $`fb^1`$ of the integrated luminosity cdf2 . Other analysis for the sensitivity to the FCN decays of top quark have been made for the LHC collider lhc . For $`tZq`$, ATLAS collaboration could achieve $`B(tZq)2\times 10^4`$ at the $`5\sigma `$ level and CMS collaboration up to $`B(tZq)1.9\times 10^4`$ at the $`5\sigma `$ level, assuming an integrated luminosity of 100 $`fb^1`$. Now, for the decay $`t\gamma q(u,c)`$, ATLAS could go up to $`B(t\gamma q)10^4`$ while CMS about $`B(t\gamma q)3.4\times 10^5`$ at $`5\sigma `$ level.

Our aim is to calculate the branching ratio for $`tc\gamma `$ in the framework of the most general CP-conserving 2HDM type III, and it shall be shown that in this model is possible to reach the limits of the sensitivity expected in near future top quark experiments. This branching ratio was considered by Atwood, et. al. ARS  in the 2HDM type III but they choose an specific set of values for the VEV, so that the parameter $`\mathrm{tan}\beta `$ does not appear in the analysis, which we call case (b). We have shown previously us  that the CP-conserving 2HDM type III with both VEV different from zero has two interesting features. First, it is possible to generate the Lagrangians for the models type I and II choosing conveniently the basis for the rotation of the mixing matrices in the Higgs sector. Second, in this model appears explicitly the parameter $`\mathrm{tan}\beta `$ which is the ratio of the vacuum expectation values. Such parameter plays an important role in the phenomenology of the model as we will show in next sections. In fact, this parameter appears as an important one in SUSY versions of the SM in a natural way, and there is not any fundamental reason to avoid its appearence in the 2HDMs type III. In order to give a better presentation we reviewed the general 2HDM type III in section II and then we present the analysis and concluding remarks for the FCN decays of the top quark in section III.

## II The General 2HDM type III

The most general Yukawa’s Lagrangian which includes a new Higgs doublet, is as follow

$`\mathrm{\pounds }_Y`$ $`=`$ $`\eta _{ij}^U\overline{Q}_{iL}\stackrel{~}{\mathrm{\Phi }}_1U_{jR}+\eta _{ij}^D\overline{Q}_{iL}\mathrm{\Phi }_1D_{jR}+\eta _{ij}^E\overline{l}_{iL}\mathrm{\Phi }_1E_{jR}`$
$`+`$ $`\xi _{ij}^U\overline{Q}_{iL}\stackrel{~}{\mathrm{\Phi }}_2U_{jR}+\xi _{ij}^D\overline{Q}_{iL}\mathrm{\Phi }_2D_{jR}+\xi _{ij}^E\overline{l}_{iL}\mathrm{\Phi }_2E_{jR}+h.c.`$
where $`\mathrm{\Phi }_{1,2}`$are the Higgs doublets,$`\eta _{ij}^0`$and $`\xi _{ij}^0`$ are non-diagonal $`3\times 3`$ matrices and $`i`$, $`j`$ are family indices.

In this kind of model (type III), we consider two cases. In the case (a) we assume the VEV as

$$\mathrm{\Phi }_1_0=\left(\begin{array}{c}0\\ v_1/\sqrt{2}\end{array}\right),\mathrm{\Phi }_2_0=\left(\begin{array}{c}0\\ v_2/\sqrt{2}\end{array}\right).$$

The mass eigenstates of the scalar fields are given in ref.moda . The case (b) corresponds to the case in which one of the VEV is equal to zero,

$$\mathrm{\Phi }_1_0=\left(\begin{array}{c}0\\ v_1/\sqrt{2}\end{array}\right),\mathrm{\Phi }_2_0=0.$$

The mass eigenstate scalar fields in this case are given in ref. modb . An important difference between both models is that the would-be Goldstone bosons $`G_{Z\left(W\right)}`$are a linear combination of components of $`\mathrm{\Phi }_1`$and $`\mathrm{\Phi }_2`$ in the model (a), meanwhile in the model (b) are a component of the doublet $`\mathrm{\Phi }_1`$ us .

To convert the Lagrangian (II) into mass eigenstates we make the unitary transformations from which we obtain the mass matrices. In the context of case (a)

$`M_D^{diag}`$ $`=`$ $`V_L\left[{\displaystyle \frac{v_1}{\sqrt{2}}}\eta ^{D,0}+{\displaystyle \frac{v_2}{\sqrt{2}}}\xi ^{D,0}\right]V_R^{},`$
$`M_U^{diag}`$ $`=`$ $`T_L\left[{\displaystyle \frac{v_1}{\sqrt{2}}}\eta ^{U,0}+{\displaystyle \frac{v_2}{\sqrt{2}}}\xi ^{U,0}\right]T_R^{}.`$ (2)
Where $`V_{L(R)}`$ and $`T_{L(R)}`$ rotate down-type quarks and up-type quarks, respectively.

We can solve for $`\xi ^{D,0},\xi ^{U,0}`$ (rotations of type I). Replacing them into (II) the expanded neutral Lagrangian for up and down sectors are

$`\mathrm{\pounds }_{Y\left(Q=U,D\right)}^{(a,I)}`$ $`=`$ $`{\displaystyle \frac{g}{2M_W\mathrm{sin}\beta }}\overline{Q}M_Q^{diag}Q\left(\mathrm{sin}\alpha H^0+\mathrm{cos}\alpha h^0\right)`$ (3)
$``$ $`{\displaystyle \frac{1}{\sqrt{2}\mathrm{sin}\beta }}\overline{Q}\eta ^QQ\left[\mathrm{sin}\left(\alpha \beta \right)H^0+\mathrm{cos}\left(\alpha \beta \right)h^0\right]`$
$``$ $`{\displaystyle \frac{ig\mathrm{cot}\beta }{2M_W}}\overline{Q}M_Q^{diag}\gamma _5QA^0\pm {\displaystyle \frac{i}{\sqrt{2}\mathrm{sin}\beta }}\overline{Q}\eta ^Q\gamma _5QA^0`$
$``$ $`{\displaystyle \frac{ig}{2M_W}}\overline{Q}M_Q^{diag}\gamma _5QG_Z^0`$
where $`\eta ^{U(D)}=T_L(V_L)\eta ^{U(D),0}T_R^{}(V_R)^{}`$ and similarly for $`\xi ^{U(D)}`$. The superindex $`(a,I)`$ refers to the case (a) and rotation type I. It is easy to check that the Lagrangian constructed as $`\mathrm{\pounds }_{Y(U)}^{(a,I)}+\mathrm{\pounds }_{Y(D)}^{(a,I)}`$ is the Lagrangian of the 2HDM type I moda , plus FC interactions. Moreover,, it is observed that the case (b) in both up and down sectors can be calculated just taking the limit $`\mathrm{tan}\beta \mathrm{}`$.

On the other hand, from (II), we can also solve for $`\eta ^{D,0},\eta ^{U,0}`$ instead of $`\xi ^{D,0},\xi ^{U,0}`$ which we call rotations of type II. Replacing them into (II) the expanded neutral Lagrangian for up and down sectors become

$`\mathrm{\pounds }_{Y\left(Q=U,D\right)}^{(a,II)}`$ $`=`$ $`{\displaystyle \frac{g}{2M_W\mathrm{cos}\beta }}\overline{Q}M_Q^{diag}Q(\mathrm{cos}\alpha H^0\mathrm{sin}\alpha h^0)`$ (4)
$`+`$ $`{\displaystyle \frac{1}{\sqrt{2}\mathrm{cos}\beta }}\overline{Q}\xi ^QQ\left[\mathrm{sin}(\alpha \beta )H^0+\mathrm{cos}(\alpha \beta )h^0\right]`$
$`\pm `$ $`{\displaystyle \frac{ig\mathrm{tan}\beta }{2M_W}}\overline{Q}M_Q^{diag}\gamma _5QA^0{\displaystyle \frac{i}{\sqrt{2}\mathrm{cos}\beta }}\overline{Q}\xi ^Q\gamma _5QA^0`$
$``$ $`{\displaystyle \frac{ig}{2M_W}}\overline{Q}M_Q^{diag}\gamma _5QG_Z^0.`$
In this situation the case (b) is obtained in the limit $`\mathrm{tan}\beta 0`$. Moreover, the Lagrangian built as $`\mathrm{\pounds }_{Y(U)}^{(a,I)}+\mathrm{\pounds }_{Y(D)}^{(a,II)}`$ coincides with the one of the 2HDM type II moda , plus some FC interactions. In addition, we can build two additional lagrangians by adding $`\mathrm{\pounds }_{Y(U)}^{(a,II)}+\mathrm{\pounds }_{Y(D)}^{(a,II)}`$and $`\mathrm{\pounds }_{Y(U)}^{(a,II)}+\mathrm{\pounds }_{Y(D)}^{(a,I)}`$. So four different Lagrangians are generated from the case (a). Therefore, with the two Lagrangians (3) and (4) is possible to describe all possible CP conserving neutral current interactions in the quark sector of any model that includes only one extra Higgs doublet as Physics beyond the SM.

## III The rare decay $`tc\gamma `$

In the present work, we calculate the branching fraction for the rate of the top quark decay $`tc\gamma `$ using the model presented in section II. The widths of $`tc\gamma `$ for the rotations type I and II are

$`\mathrm{\Gamma }^{a,I}(tc\gamma )`$ $`=`$ $`{\displaystyle \frac{m_t^3G_F\alpha _{em}}{1152\sqrt{2}\pi ^4}}{\displaystyle \frac{|\eta _{tc}|^2}{\mathrm{sin}^4\beta }}|(\mathrm{sin}\alpha \mathrm{sin}(\alpha \beta )+\widehat{m}_w\eta _{tt}^{}\mathrm{sin}^2(\alpha \beta ))F_2(\widehat{m}_{H^0}^2)`$
$`+`$ $`(\mathrm{cos}\alpha \mathrm{cos}(\alpha \beta )+\widehat{m}_w\eta _{tt}^{}\mathrm{cos}^2(\alpha \beta ))F_2(\widehat{m}_{h^0}^2)(\mathrm{cos}\beta \widehat{m}_w\eta _{tt}^{})F_2(\widehat{m}_{A^0}^2)|^2,`$
$`\mathrm{\Gamma }^{a,II}(tc\gamma )`$ $`=`$ $`{\displaystyle \frac{m_t^3G_F\alpha _{em}}{1152\sqrt{2}\pi ^4}}{\displaystyle \frac{|\xi _{tc}|^2}{\mathrm{cos}^4\beta }}|(\mathrm{cos}\alpha \mathrm{sin}(\alpha \beta )+\widehat{m}_w\xi _{tt}^{}\mathrm{sin}^2(\alpha \beta ))F_2(\widehat{m}_{H^0}^2)`$
$`+`$ $`(\mathrm{sin}\alpha \mathrm{cos}(\alpha \beta )+\widehat{m}_w\xi _{tt}^{}\mathrm{cos}^2(\alpha \beta ))F_2(m_{\widehat{h}^0}^2)(\mathrm{sin}\beta +\widehat{m}_w\xi _{tt}^{})F_2(\widehat{m}_{A^0}^2)|^2`$
where

$$F_2(\widehat{m}_h^2)=C_{21}(\widehat{m}_h^2)+C_{23}(\widehat{m}_h^2)C_{12}(\widehat{m}_h^2)2C_{11}(\widehat{m}_h^2)$$

with $`C_{ij}(\widehat{m}_h^2)=C_{ij}(P^2=\widehat{m}_t^2,0,0,\widehat{m}_t^2,\widehat{m}_h^2)`$ the Veltman-Passarino integrals and the definitions $`\widehat{m}_i^2=m_i^2/m_t^2`$, $`\eta ^{}(\xi ^{})_{ij}=\eta (\xi )_{ij}/g`$.

We have taken into account the contributions due to the diagonal and nondiagonal terms of $`\eta (\xi )^{U,D}`$ but we parametrize them using the prescription of Cheng and Sher sher , where the FC couplings are proportional to the geometric average of the quark masses, $`\eta (\xi )_{ij}g\sqrt{m_im_j}/2m_W`$.

In figure 1, we have plotted the branching ratio $`B(tc\gamma )`$ for the rotations type I and II. We set $`\alpha =0`$, $`m_H=m_h`$ and $`m_{A^0}`$ high enough to have only one scalar mass varying in the plot. The plot shows the behaviour with $`\mathrm{tan}\beta =0.1`$ (dashed line) for rotation type I and $`\mathrm{tan}\beta =15`$ (solid line) for rotation type II. We see that for reasonable values of $`\mathrm{tan}\beta `$ is possible to reach values in the range of the sensitiviy of future top quark experiments. This is because the expressions for the widths in rotation type I and II have a strong dependence on the $`\mathrm{tan}\beta `$ parameter in such way that in rotation type I we can enhance the branching ratio choosing $`\mathrm{tan}\beta `$ small and in rotation type II we can do it setting $`\mathrm{tan}\beta `$ large enough. We also point out that in the case of rotation type II with $`\mathrm{tan}\beta =0`$ and $`\alpha =0`$ we reproduce the results already obtained by Atwood, et.al. ARS .

In figure 2, we have plotted the branching ratio $`B(tc\gamma )`$ versus $`\alpha `$ for rotation type II with $`m_H=m_h=150`$ GeV and $`\mathrm{tan}\beta =8`$. We can see that it is possible to get values of the order of the sensitivity expected by coming soon experiments in top decays $`10^5`$ with even lower values of $`\mathrm{tan}\beta `$. And, in fig. 3 we are showing the lower bounds of $`\mathrm{tan}\beta `$ in rotation type II according with the value of the scalar Higgs boson mass which are enough to reach the experimental discovery range, i.e. for $`B(tc\gamma )=10^5(10^4)`$ solid line (dashed line).

In conclusion, we can observe that the bounds obtained in this paper are highly sensitive to the type of rotation. Also we might notice that in the 2HDMs type III $`\mathrm{tan}\beta `$ should be taken into account due to its importance as a free parameter related with the symmetry breaking sector. Instead of that, we are showing that by driving $`\mathrm{tan}\beta `$ properly, it is possible to enhance the rare branching $`B(tc\gamma )`$. Further, it is worth to remark that this enhancement in 2HDM type III could also appear in other rare top quark decays $`tc(Zorg)`$, to achieve the corresponding experimental thresholds. Specifically, for values of $`B(tc\gamma )19\times 10^5`$ which are the orders of sensitivity of the experiments, we get the allowed region $`8\mathrm{tan}\beta 15`$ for rotation type II.

This work was supported by COLCIENCIAS, DIB and DINAIN.


text:
# Neutrino Oscillations Induced by Space-Time Torsion

## 1 Introduction

Neutrinos are the only particles for which a poem has been written; ”Cosmic Gall” by John Updike, quoted in . Recently neutrinos have attracted a lot of attention in high energy physics; see  and references therein. A major problem at present is the problem of solar neutrinos. Thermonuclear reactions in the sun produce large number of electron neutrinos; the standard solar model allows us to predict the flux of these neutrinos approximately as $`7.3\pm 2.3`$ SNU111SNU stands for solar neutrino unit and 1 SNU = 1 capture/sec/$`10^{36}`$ target atoms . On the other hand experiments which have been done for decades give a measured flux of $`2.55\pm 0.17\pm 0.18`$ SNU – approximately one third of the theoretical rate. This is the solar neutrino problem. A much-publicised solution is ”neutrino oscillation (or mixing)”. According to neutrino oscillation, in contrast to the conventional view that neutrinos are massless particles, they actually have a mass different from zero, although very small. The essential idea belongs to B.Pontecorvo . Technically speaking it is that the neutrino states of definite mass do not coincide with the weak interaction eigenstates. In other words, an electron neutrino emitted by sun becomes a linear combination of all three neutrinos (electron, muon and tau) while travelling in space towards the earth. Therefore the probability of detecting $`\nu _e`$ emitted by sun as $`\nu _\mu `$ or $`\nu _\tau `$ on the earth is different from zero. This situation, of course, invokes physics beyond the standard model. Although there are three neutrino families, the essential physics of neutrino oscillations is illustrated very well by considering the interactions of only two of these, which we choose to be $`\nu _e`$ and $`\nu _\mu `$.

All the above arguments have been cast in Minkowski spacetime. However we know that we live in a curved spacetime – perhaps even in a curved spacetime with torsion. Therefore, in more recent years, physicists have turned their attention to specifically gravitational contributions to neutrino oscillations – see  and references therein.

In this paper we reinvestigate the effects of gravitation, but we consider a spacetime with torsion, as well as simply curvature; this is the Einstein-Cartan theory, or more specifically the Einstein-Cartan-Dirac theory in our case, since we will be considering Dirac particles in the spacetime. Some work along these lines has already been done , but our results do not entirely replicate these findings. The essence of our work is to calculate the dynamical phase of neutrinos, by finding the form of the Hamiltonian, H, from the Dirac equation in Riemann-Cartan spacetime. The phase then follows from the formula

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}=H\psi \psi (t)=e^{\frac{i}{\mathrm{}}{\scriptscriptstyle H𝑑t}}\psi (0).`$ (1)
In what follows, $`\psi `$ is a Dirac spinor and $`H`$ is a $`4\times 4`$ matrix; the exponential function is therefore well-defined, as usual, by the series expansion. The Hamiltonian $`H`$ will depend, for example, on momentum $`\stackrel{}{p}`$, and this is expressed not as a differential operator but simply as a vector.

After we have introduced briefly the Einstein-Cartan-Dirac theory in Sec.2 we find the hamiltonian of a Dirac particle in subsec.2.1. In sec.3 we investigate, in turn, the neutrino oscillations for azimuthal and radial motions. The conclusions are given in sec.4.

## 2 Einstein-Cartan-Dirac Theory

A Riemann-Cartan space-time is defined by the triple $`\{M,g,\}`$ where M is a 4-dimensional differentiable manifold, equipped with a Lorentzian metric $`g`$ and a metric compatible connection $``$. The metric tensor may be written in terms of orthonormal basis 1-forms $`\{e^a\}`$ as

$`g=\eta _{ab}e^ae^b`$ $`=`$ $`e^0e^0+e^1e^1+e^2e^2+e^3e^3,`$ (2)
while the connection is specified by a set of connection 1-forms $`\{\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}\}`$. Metric compatibility requires

$`D\eta _{ab}d\eta _{ab}\mathrm{\Gamma }_{}^{c}{}_{a}{}^{}\eta _{cb}\mathrm{\Gamma }_{}^{c}{}_{b}{}^{}\eta _{ac}=0,`$ (3)
so that $`\mathrm{\Gamma }_{ab}=\mathrm{\Gamma }_{ba}`$. The Cartan structure equations below define the torsion 2-forms $`\{T^a\}`$ and the curvature 2-forms $`\{R_{}^{a}{}_{b}{}^{}\}`$ of the space-time:

$`De^a`$ $``$ $`de^a+\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}e^b=T^a,`$ (4)
$`D\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}`$ $``$ $`d\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}+\mathrm{\Gamma }_{}^{a}{}_{c}{}^{}\mathrm{\Gamma }_{}^{c}{}_{b}{}^{}=R_{}^{a}{}_{b}{}^{}.`$ (5)
In this study d, D, $`ı_a`$, * denote the exterior derivative, the covariant exterior derivative, the interior derivative and the Hodge star operator, respectively. The local orthonormal frame $`\{X_a\}`$ is dual to the coframe $`\{e^a\}`$;

$`ı_{X_a}e^bı_ae^b=e^b(X_a)=\delta _a^b.`$ (6)
The space-time orientation is set by the choice $`ϵ_{0123}=+1`$. In addition, the Dirac matrices $`\gamma ^a`$ are in the standard Bjorken and Drell form  but for our metric convention they satisfy

$`\gamma ^a\gamma ^b+\gamma ^b\gamma ^a=2\eta ^{ab}1_{4\times 4}.`$ (7)
We use the formalism of Clifford algebra-valued exterior forms for the Dirac particles. In this formalism $`\gamma =\gamma ^ae_a`$ and $`\gamma =\gamma ^ae_a`$. Furthermore, the covariant exterior derivative of spinor fields in the Riemann-Cartan space-time is

$`D\psi =d\psi {\displaystyle \frac{1}{8}}\mathrm{\Gamma }_{ab}[\gamma ^a,\gamma ^b]\psi .`$ (8)
The connection 1-forms, $`\{\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}\},`$ are decomposed as follows

$`\mathrm{\Gamma }_{}^{a}{}_{b}{}^{}=\omega _{}^{a}{}_{b}{}^{}+K_{}^{a}{}_{b}{}^{}.`$ (9)
Here, $`\{\omega _{}^{a}{}_{b}{}^{}\}`$ are Levi-Civita connection 1-forms which are given by

$`de^a+\omega _{}^{a}{}_{b}{}^{}e^b=0`$ (10)
and $`\{K_{}^{a}{}_{b}{}^{}\}`$ are kontorsion 1-forms which are given by

$`K_{}^{a}{}_{b}{}^{}e^b=T^a.`$ (11)
The Einstein-Cartan-Dirac Lagrangian density 4-form is

$`L={\displaystyle \frac{1}{2\mathrm{}^2}}R_{ab}(e^ae^b)+{\displaystyle \frac{i}{2}}\{\overline{\psi }\gamma D\psi +\overline{D\psi }\gamma \psi \}+{\displaystyle \frac{mc}{\mathrm{}}}\overline{\psi }\psi .`$ (12)
The field equations are then  ,

$`i\mathrm{}\gamma (D\psi {\displaystyle \frac{1}{2}}T\psi )+mc\psi `$ $`=`$ $`0,`$ (13)
$`R^{bc}(e_ae_be_c)`$ $`=`$ $`i\mathrm{}^2(\overline{\psi }\gamma ı_aD\psi (ı_a\overline{D\psi })\gamma \psi ),`$ (14)
$`T^c(e_ae_be_c)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2}}e_ae_b\overline{\psi }\gamma \gamma _5\psi ,`$ (15)
where $`T=ı_aT^a`$ is the torsion trace 1-form and $`\mathrm{}=\sqrt{\frac{Gh}{c^3}}10^{31}m`$ is the Planck length. As is seen from eqn.(15) the source of torsion is spin and magnitude of torsion is proportional to square of the Planck length. In this work, the spinorial contributions to the curvature will be assumed small, and therefore neglected. That is both the Dirac stress-energy tensor on the right hand side of eqn.(14) and the torsion dependent terms on the left hand side are dropped.

The Dirac equation (13) can be recast as

$`i\mathrm{}\gamma ^aı_a(D\psi {\displaystyle \frac{1}{2}}T\psi )mc\psi =0.`$ (16)
Here the torsion trace 1-form $`T`$ is cancelled by a similar term in the covariant exterior derivative and what is left behind are the totally anti-symmetric axial components of the torsion. These are going to be treated as given background fields.

### 2.1 Hamiltonian of a Dirac Particle

As mentioned above, our aim is to calculate the phase of a neutrino beam in Riemann-Cartan spacetime. We do this by finding the Hamiltonian for a Dirac particle. The spacetime we consider is Schwarzschild spacetime and the metric is therefore given by

$`g=f^2dtdt+{\displaystyle \frac{1}{f^2}}drdr+r^2(d\theta d\theta +\mathrm{sin}^2\theta d\phi d\phi ).`$ (17)
The orthonormal basis coframe is (from (2))

$`e^0=fcdte^1={\displaystyle \frac{1}{f}}dre^2=rd\theta e^3=r\mathrm{sin}\theta d\phi `$ (18)
where $`f^2=1\frac{2MG}{rc^2}`$. Using (18) and (10) we can calculate the Levi-Civita connection 1-forms, $`\omega _{}^{a}{}_{b}{}^{}`$. We then substitute eqn.(9) and the Levi-Civita 1-forms into eqn.(8) and then into eqn.(16). Now we use the identity

$`\gamma ^a[\gamma ^b,\gamma ^c]=2\eta ^{ab}\gamma ^c+2\eta ^{ac}\gamma ^b+2iϵ^{abcd}\gamma _5\gamma _d`$ (19)
and eqn.(11), and then write the contorsion 1-forms in component form as $`K_{}^{a}{}_{b}{}^{}=K_{}^{a}{}_{bc}{}^{}e^c`$. The resulting Dirac equation is

$`\{{\displaystyle \frac{i\mathrm{}}{cf}}\gamma ^0{\displaystyle \frac{}{t}}+i\mathrm{}f\gamma ^1({\displaystyle \frac{}{r}}+{\displaystyle \frac{1}{r}}+{\displaystyle \frac{f^{}}{2f}})+{\displaystyle \frac{i\mathrm{}}{r}}\gamma ^2({\displaystyle \frac{}{\theta }}+{\displaystyle \frac{\mathrm{cot}\theta }{2}})`$
$`+{\displaystyle \frac{i\mathrm{}}{r\mathrm{sin}\theta }}\gamma ^3{\displaystyle \frac{}{\phi }}mc+{\displaystyle \frac{\mathrm{}}{4}}K_{abc}ϵ^{abcd}\gamma _5\gamma _d\}\psi =0`$ (20)
where $`\gamma _5=i\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$. From the last term we see that it is only the totally anti-symmetric part of the contorsion tensor (in other words, the purely axial component) which contributes to the Dirac equation. \[Actually, in theories in which spin is the source of torsion, the contorsion tensor is in any case totally anti-symmetric .\] Relative to the orthonormal basis 1-forms $`\{e^a\}`$, we may write

$`K_{abc}=ϵ_{abcf}A^f`$ (21)
for some vector field222 From eqn.(18) $`A_0=\frac{1}{fc}A_t,A_1=fA_r,A_2=\frac{1}{r}A_\theta ,A_3=\frac{1}{r\mathrm{sin}\theta }A_\phi `$. This means $`\stackrel{}{A}=(A_1,A_2,A_3)`$ are not the cartesian components, rather the spherical components in the orthonormal basis.  $`A=A^aX_a`$. Noting that $`ϵ_{abcf}ϵ^{abcd}=3!\delta _f^d`$ and defining the components of the momentum operator by

$`p_r`$ $`=`$ $`i\mathrm{}({\displaystyle \frac{}{r}}+{\displaystyle \frac{1}{r}}+{\displaystyle \frac{f^{}}{2f}})`$ (22)
$`p_\theta `$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{r}}({\displaystyle \frac{}{\theta }}+{\displaystyle \frac{\mathrm{cot}\theta }{2}})`$ (23)
$`p_\phi `$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{r\mathrm{sin}\theta }}{\displaystyle \frac{}{\phi }},`$ (24)
the eqn.(20) becomes

$`i\mathrm{}\gamma ^0{\displaystyle \frac{\psi }{t}}`$ $`=`$ $`\{f^2\gamma ^1p_rc+f\gamma ^2p_\theta c+f\gamma ^3p_\phi c`$ (25)
$`+fmc^2+{\displaystyle \frac{3}{2}}\mathrm{}cf\gamma _5\gamma _aA^a\}\psi .`$
Using $`\beta \gamma ^0`$, $`\stackrel{}{\alpha }=\beta \stackrel{}{\gamma }`$ and $`\gamma ^0(\gamma _5\gamma _a)A^a=\gamma _5A^0+\stackrel{}{\mathrm{\Sigma }}.\stackrel{}{A}`$ where $`A^a=(A^0,\stackrel{}{A})`$ and $`\stackrel{}{\mathrm{\Sigma }}=\left(\begin{array}{cc}\stackrel{}{\sigma }& 0\\ 0& \stackrel{}{\sigma }\end{array}\right)`$. Comparing the above with the Schrödinger equation

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}=H\psi `$ (26)
we find that

$`H=f^2\alpha ^1p_rc+f\alpha ^2p_\theta c+f\alpha ^3p_\phi c+f\beta mc^2+{\displaystyle \frac{3}{2}}\mathrm{}cfA^0\gamma _5+{\displaystyle \frac{3}{2}}\mathrm{}cf\stackrel{}{\mathrm{\Sigma }}.\stackrel{}{A}.`$ (27)
## 3 Neutrino Oscillations

Neutrinos that are produced in the sun arrive on the earth. According to the theory of neutrino oscillations, a neutrino with fixed flavour which is emitted by sun becomes a superposition (in our model) of two neutrinos, i.e., they mix during the journey. In this case, even if only one kind of neutrino were emitted by the sun, the probability of measuring the other kind of neutrino on the earth would not be zero. This is called ”neutrino oscillation” and is measured by the interference pattern observed on the earth. The starting point of neutrino oscillation theory is that flavour eigenstates are different from mass eigenstates; each flavour eigenstate is a linear superposition of mass eigenstates;

$`\nu _e`$ $`=`$ $`\mathrm{cos}\mathrm{\Theta }\nu _1+\mathrm{sin}\mathrm{\Theta }\nu _2`$ (28)
$`\nu _\mu `$ $`=`$ $`\mathrm{sin}\mathrm{\Theta }\nu _1+\mathrm{cos}\mathrm{\Theta }\nu _2`$ (29)
where $`\nu _e`$ and $`\nu _\mu `$ denotes electron and muon states, respectively, and $`\nu _1`$ and $`\nu _2`$ mass eigenstates, and $`\mathrm{\Theta }`$ is a mixing angle. It should be clear that this scheme only works if the neutrino masses are different from each other, and therefore in general are non-zero; this means that in general there are right-handed neutrinos as well as left-handed ones. The right-handed neutrinos, however, interact with matter only very weakly. That is why they have not yet been observed.

Now how does the interference work? The wave packets which represent electron and muon neutrinos are constituted by the sum of the wave packets which represent mass eigenstates. These wave packets will travel with different speeds, hence, even though they may be centered at the same point in space at the point of production, there will be a seperation between the packets at the point of detection. In order for interference to be observed, however, the spatial separation between them must not be too large. The waves will then interfere at the detector and we can observe the phase difference. For more detailed discussion see Refs..

The model which explains this phenomenon realistically is the radial motion analysis. We begin the investigation, however, by studying azimuthal motion. Although not realistic, this has the virtue of simplicity; and it has also been considered in Ref., though in a different formalism. Our method of procedure will be to write down the Dirac equation and find phases corresponding to mass eigenstates, then finally calculate phase differences.

### 3.1 The Azimuthal Motion

The Hamiltonian for the azimuthal motion, that is, $`\stackrel{}{p}=(p_r,p_\theta ,p_\phi )=(0,0,p)`$ could be written from eqn.(27) as

$`H=f\alpha ^3p_\phi c+f\beta mc^2+{\displaystyle \frac{3}{2}}\mathrm{}cfA^0\gamma _5+{\displaystyle \frac{3}{2}}\mathrm{}cf\stackrel{}{\mathrm{\Sigma }}.\stackrel{}{A}.`$ (30)
This Hamiltonian couples the positive energy states to the negative energy states because of $`\alpha ^3`$ and $`\gamma _5`$ :

$`H\psi =E\psi `$ (31)
where $`\psi `$ is a four-component spinor whose first two components correspond to positive energy and second two to negative energy states. Firstly we decouple these states (or block diagonalize the Hamiltonian) by writing

$`\mathrm{det}(\mathrm{H}\mathrm{E})=0`$ (32)
where $`E=\left(\begin{array}{cc}E_+& 0\\ 0& E_{}\end{array}\right)`$. Therefore the upper block $`E_+`$ and lower block $`E_{}`$ satisfy

$`E_+`$ $`=`$ $`fpc+{\displaystyle \frac{fm^2c^3}{2p}}+{\displaystyle \frac{3}{2}}\mathrm{}cf(A^0\sigma ^3+\stackrel{}{\sigma }.\stackrel{}{A})`$ (33)
$`E_{}`$ $`=`$ $`fpc{\displaystyle \frac{fm^2c^3}{2p}}+{\displaystyle \frac{3}{2}}\mathrm{}cf(A^0\sigma ^3+\stackrel{}{\sigma }.\stackrel{}{A}).`$ (34)
Here we assume that the torsional terms are small compared with the momentum and mass terms333 Numerical values are given at the end of next subsection . With a unitary transformation

$`\psi =U\xi `$ (35)
where U is a unitary matrix and $`\xi =\left(\begin{array}{c}\xi _+\\ \xi _{}\end{array}\right)`$, our equation system decouples into two equations

$`E_+\xi _+`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}\xi _+`$ (36)
$`E_{}\xi _{}`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}\xi _{}.`$ (37)
We will be interested in only the positive energy states for simplicity. Similar calculations may be done for the negative energy states. At the end these states will correspond to the mass eigenstates whose linear superpositions will define, for example, electron and muon neutrinos as given by (28) and (29). Now eqn.(36) can be cast easily as two separate equations by diagonalizing, one for spin up, one for spin down444 Spin operator is chosen as $`\stackrel{}{\mathrm{\Sigma }}`$. Besides, it is only in the MASSLESS limit that the neutrino helicity is $`\pm 1`$, so the state we are discussing is actually allowed, since we are discussing MASSIVE Dirac neutrinos. ;

$`\{fpc+{\displaystyle \frac{fm^2c^3}{2p}}+{\displaystyle \frac{3}{2}}\mathrm{}cfA\}\xi _+^{}`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}\xi _+^{}`$ (38)
$`\{fpc+{\displaystyle \frac{fm^2c^3}{2p}}{\displaystyle \frac{3}{2}}\mathrm{}cfA\}\xi _+^{}`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}\xi _+^{},`$ (39)
whose time evolutions are given by

$`\xi _+^{}(t)`$ $`=`$ $`e^{i\mathrm{\Phi }^{}(t)}\xi _+^{}(0)`$ (40)
$`\xi _+^{}(t)`$ $`=`$ $`e^{i\mathrm{\Phi }^{}}(t)\xi _+^{}(0)`$ (41)
where $`A=\sqrt{(A^0A^3)^2+(A^1)^2+(A^2)^2}`$. The phase of the spin up state is given by

$`\mathrm{\Phi }^{}(t)={\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \{fpc+\frac{fm^2c^3}{2p}+\frac{3}{2}\mathrm{}cfA\}𝑑t}.`$ (42)
For ultrarelativistic neutrinos $`pcE`$ and $`cdtRd\phi ,`$ so the above becomes

$`\mathrm{\Phi }^{}={\displaystyle \frac{1}{\mathrm{}}}\{{\displaystyle \frac{fER}{c}}\mathrm{\Delta }\phi +{\displaystyle \frac{fm^2c^3R}{2E}}\mathrm{\Delta }\phi +{\displaystyle \frac{3}{2}}\mathrm{}fAR\mathrm{\Delta }\phi \}.`$ (43)
Similarly the phase of the spin down state is

$`\mathrm{\Phi }^{}={\displaystyle \frac{1}{\mathrm{}}}\{{\displaystyle \frac{fER}{c}}\mathrm{\Delta }\phi +{\displaystyle \frac{fm^2c^3R}{2E}}\mathrm{\Delta }\phi {\displaystyle \frac{3}{2}}\mathrm{}fAR\mathrm{\Delta }\phi \}.`$ (44)
These phases alone do not have an absolute meaning; the quantities relevant for the interference pattern at the observation point are discussed below.

Now consider electron and muon neutrinos given by (28) and (29) above, where $`\nu _1`$ and $`\nu _2`$ correspond to $`\xi _+`$’s. If electron neutrinos are produced at $`t=0`$, the time evolution of $`\nu _e`$ is

$`\nu _e(t)`$ $`=`$ $`\mathrm{cos}\mathrm{\Theta }\nu _1(t)+\mathrm{sin}\mathrm{\Theta }\nu _2(t)`$ (45)
$`=`$ $`\mathrm{cos}\mathrm{\Theta }e^{i\mathrm{\Phi }_1(t)}\nu _1+\mathrm{sin}\mathrm{\Theta }e^{i\mathrm{\Phi }_2(t)}\nu _2.`$
The probability of measuring muon neutrino at a later time is given by

$`P(\nu _e\nu _\mu )=\mathrm{sin}^22\mathrm{\Theta }\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Phi }}{2}}`$ (46)
where $`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2\mathrm{\Phi }_1`$. As is seen from (46) the phase difference has an absolute meaning. Now there are four possibilities.

* Both mass eigenstates are spin up, in which case eqn.(45) becomes

$`\nu _e(t)=\mathrm{cos}\mathrm{\Theta }e^{i\mathrm{\Phi }_1^{}(t)}\nu _1+\mathrm{sin}\mathrm{\Theta }e^{i\mathrm{\Phi }_2^{}(t)}\nu _2.`$ (47)
From (43) and (44), the phase difference is

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}={\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2\frac{E}{f}\mathrm{}}}R\mathrm{\Delta }\phi .`$ (48)
* Both mass eigenstates are spin down. In this case an analogous equation to (47) holds and the phase difference $`\mathrm{\Delta }\mathrm{\Phi }`$ is

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}={\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2\frac{E}{f}\mathrm{}}}R\mathrm{\Delta }\phi .`$ (49)
* The first mass eigenstate is spin up, the second spin down, then

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}=\{{\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2\frac{E}{f}\mathrm{}}}3fA\}R\mathrm{\Delta }\phi .`$ (50)
* The first one is spin down, the second spin up, then

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}=\{{\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2\frac{E}{f}\mathrm{}}}+3fA\}R\mathrm{\Delta }\phi `$ (51)

where $`\mathrm{\Delta }m^2=m_2^2m_1^2`$. From (48) and (49) it is seen that if both mass eigenstates have the same spin projection there is no contribution to the oscillation coming from the torsion, but if the mass eigenstates have opposite spin it is seen from (50) and (51) that there is a contribution to the neutrino oscillation coming from the torsion.

When we do an order of magnitude analysis we see that the contribution coming from the torsion is really very small compared with that coming from the mass difference555 Numerical values are given at the end of the section 3.2. Here we also emphasise that $`\frac{E}{f}`$ is the locally measured energy . If there is no torsion, i.e., $`A^a=0`$, our result becomes the same as the result of Ref., and in this case there is no gravitational contribution to the phase difference, since the locally measured energy is $`\frac{E}{f}`$ rather than $`E`$. We show below that in the case of radial motion the situation is different.

### 3.2 The Radial Motion

The Hamiltonian for the radial motion, that is, $`\stackrel{}{p}=(p,0,0)`$ may be written from eqn.(27) as follows

$`H=f^2\alpha ^1p_rc+f\beta mc^2+{\displaystyle \frac{3}{2}}\mathrm{}cfA^0\gamma _5+{\displaystyle \frac{3}{2}}\mathrm{}cf\stackrel{}{\mathrm{\Sigma }}.\stackrel{}{A}.`$ (52)
Applying the same arguments as in the previous subsection we find the upper and lower blocks of the Hamiltonian

$`E_+`$ $`=`$ $`f^2pc+{\displaystyle \frac{m^2c^3}{2p}}+{\displaystyle \frac{3}{2}}\mathrm{}cf(A^0\sigma ^1+\stackrel{}{\sigma }.\stackrel{}{A})`$ (53)
$`E_{}`$ $`=`$ $`f^2pc{\displaystyle \frac{m^2c^3}{2p}}{\displaystyle \frac{3}{2}}\mathrm{}cf(A^0\sigma ^1+\stackrel{}{\sigma }.\stackrel{}{A}).`$ (54)
As before, we are only interested in the positive energy states which obey

$`E_+\xi _+`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{}{t}}\xi _+.`$ (55)
Again we first diagonalize eqn.(55) in order to find eigenvalues of spin up and down states. In this case the phases appropriate to the spin up and spin down particles are, respectively,

$`\mathrm{\Phi }^{}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}\{{\displaystyle \frac{E}{c}}\mathrm{\Delta }r{\displaystyle \frac{2MGE}{c^3}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}}+{\displaystyle \frac{m^2c^3}{2E}}\mathrm{\Delta }r`$ (56)
$`+{\displaystyle \frac{3}{2}}\mathrm{}A(\mathrm{\Delta }r{\displaystyle \frac{MG}{c^2}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}})\}`$
$`\mathrm{\Phi }^{}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}\{{\displaystyle \frac{E}{c}}\mathrm{\Delta }r{\displaystyle \frac{2MGE}{c^3}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}}+{\displaystyle \frac{m^2c^3}{2E}}\mathrm{\Delta }r`$ (57)
$`{\displaystyle \frac{3}{2}}\mathrm{}A(\mathrm{\Delta }r{\displaystyle \frac{MG}{c^2}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}})\}`$
where $`A=\sqrt{(A^0A^1)^2+(A^2)^2+(A^3)^2}`$. Here we made the assumptions $`pcE`$, $`cdtdr`$, $`f=(1\frac{2MG}{rc^2})^{1/2}1\frac{MG}{rc^2}`$, $`\mathrm{\Delta }r=r_Br_A`$, where A and B refer respectively to the points of production and detection; hence $`r_A`$ is the radius of the sun and $`r_B`$ is the distance from the center of the sun to the surface of the earth.

If we apply the same procedure as in the previous subsection, from eqn.(56) and (57) there are three different phase differences;

* both mass eigenstates have the same spin state

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}={\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2E\mathrm{}}}\mathrm{\Delta }r`$ (58)
* the first is spin up, the second spin down

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}={\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2E\mathrm{}}}\mathrm{\Delta }r3A(\mathrm{\Delta }r{\displaystyle \frac{MG}{c^2}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}})`$ (59)
* the first is down, the second up

$`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_2^{}\mathrm{\Phi }_1^{}={\displaystyle \frac{\mathrm{\Delta }m^2c^3}{2E\mathrm{}}}\mathrm{\Delta }r+3A(\mathrm{\Delta }r{\displaystyle \frac{MG}{c^2}}\mathrm{ln}{\displaystyle \frac{r_B}{r_A}})`$ (60)

where the mass difference is $`\mathrm{\Delta }m^2=m_2^2m_1^2`$.

If an order of magnitude analysis is done, for vacuum $`\mathrm{\Delta }m^2c^410^{10}eV^2`$, $`E10MeV`$, $`r_AR_{}=7\times 10^8m`$, the earth-sun distance $`\mathrm{\Delta }r1.5\times 10^{11}m`$ and $`r_B=\mathrm{\Delta }r+r_A1.5\times 10^{11}m`$ and for the sun $`\frac{MG}{c^2}1.5km`$. As for the magnitude of the vector $`A^a`$, it follows from (15) that

$`A^a\mathrm{}^2(\overline{\psi }\psi ).`$ (61)
Then $`A^a10^{62}m^1`$ and hence the torsional term in (59), (60) is much smaller than the mass difference term. When we put $`A^a=0`$ we recover the same result as Ref., but we see that in general there is a contribution to neutrino oscillation coming from the torsion of spacetime.

## 4 Result

We have treated the gravitational neutrino oscillation problem in a different way from previous authors by calculating the Hamiltonian of a Dirac particle in Riemann-Cartan spacetime and finding the dynamical phase. We began for simplicity by treating the (unrealistic) case of azimuthal motion; this gave us the opportunity to compare with , in which azimuthal motion is considered in Riemann spacetime. We have found that there is a torsional contribution which depends on the spin directions of mass eigenstates, eqn.(50) and (51). If we put the torsion equal to zero in these equations we recover the same result as Ref.. Next we studied the radial motion which is a more realistic model for our problem. We found here also a torsional contribution to the neutrino oscillation that is dependent on the spin polarizations of the mass eigenstates. Again, when the torsion vanishes we recover exactly the same result as in ref.

Acknowledgement

One of the authors (M.A.) thanks TUBITAK (Scientific and Technical Research Council of Turkey) for a grant through BAYG-BDP that makes his stay at UK possible and the Department of Physics, University of Kent at Canterbury for hospitality.