text:
# Split of zero-bias conductance peak in normal-metal / 𝒅 wave superconductor junctions

## I introduction

The zero-energy state (ZES) hu  formed at surfaces of superconductors is a consequence of the unconventional symmetry of Cooper pairs. Since the ZES appears just on the Fermi energy, it drastically affects transport properties through the interface of junctions consist of unconventional superconductors. RPP  For instance in normal-metal / high-$`T_c`$ superconductor junctions, a large peak is observed in the differential conductance at the zero bias voltage. tanaka0 ; kashiwaya ; Alff ; Wang ; wei ; iguchi ; geerk ; mao  The ZES is also responsible for the low-temperature anomaly of the Josephson current between the two unconventional superconductors. barash ; tanaka1 ; tanaka2 ; tanaka3 ; T1 ; asano ; arie ; ilichev1 ; ilichev2 ; shirai ; testa

An electron incident into a normal-metal / superconductor (NS) interface suffers the Andreev reflection andreev  by the pair potential in the superconductor. As a result, a hole traces back the original propagation path of the incident electron. This is called the retro property of a quasiparticle which supports the formation of the ZES. Strictly speaking, the electron-hole pairs just on the Fermi energy hold the retro property in the presence of the time reversal symmetry (TRS). Thus the ZES is sensitive to the TRS of the system. Actually, the zero-bias conductance peak (ZBCP) in NS junctions splits into two peaks under magnetic fields. fogelstrom ; YT021 ; YT022 ; YT023  The peak splitting is also discussed kashi95 ; matsumoto2 ; TJ1 ; TJ2 ; Tanuma2001 ; lubimova ; kitaura  when the broken time reversal symmetry state (BTRSS) is formed at the interface. Theoretical studies showed that such BTRSS’s are characterized by the $`s+id_{xy}`$ matsumoto2  or $`d_{xy}+id_{x^2y^2}`$ laughlin  wave pairing symmetry. Experimental results, however, are still controversial. Some experiments reported the split of the ZBCP at the zero magnetic field, covington ; biswas ; dagan ; sharoni ; kohen ; greene ; aprili2 ; krupke  other did not observe the splitting. Ekin ; Alff ; wei ; iguchi ; Sawa1 ; Aubin  The ZBCP is also sensitive to the exchange potential in ferromagnets attaching to unconventional superconductors. H2 ; Y2

In previous papers, we numerically showed that random potentials at the NS interface cause the split of the ZBCP at the zero magnetic field by using the recursive Green function method. asano2 ; asano3  We also showed that the splitting due to the impurity scattering can be seen more clearly when realistic electronic structures of high-$`T_c`$ materials are taken into account. asano4  Unfortunately, we could not make clear a mechanism of splitting. Our conclusion, however, contradicts to those of a number of theories barash2 ; golubov ; poenicke ; yamada ; tanaka5 ; luck  based on the quasiclassical Green function method. eilenberger ; larkin ; zaitsev ; shelankov ; bruder  The drastic suppression of the ZBCP by the interfacial randomness is the common conclusion of all the theories. The theories of the quasiclassical Green function method, however, concluded that the random potentials do not split the ZBCP. Thus this issue has not been fixed yet. There are mainly two reasons for the disagreement in the two theoretical approaches (i.e., the recursive Green function method and the quasiclassical Green function method). One is the treatment of the random potentials, the other is the effects of the rapidly oscillating wave functions on the conductance. In our simulations, we calculate the conductance without any approximation to the random potentials and the wave functions; this is an advantage of the recursive Green function method. lee ; asano3

In this paper, we discuss effects of the impurity scattering on the conductance in normal-metal/$`d`$ wave superconductor junctions by using the Lippmann-Schwinger equation. We assume that impurities are near the NS interface on the superconductor side. The differential conductance is analytically calculated within the single-site approximation based on the conductance formula. blonder ; takane  The split of the ZBCP due to the impurity scattering is the main conclusion of this paper. The impurity scattering affects the conductance in two ways: (i) drastically suppressing the conductance around the zero bias voltage and (ii) making the conductance peak wider. The split of the ZBCP is a consequence of the interplay between the two effects. In the present theory, we successfully explain typical conductance shapes observed in several experiments. We also show that the splitting peaks are merged into a single conductance peak for sufficiently high temperatures and that the peak splitting width increases with increasing external magnetic fields.

This paper is organized as follows. In Sec. II, we derive the reflection coefficients in NS junctions within the single-site approximation based on the Lippmann-Schwinger equation. The split of the ZBCP is discussed in Sec. III. In Sec. IV, calculated results are compared with experiments and another theories. In Sec. V, we summarize this paper.

## II Lippmann-Schwinger equation

Let us consider two-dimensional NS junctions as shown in Fig. 1, where normal metals ($`x<0`$) and $`d`$ wave superconductors ($`x>0`$) are separated by the potential barrier $`V_B(𝒓)=V_0\delta (x)`$. We assume the periodic boundary condition in the $`y`$ direction and the width of the junction is $`W`$.

The $`a`$ axis of high-$`T_c`$ superconductors is oriented by the 45 degrees from the interface normal. The pair potential of a high-$`T_c`$ superconductor is described by

$$\mathrm{\Delta }_𝒌=2\mathrm{\Delta }_0\overline{k}_x\overline{k}_y,$$
(1)
in the momentum space, where $`\mathrm{\Delta }_0`$ is the amplitude of the pair potential at the zero temperature, $`\overline{k}_x=\mathrm{cos}\gamma =k_x/k_F`$ and $`\overline{k}_y=\mathrm{sin}\gamma =k_y/k_F`$ are the normalized wave number on the Fermi surface in the $`x`$ and the $`y`$ direction, respectively. The Fermi wave number $`k_F`$ satisfies $`\mathrm{}^2k_F^2/2m=\mu _F`$, where $`\mu _F`$ is the Fermi energy. The schematic figure of the pair potential is shown in Fig. 1. The NS junctions are described by the Bogoliubov-de Gennes equation, degennes

$`{\displaystyle 𝑑𝒓^{}}`$ $`\left(\begin{array}{cc}\delta (𝒓𝒓^{})h_0(𝒓^{})& \mathrm{\Delta }(𝒓,𝒓^{})\\ \mathrm{\Delta }^{}(𝒓,𝒓^{})& \delta (𝒓𝒓^{})h_0(𝒓^{})\end{array}\right)`$ (4)
$`\times \left(\begin{array}{c}u(𝒓^{})\\ v(𝒓^{})\end{array}\right)=E\left(\begin{array}{c}u(𝒓)\\ v(𝒓)\end{array}\right),`$ (9)
$`h_0(𝒓)=`$ $`{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V_{pot}(𝒓)\mu _F,`$ (10)
$`V_{pot}(𝒓)=`$ $`V_B(𝒓)+V_I(𝒓),`$ (11)
$`\mathrm{\Delta }(𝑹_𝒄,𝒓_r)=`$ $`\{\begin{array}{ccc}\frac{1}{V_{vol}}_𝒌\mathrm{\Delta }_𝒌e^{i𝒌𝒓_r}& :& X_c>0\\ 0& :& X_c<0\end{array},`$ (14)
where $`𝑹_c=(X_c,Y_c)=(𝒓+𝒓^{})/2`$ and $`𝒓_r=𝒓𝒓^{}`$. Throughout this paper, we neglect the spatial dependence of the pair potential near the junction interface. This is a reasonable approximation when we consider the conductance around the zero bias voltage. tanaka4  The spatial dependence of the pair potential should be determined in a self-consistent way when we discussed the conductance far from the zero bias such as $`eV\mathrm{\Delta }_0`$. Here $`V`$ is the bias voltage applied to junctions.

We consider impurities near the interface on the superconductor side as indicated by crosses in Fig. 1. The potential of impurities is given by

$$V_I(𝒓)=V_i\underset{j=1}{\overset{N_i}{}}\delta (𝒓𝒓_j),$$
(15)
where $`N_i`$ is the number of impurities. In the absence of impurities, the transmission and the reflection coefficients are calculated from boundary conditions of wave functions at the junction interface as shown in Appendix A. By using these coefficients, four retarded Green functions are obtained as shown in Appendix B. The normal conductance of the junction is given by

$`G_N=`$ $`{\displaystyle \frac{2e^2}{h}}N_cT_B,`$ (16)
$`T_B=`$ $`{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{\mathrm{cos}^3\gamma }{z_0^2+\mathrm{cos}^2\gamma }},`$ (17)
where $`T_B`$ is the transmission probability of the junction, $`N_c=2W/\lambda _F`$ is the number of the propagating channels on the Fermi surface, $`\lambda _F=2\pi /k_F`$ is the Fermi wave length and $`z_0=mV_0/(\mathrm{}^2k_F)`$ represents the strength of the potential barrier at the NS interface. In the limit of $`z_0^21`$, $`T_B`$ is proportional to $`1/z_0^2`$.

Effects of impurities on the wave functions are taken into account by using the Lippmann-Schwinger equation,

$`\psi ^{(l)}(𝒓)=`$ $`\psi _0^{(l)}(𝒓)+{\displaystyle 𝑑𝒓^{}\widehat{G}_0(𝒓,𝒓^{})V_I(𝒓^{})\widehat{\sigma }_3\psi ^{(l)}(𝒓^{})},`$ (18)
$`=`$ $`\psi _0^{(l)}(𝒓)+{\displaystyle \underset{j=1}{\overset{N_i}{}}}\widehat{G}_0(𝒓,𝒓_j)V_i\widehat{\sigma }_3\psi ^{(l)}(𝒓_j),`$ (19)
where $`l`$ indicates a propagating channel characterized by the transverse wave number $`k_y^{(l)}`$. Here $`\psi _0^{(l)}(𝒓)`$ is the wave function in which an electronlike quasiparticle with $`k_y^{(l)}`$ is incident into the clean NS interface from normal metals and is described as

$`\psi _0^{(l)}(𝒓)=`$ $`\chi _l(y)[\left(\begin{array}{c}1\\ 0\end{array}\right)e^{iq_l^+x}+\left(\begin{array}{c}0\\ 1\end{array}\right)e^{iq_l^{}x}r_{NN}^{he}(l)`$ (24)
$`+\left(\begin{array}{c}1\\ 0\end{array}\right)e^{iq_l^+x}r_{NN}^{ee}(l)],`$ (27)
$`\chi _l(y)=`$ $`{\displaystyle \frac{e^{ik_y^{(l)}y}}{\sqrt{W}}},`$ (28)
for $`x<0`$, where $`q_l^\pm =\sqrt{k_l^2\pm k_F^2E/\mu _F}`$ is the wave number of a quasiparticle in normal metals and $`k_l^2+k_{y}^{(l)}{}_{}{}^{2}=k_F^2`$. For $`x>0`$, the wave function in clean junctions is given by

$`\psi _0^{(l)}(𝒓)=`$ $`\chi _l(y)\widehat{\mathrm{\Phi }}[\left(\begin{array}{c}u_l\\ v_l\end{array}\right)e^{ik_l^+x}t_{SN}^{ee}(l)`$ (31)
$`+\left(\begin{array}{c}v_l\\ u_l\end{array}\right)e^{ik_l^{}x}t_{SN}^{he}(l)],`$ (34)
$`u_l=`$ $`\sqrt{{\displaystyle \frac{E+\mathrm{\Omega }_l}{2E}}},`$ (35)
$`v_l=`$ $`\text{sgn}(k_y^{(l)})\sqrt{{\displaystyle \frac{E\mathrm{\Omega }_l}{2E}}},`$ (36)
$`\widehat{\mathrm{\Phi }}=`$ $`\left(\begin{array}{cc}e^{i\phi }& 0\\ 0& e^{i\phi }\end{array}\right),`$ (39)
where $`\phi `$ is a macroscopic phase of a superconductor, $`k_l^\pm =\left(k_l^2\pm k_F^2\mathrm{\Omega }_l/\mu _F\right)^{1/2}`$ is the wave number of a quasiparticle in superconductors, $`\mathrm{\Omega }_l=\sqrt{E^2\mathrm{\Delta }_l^2}`$ and $`\mathrm{\Delta }_l=2\mathrm{\Delta }_0\overline{k}_l\overline{k}_y^{(l)}`$. The wave function at an impurity $`\psi ^{(l)}(𝒓_j^{})`$ can be obtaind by $`𝒓𝒓_j^{}`$ in Eq. (19)

$`\psi _0^{(l)}(𝒓_j^{})=`$ $`{\displaystyle \underset{j=1}{\overset{N_i}{}}}\left[\widehat{\sigma }_0\delta _{j,j^{}}\widehat{G}_0^{SS}(𝒓_j^{},𝒓_j)V_i\widehat{\sigma }_3\right]\psi ^{(l)}(𝒓_j).`$ (40)
It is possible to calculate the exact conductance if we obtain $`\psi ^{(l)}(𝒓_j)`$ for all impurities by solving Eq. (40). Actually it was confirmed that the conductance calculated from the numerical solution of Eq. (40) is exactly identical to that computed in another numerical methods such as the recursive Green function method. nonoyama ; lee  In this paper, we solve Eq. (40) within the single-site approximation, where the multiple scattering effect involving many impurities (Anderson localization) are neglected. However the multiple scattering by an impurity is taken into account up to the infinite order of the scattering events. In the summation of $`j`$ in Eq. (40), only the contribution with $`j=j^{}`$ is taken into account in the single-site approximation. asano5  In this way, the wave function at $`𝒓_j`$ is approximately given by

$$\psi ^{(l)}(𝒓_j)\left[\widehat{\sigma }_0\widehat{G}_0^{SS}(𝒓_j,𝒓_j)V_i\widehat{\sigma }_3\right]^1\psi _0^{(l)}(𝒓_j).$$
(41)
We note that the single-site approximation yields the exact conductance when $`N_i=1`$. Within the single-site approximation, Eq. (19) can be solved as

$`\psi _{\text{SSA}}^{(l)}(𝒓)=`$ $`\psi _0^{(l)}(𝒓)+{\displaystyle \underset{j}{\overset{N_i}{}}}\widehat{G}_0^{NS}(𝒓,𝒓_j)V_i\widehat{\sigma }_3`$
$`\times \left[\widehat{\sigma }_0\widehat{G}_0^{SS}(𝒓_j,𝒓_j)V_i\widehat{\sigma }_3\right]^1\psi _0^{(l)}(𝒓_j),`$ (42)
for $`x<0`$. On the right hand side of Eq. (42), all functions have been given by analytical expressions.

In the presence of the impurity scattering, the wave function Eq. (42) can be expressed as

$`\psi _{\text{SSA}}^{(l)}(𝒓)=\left(\begin{array}{c}1\\ 0\end{array}\right)\chi _l(y)e^{iq_l^+x}`$ (45)
$`+`$ $`{\displaystyle \underset{l^{}}{}}\chi _l^{}(y)\left[\left(\begin{array}{c}0\\ 1\end{array}\right)e^{iq_l^{}^{}x}A_{l^{},l}+\left(\begin{array}{c}1\\ 0\end{array}\right)e^{iq_l^{}^+x}B_{l^{},l}\right],`$ (50)
for $`x<0`$, where $`A_{l^{},l}`$ and $`B_{l^{},l}`$ are the Andreev and the normal reflection coefficients in the presence of impurities, respectively. These coefficients are obtained from relations

$`{\displaystyle _{W/2}^{W/2}}`$ $`dy\chi _m^{}(y)(0,1)\psi _{\text{SSA}}^{(l)}(𝒓)=e^{iq_m^{}x}A_{m,l},`$ (51)
$`{\displaystyle _{W/2}^{W/2}}`$ $`dy\chi _m^{}(y)(1,0)\psi _{\text{SSA}}^{(l)}(𝒓)=e^{iq_l^+x}\delta _{l,m}+e^{iq_m^+x}B_{m,l}.`$ (52)
The scattering theory based on the Lippmann-Schwinger equation requires complicated algebra as shown below because the perturbation expansion is carried out in the real space. In return, effects of the impurity scattering can be taken into account up to the infinite order of the perturbation expansion without using any self-consistent treatments. In addition, the reflection coefficients are explicit functions of the impurity positions in a single disordered sample. These are advantages of the present method.

In what follows, we consider low transparent junctions, (i.e., $`z_0^21`$). From the reflection coefficients in Appendix A, the Green function in the superconductor is given by

$`\widehat{G}_0^{SS}`$ $`(𝒓,𝒓)=i\pi N_0{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}d\gamma [{\displaystyle \frac{E}{2\mathrm{\Omega }}}\widehat{\sigma }_0{\displaystyle \frac{z_0^2\mathrm{\Delta }_𝒌^2e^{2ipx}}{2\mathrm{\Xi }\mathrm{\Omega }}}\widehat{\sigma }_0+{\displaystyle \frac{z_0^2\mathrm{\Delta }_𝒌^2\mathrm{cos}(2k_Fx\mathrm{cos}\gamma )e^{2ipx}}{2\mathrm{\Xi }\mathrm{\Omega }}}\widehat{\sigma }_0`$
$`{\displaystyle \frac{z_0^2E}{2\mathrm{\Xi }}}\left\{{\displaystyle \frac{E}{\mathrm{\Omega }}}\mathrm{cos}(2k_Fx\mathrm{cos}\gamma )\widehat{\sigma }_0+i\mathrm{sin}(2k_Fx\mathrm{cos}\gamma )\widehat{\sigma }_3\right\}e^{2ipx}{\displaystyle \frac{\mathrm{\Delta }_𝒌^2\mathrm{cos}^2\gamma e^{2ipx}}{4\mathrm{\Omega }\mathrm{\Xi }}}\widehat{\sigma }_0`$
$`{\displaystyle \frac{z_0\mathrm{cos}\gamma }{2\mathrm{\Xi }}}e^{2ipx}\{E\mathrm{cos}(2k_Fx\mathrm{cos}\gamma )\widehat{\sigma }_3+i\mathrm{\Omega }\mathrm{sin}(2k_Fx\mathrm{cos}\gamma )\widehat{\sigma }_0\}],`$ (53)
where $`p\frac{k_F}{2\mathrm{cos}\gamma }\frac{\mathrm{\Omega }}{\mu _F}.`$ The local density of states tanuma98  at $`𝒓`$ is defined by

$$N_s(E,x)=\frac{1}{\pi }\text{Im}\text{Tr}\widehat{G}_0^{SS}(𝒓,𝒓).$$
(54)
The first term in Eq. (53) contributes to the bulk density of states. Since $`2p`$ is roughly estimated to be $`i/\xi _0`$, another terms contribute to the local density of states near the interface, where $`\xi _0=\mathrm{}v_F/\pi \mathrm{\Delta }_0`$ is the coherence length and $`v_F=\mathrm{}k_F/m`$ is the Fermi velocity. In low transparent junctions, 6th, 7th and 8th terms are negligible. The 4th and the 5th terms are also negligible because integrals of such rapidly oscillating functions become very small. The 2nd and 3rd terms are dominant for $`E\mathrm{\Delta }_0`$. In Fig. 2 (a), we show the trace of the Green function in Eq. (53) as a function of $`E`$, where $`z_0`$=10, $`\mathrm{\Delta }_0=0.1\mu _F`$, $`xk_F=6`$, and $`N_0=m/(\pi \mathrm{}^2)`$ is the normal density of states in the unit area. Since $`\xi _0k_F6.3`$, the results correspond to the Green function at a distance $`\xi _0`$ away from the interface. The horizontal axis in Fig. 2 (a) is normalized by

$$E_{\text{ZEP}}\frac{\mathrm{\Delta }_0}{z_0^2}.$$
(55)
The solid and the broken lines represent negative of the imaginary part and the real part of Eq. (53), respectively.

As shown in Fig. 2 (a), the imaginary part of the Green function has a large peak around $`E=0`$ reflecting the ZES formed at the junction interface. The energy scale, $`E_{\text{ZEP}}`$, characterizes the width of the zero-energy peak. The real part of the Green function first increases with decreasing $`E`$ then suddenly decreases to zero for $`E0`$. The detail analysis indicates that the real part of the Green function has its maximum around an energy

$$E_{\text{dip}}\frac{\mathrm{\Delta }_0}{z_0^2(xk_F)^3}.$$
(56)
The Green function for $`E<E_{\text{ZEP}}`$ is approximately calculated from the 2nd and the 3rd terms of Eq. (53)

$`\widehat{G}_0^{SS}`$ $`(𝒓,𝒓)2\pi N_0z_0^2e^{x/\xi _0}(g_2ig_1)\widehat{\sigma }_0,`$ (57)
$`g_1=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{\mathrm{\Delta }_0^2\mathrm{cos}^4\gamma \mathrm{sin}^2\gamma \mathrm{sin}^2(xk_F\mathrm{cos}\gamma )}{E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }},`$ (58)
$`g_2=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{Ez_0^2\mathrm{\Delta }_0\mathrm{cos}\gamma \mathrm{sin}\gamma \mathrm{sin}^2(xk_F\mathrm{cos}\gamma )}{E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }},`$ (59)
where we use $`\mathrm{\Omega }i|2\mathrm{\Delta }_0\mathrm{cos}\gamma \mathrm{sin}\gamma |`$. The imaginary part of the Green function, $`g_1`$, is of the order of unity when $`EE_{\text{ZEP}}`$. However, $`g_1`$ at $`E=0`$ becomes much larger than unity for $`xk_F1`$ because

$$g_1(E=0)=\frac{2}{\pi }_0^{\pi /2}𝑑\gamma \frac{\mathrm{sin}^2\left\{xk_F\mathrm{cos}\gamma \right\}}{\mathrm{cos}^2\gamma }xk_F.$$
(60)
Thus the energy scale, $`E_{\text{dip}}`$, characterizes the drastic increase of $`g_1`$ and the drastic decrease of $`g_2`$. The local density of states at $`E=0`$ calculated from Eqs. (53) and (54) is plotted as a function of $`xk_F`$ in Fig. 2(b) with the solid line. For comparison, we also show the analytical results represented by

$$\text{Im}\widehat{G}_0^{SS}(𝒓,𝒓)|_{E=0}2\pi N_0z_0^2e^{x/\xi _0}xk_F\widehat{\sigma }_0,$$
(61)
with the broken line. The results show the remarkable enhancement of the local density of states around $`x\xi _0`$. This implies that the ZES is formed around $`x\xi _0`$.

Here we note several remarks as follows. To calculate the Green function, we consider the 3rd term in Eq.(53) which rapidly oscillates as $`\mathrm{cos}(2xk_F\mathrm{cos}\gamma )`$. Such rapidly oscillating terms are usually neglected in the quasiclassical Green function method. We, however, cannot neglect the 3rd term because it removes the divergence of the local density of state at $`E=0`$tanaka6 ; matsumoto  The 3rd term also becomes important when we calculate the local density of states just at the surface, (i.e., $`x=0`$),

$`{\displaystyle \frac{N_s(E,0)}{N_0}}=`$ $`\text{Re}{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma \left[{\displaystyle \frac{E}{\mathrm{\Omega }}}{\displaystyle \frac{2E^2z_0^2+\mathrm{\Delta }_𝒌^2\mathrm{cos}^2\gamma }{2\mathrm{\Xi }\mathrm{\Omega }}}\right],`$ (62)
$`{\displaystyle \frac{2}{\pi }}K\left({\displaystyle \frac{\mathrm{\Delta }_0}{E}}\right)`$ $`+{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{E^2z_0^2+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }{E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }},`$ (63)
where $`K(x)`$ is the complete elliptic integral of the first kind and describes the bulk density of states. Another terms come from the 4th and 6th terms in Eq. (53). The first equation is the exact expression and we use $`E<E_{\text{ZEP}}`$ in the second line. We exactly obtain $`N(E=0,x=0)=N_0`$. Thus there is no remarkable enhancement in the zero energy local density of states just at the interface. The 2nd and the 3rd terms in Eq. (53) do not contribute to the $`N_s(E,0)`$ because they exactly cancel with each other at $`x=0`$.

In the next section, the conductance for $`E<E_{\text{ZEP}}`$ will be discussed. It is possible to rewrite a part of Eq. (42) as

$`\left[\widehat{\sigma }_0V_i\widehat{G}_0^{SS}(𝒓,𝒓)\widehat{\sigma }_3\right]^1={\displaystyle \frac{\widehat{\sigma }_0+s\widehat{\sigma }_3}{1s^2}},`$ (64)
$`s_1`$ $`=Q_0g_1,`$ (65)
$`s_2`$ $`=Q_0g_2,`$ (66)
$`Q_0`$ $`=2\pi V_iN_0z_0^2e^{x/\xi _0},`$ (67)
$`s`$ $`=s_2is_1,`$ (68)
where $`s`$ corresponds to the self-energy of a quasiparticle scattered by an impurity once and describes two important features. First, $`s`$ becomes large in low transparent junctions even if $`V_iN_0`$ is fixed at a small constant because $`Q_0`$ in Eq.(67) is proportional to $`V_iN_0z_0^2V_iN_0/T_B`$. Thus $`Q_0`$ represents the normalized strength of the impurity scattering. This behavior explains the previous numerical simulation. asano2  Second, effects of impurity scattering far away from the interface on the conductance is negligible because $`Q_0`$ decreases exponentially with the increase of $`x`$. Impurities around the ZES, (i.e., $`x\xi _0`$) seriously affect the conductance for $`E<E_{\text{ZEP}}`$.

## III conductance

The differential conductance in NS junctions is calculated from the normal and the Andreev reflection coefficients, blonder ; takane

$`G_{NS}`$ $`(eV)={\displaystyle \frac{2e^2}{h}}{\displaystyle \underset{l,m}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\left({\displaystyle \frac{f_{\text{FD}}(EeV)}{(eV)}}\right)`$
$`\times \left[\delta _{l,m}|\overline{B}_{m,l}(E)|^2+|\overline{A}_{m,l}(E)|^2\right],`$ (69)
$`\overline{A}_{m,l}=`$ $`\sqrt{{\displaystyle \frac{k_m}{k_l}}}A_{m,l},`$ (70)
$`\overline{B}_{m,l}=`$ $`\sqrt{{\displaystyle \frac{k_m}{k_l}}}B_{m,l},`$ (71)
where $`f_{\text{FD}}(E)`$ is the Fermi-Dirac distribution function. When $`z_0^21`$ and $`E<E_{\text{ZEP}}`$, the reflection coefficients are calculated as

$`\overline{A}_{m,l}=`$ $`\delta _{l,m}r_{NN}^{he}(l)`$
$`+`$ $`e^{i\phi }L_{m,l}\left[\mathrm{\Delta }_m|\mathrm{\Delta }_l|(s+1)+\mathrm{\Delta }_l|\mathrm{\Delta }_m|(s1)\right],`$ (72)
$`\overline{B}_{m,l}=`$ $`\delta _{l,m}r_{NN}^{ee}(l)`$
$`+iL_{m,l}\left[|\mathrm{\Delta }_m||\mathrm{\Delta }_l|(s+1)+\mathrm{\Delta }_l\mathrm{\Delta }_m(s1)\right],`$ (73)
$`L_{m,l}=`$ $`{\displaystyle \frac{\pi N_0V_iz_0^2}{(1s^2)k_F}}{\displaystyle \underset{j=1}{\overset{N_i}{}}}\chi _l(y_j)\chi _m^{}(y_j)`$
$`\times {\displaystyle \frac{\sqrt{\overline{k}_l\overline{k}_m}}{\mathrm{\Xi }_l\mathrm{\Xi }_m}}e^{i(p_m+p_l)x_j}\mathrm{sin}(k_mx_j)\mathrm{sin}(k_lx_j).`$ (74)
The conductance is then given by

$`G_{NS}={\displaystyle \frac{2e^2}{h}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\left({\displaystyle \frac{f_{\text{FD}}(EeV)}{(eV)}}\right)`$
$`\times \left[N_cg^{(0)}4{\displaystyle \underset{j=1}{\overset{N_i}{}}}\mathrm{\Gamma }_j\right],`$ (75)
$`g^{(0)}=2{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{\mathrm{\Delta }_0^2\mathrm{cos}^7\gamma \mathrm{sin}^2\gamma }{E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }},`$ (76)
$`\mathrm{\Gamma }_j=`$ $`\text{Re}{\displaystyle \frac{sQ_0(2I_2iI_1+iI_3)}{s^21}},`$ (77)
$`I_1=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{\mathrm{\Delta }_0^4\mathrm{cos}^{10}\gamma \mathrm{sin}^4\gamma \mathrm{sin}^2(x_jk_F\mathrm{cos}\gamma )}{(E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma )^2}},`$ (78)
$`I_2=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{Ez_0^2\mathrm{\Delta }_0^3\mathrm{cos}^7\gamma \mathrm{sin}^3\gamma \mathrm{sin}^2(x_jk_F\mathrm{cos}\gamma )}{(E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma )^2}},`$ (79)
$`I_3=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\gamma {\displaystyle \frac{E^2z_0^4\mathrm{\Delta }_0^2\mathrm{cos}^4\gamma \mathrm{sin}^2\gamma \mathrm{sin}^2(x_jk_F\mathrm{cos}\gamma )}{(E^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma )^2}}.`$ (80)
The first term of Eq. (75), $`N_cg^{(0)}`$, is the conductance in clean junctions and $`\mathrm{\Gamma }_j`$ represents effects of the impurity scattering on the conductance. When we calculate $`|\overline{A}_{m,l}|^2`$ and $`|\overline{B}_{m,l}|^2`$, the summation with respect to impurities $`_j^{N_i}_j^{}^{N_i}`$ must be carried out only for $`j^{}=j`$ in the single-site approximation. asano5  As a consequence, the current conservation law is satisfied for $`E<E_{\text{ZEP}}`$.

To study effects of impurities on the conductance, we first assume $`x_j=x_0`$ for all impurities. The conductance is rewritten as

$`G_{NS}=`$ $`{\displaystyle \frac{2e^2}{h}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\left({\displaystyle \frac{f_{\text{FD}}(EeV)}{(eV)}}\right)`$
$`\times N_c\left[g^{(0)}2n_i\mathrm{\Gamma }_j\right],`$ (81)
where $`n_i=N_i\lambda _F/W`$ is the dimensionless line density of impurities less than unity. When scattering effects are strong, $`|s|1`$, $`N_i`$ cannot be much larger than $`W/\lambda _F`$. This limits the applicability of the single-site approximation.

We show conductance for several choices of $`x_0k_F`$ and $`V_iN_0`$ in Fig. 3, where $`z_0=10`$ and $`n_i=0.9`$. The two parameters are chosen as $`x_0k_F=10`$, $`V_iN_0=0.01`$ in (a), $`x_0k_F=2.0`$, $`V_iN_0=0.005`$ in (b), $`x_0k_F=26.0`$, $`V_iN_0=0.1`$ in (c) and $`x_0k_F=12.0`$, $`V_iN_0=0.1`$ in (d). The broken line is the conductance in clean junctions. The temperature is fixed at a very low temperature $`T=0.01E_{\text{ZEP}}`$ which is estimated to be 0.05 K by using $`\mathrm{\Delta }_0=50`$ meV for $`z_0=10`$. As shown in (a)-(c), the ZBCP is splitting into two peaks by the impurity scattering. While the results in (d) shows the single ZBCP.

Roughly speaking, the impurity scattering affects the ZBCP in two ways: (i) it decreases the conductance around the zero bias voltage and (ii) it makes the ZBCP wider. The two effects (i) and (ii) are well characterized by the $`E`$ dependence of $`\mathrm{\Gamma }_j`$ and the sign change of $`\mathrm{\Gamma }_j`$ in Eq. (81), respectively. In Fig. 4, $`\mathrm{\Gamma }_j`$ is plotted as a function of $`E`$ for several $`x_0k_F`$, where $`V_iN_0`$ is fixed at 0.01. We note that $`\mathrm{\Gamma }_j`$ for $`x_0k_F=10`$ yields the conductance in Fig. 3 (a). When $`x_0k_F\xi _0`$, $`\mathrm{\Gamma }_j0`$ and impurity scattering is negligible as shown for $`x_0k_F=50`$ because $`Q_0`$ in Eq. (67) becomes almost zero. For $`x_0k_F=`$ 2, 5 and 10, $`\mathrm{\Gamma }_j`$ increases with decreasing $`E`$, which indicates the enhancement of the impurity scattering around $`E=0`$. The suppression of the conductance around the zero bias is explained in terms of the drastic increase of the local density of states with decreasing $`E`$ as shown in Fig. 2 (a). Therefore the suppression of the zero-bias conductance happens irrespective of $`x_0k_F`$ and $`V_iN_0`$. In addition, a nonmonotonic $`E`$ dependence of $`\mathrm{\Gamma }_j`$ for $`x_0k_F=10`$ is a source of small conductance peak at the zero-bias in Fig. 3 (a). The same small peak is also found in Fig. 3 (b).

The widening of the ZBCP can be explained by the sign change of $`\mathrm{\Gamma }_j`$. When $`E>0.15E_{\text{ZEP}}`$ in Fig. 4, $`\mathrm{\Gamma }_j`$ for $`x_0k_F=10`$ becomes negative and impurities enhance the conductance. As a consequence, the conductance peak becomes wider than $`g^{(0)}`$. The split of the ZBCP is a consequence of the interplay between the suppression of the conductance around the zero bias voltage and the widening of the ZBCP as shown in Fig. 3 (a)-(c). Therefore the sign change of $`\mathrm{\Gamma }_j`$ explains the split of the ZBCP.

For $`E_{\text{dip}}<E<E_{\text{ZEP}}`$, $`\mathrm{\Gamma }_j`$ is almost a decreasing function of $`E`$ and is positive at $`E=E_{\text{dip}}`$ as shown in Fig. 4. It is followed from Eq. (77) that

$$\mathrm{\Gamma }_j|s|^2\left\{s_1(I_1I_3)+2I_2s_2\right\}+\left\{s_1(I_1I_3)2I_2s_2\right\}.$$
(82)
Within our study, $`s_1(I_1I_3)`$ tends to be much smaller than $`2I_2s_2`$ for $`E_{\text{dip}}<E<E_{\text{ZEP}}`$, which implies an importance of the real part of the self-energy, $`s_2`$, for the splitting. The sign change of $`\mathrm{\Gamma }_j`$ happens when the impurity scattering is sufficiently weak so that

$$|s|^2=s_1^2+s_2^21,$$
(83)
is satisfied. When $`|s(E=E_{\text{dip}})|`$ is a small value less than unity, the effects of impurities are negligible and the conductance almost remains unchanged from that in clean junctions. The split also cannot be seen when $`|s(E=E_{\text{ZEP}})|`$ is larger than unity. An example is shown in Fig. 3(d), where $`x_0k_F=12.0`$, $`V_iN_0=0.1`$ and $`|s(E=E_{\text{ZEP}})|`$ is estimated to be 1.5. In this case, the suppression of the zero-bias conductance dominates over the widening of the ZBCP. As a result, the conductance is always smaller than that in clean junctions and the ZBCP remains in a single peak.

We should pay attention to a similarity in the shapes of the conductance in the present theory and those in experiments. Amazingly, the conductance structure in Fig. 3 (a) is very similar to that observed in the experiment covington . It is possible to find a very small conductance peak at $`V=0`$ in addition to the splitting peaks around $`V\pm 1`$mV in Fig. 2 of Ref. covington, . The impurities away from the interface explains another conductance shape in the experiment iguchi . The conductance structure in Fig. 2 of Ref. iguchi,  is very similar to that in Fig. 3 (d). The present theory explains, at least, two typical conductance shapes observed in the experiments.

As shown in Fig. 3, the magnitude of the impurity potential and the position of impurities are key factors for the degree of splitting. In Fig.5, the gray area indicates sets of ($`V_iN_0,x_0k_F`$) which satisfy Eq. (83) within $`E_{\text{dip}}<E<E_{\text{ZEP}}`$. The open circles denote sets of ($`V_iN_0,x_0k_F`$), where we find the split of the ZBCP in Eq. (81).

All the circles are inside of the gray area. Although the circles and the gray region are not perfectly coincide with each other, they show qualitatively the same tendency. The parameters used in Fig. 3 are indicated by filled squares. Since $`s`$ at $`E=0`$ is proportional to the local density of states, we also show the $`N_s(E=0,x)`$ in Eq. (61). To satisfy Eq. (83), impurities around $`x_0\xi _0`$ should have sufficiently small scattering potentials because the local density of states has large values there. Thus the gray area appears for small impurity potentials near the interface, (i.e., $`x_0\xi _0`$). The gray region spreads to larger $`V_iN_0`$ as the increase of $`x_0>\xi _0`$ because the local density of states becomes smaller values. The results imply that the strong impurities are not necessary for the split of the ZBCP. It is evident that this phase diagram is valid in the limit of high impurity density and the diagram would be changed depending on the transmission probability of junctions.

In Fig. 3, all impurities are aligned at $`x_j=x_0`$. In real junctions, however, impurities may be distributed randomly near the interface as shown in Fig. 1. The conductance in such realistic junctions are shown in Fig. 6, where impurities are distributed randomly in the range of $`1<x_jk_F<L_sk_F`$, $`\rho _i=N_i\frac{\lambda _F^2}{WL_s}`$ is the dimensionless area density of impurities and $`z_0=10`$. The conductance is calculated from an expression

$`G_{NS}=`$ $`{\displaystyle \frac{2e^2}{h}}N_c{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\left({\displaystyle \frac{f_{\text{FD}}(EeV)}{(eV)}}\right)`$
$`\times \left[g^{(0)}\rho _i{\displaystyle \frac{k_FL_s}{\pi }}\mathrm{\Gamma }\right],`$ (84)
$`\mathrm{\Gamma }=`$ $`{\displaystyle \frac{1}{N_i}}{\displaystyle \underset{j=1}{\overset{N_i}{}}}\mathrm{\Gamma }_j,`$ (85)
where $`\mathrm{}`$ represents the ensemble average. Since the conductance in Eq. (75) is characterized by the number of impurities, a factor $`k_FL_s/\pi `$ appears in Eq. (84). We choose $`L_sk_F=20`$ in Fig. 6 because we focus on the impurities near the interface and $`\xi _0k_F6.3`$. We consider low density strong impurities in (a), where $`V_iN_0=0.1`$ and $`\rho _i=0.2`$. There is no peak splitting in Fig. 6 (a) because most impurities are outside of the gray region in Fig. 5. On the other hand, the results in Fig. 6 (b) show the split of the ZBCP at $`T=0`$, where we consider high density weak impurities with for $`V_iN_0=0.02`$ and $`\rho _i=0.6`$. This is because most impurities are inside of the gray region in Fig. 5. The splitting peaks merge into a single peak under finite temperatures such as $`T=0.1E_{\text{ZEP}}`$. In Fig. 6 (c), we show the temperature dependence of the zero-bias conductance in Fig. 6 (b). The results show the reentrant behavior of the zero-bias conductance, which was found in the experiment. covington  In Fig. 6 (d), the peak position ($`\delta eV`$) in (b) is plotted as a function temperatures. Since $`\delta eV`$ is about 0.15 $`E_{\text{ZEP}}`$ at $`T=0`$, peak splitting is washed out at high temperatures such as $`T=0.11E_{\text{ZEP}}`$.

High density impurities with weak random potential are responsible for the split of the ZBCP in low transparent junctions.

Several experiments covington ; greene  show a sensitivity of the conductance peaks to external magnetic fields. Here we discuss the conductance in the presence of magnetic fields. The effects of magnetic fields are taken into account phenomenologically by using the Aharonov-Bohm like phase shift ab1 ; zes  of a quasiparticle. Since the impurity scattering in magnetic fields itself is a difficult problem to solve analytically, we neglect the interplay between magnetic fields and impurity scatterings. Within the phenomenological theory, zes  effects of magnetic fields is considered by replacing $`E`$ in Eq. (76) by $`E+|\mathrm{\Delta }_0\mathrm{cos}\gamma \mathrm{sin}\gamma |\varphi _B`$ as

$`g_B^{(0)}=`$ $`{\displaystyle _{\pi /2}^{\pi /2}}𝑑\gamma {\displaystyle \frac{\mathrm{\Delta }_0^2\mathrm{cos}^7\gamma \mathrm{sin}^2\gamma }{E_B^2z_0^4+\mathrm{\Delta }_0^2\mathrm{cos}^6\gamma \mathrm{sin}^2\gamma }},`$ (86)
$`E_B=`$ $`E+2\mathrm{\Delta }_0|\mathrm{cos}\gamma \mathrm{sin}\gamma |\varphi _B,`$ (87)
$`\varphi _B=`$ $`2\pi {\displaystyle \frac{B\xi _0^2}{\varphi _0}}\mathrm{tan}\gamma =B_0\mathrm{tan}\gamma ,`$ (88)
where $`\varphi _0=2\pi \mathrm{}c/e`$ and $`B_0=1.0^3`$ corresponds to $`B=1`$ Tesra. A quasiparticle acquires the Aharonov-Bohm like phase shift $`\varphi _B`$ while moving near the NS interface. zes  In a previous paper, we found that ZBCP in clean junctions remains a single peak even in the strong magnetic fields zes  as shown in Fig. 7 (c), where $`z_0=10`$ and $`T=0.05E_{\text{ZEP}}`$. In Fig. 7 (a) and (b), we show the conductance in the presence of low density strong impurities and high density weak impurities, respectively, where $`V_i`$ and $`\rho _i`$ are same as those in Fig. 6 (a) and (b), respectively. A temperature is fixed at $`T=0.05E_{\text{ZEP}}`$. In contrast to clean junctions in Fig. 7 (c), the ZBCP in disordered junctions splits into two peaks under magnetic fields as shown in Fig. 7(a). The results within the phenomenological theory indicate that the sensitivity of the ZBCP to magnetic fields depends on the degree of impurity scatterings. In insets, peak positions $`(\delta eV)`$ are plotted with circles as a function of magnetic fields. For high density weak impurities in Fig. 7 (b), we also found that the degree of peak splitting increases with increasing magnetic fields. In the limit of the strong fields, $`\delta eV`$ tends to saturate as shown in the inset. These characteristic behavior are found in the experiment. covington

Although we have well explained the characteristic behavior of the experimental conductance peaks under magnetic fields, the applicability of the phenomenological theory in the presence of impurities is still unclear. This issue would be addressed more clearly in an exact numerical simulation.

## IV discussion

In experiments, the split of the ZBCP has been reported in overdoped high-$`T_c`$ superconductors. sharoni  The heavy carrier doping may bring a number of defects or imperfections in superconductors. It is evident that the impurity scattering is unavoidable even in underdoped high-$`T_c`$ superconductors. The split of the ZBCP would be found in underdoped superconductors if they have bad sample quality. In the same way, the split of the ZBCP would not be found even in overdoped superconductors if their sample quality are good enough. The sample quality is a key factor for the spilt of the ZBCP. This argument is consistent with an experiment, aprili  where the potential disorder is artificially introduced to the NS junctions by the ion-irradiation and washes out the ZBCP in the limit of strong disorder.

In order to compare the theoretical results in this paper with experiments, we may consider the impurity scattering in normal metals. The total resistance ($`R`$) in the dirty normal-metal / (110) $`d`$-wave junction can be described simply by $`R=R_D+R_{NS}`$circuit ; h-kashiwaya  where $`R_D=1/G_D`$ is the resistance of the dirty normal metal and $`G_{NS}=1/R_{NS}`$ is the conductance discussed in this paper. The equation indicates the absence of the proximity effect in the dirty normal metal. asano7 ; asano6  The height of the ZBCP can be reduced to reasonable values in the presence of the impurity scattering in normal metals because the total conductance is given by $`1/(R_D+R_{NS})`$. The degree of splitting depends on parameters such as the potential of impurities ($`V_j`$), the position of impurities ($`x_j`$) and the transparency of the junction ($`T_B1/z_0^2`$). In particular, $`T_B`$ is the most important parameter. As shown in Figs. 3 and 6, the degree of splitting is roughly given by $`E_{\text{ZEP}}=\mathrm{\Delta }_0/z_0^2\mathrm{\Delta }_0T_B`$. Thus it is possible to choose $`T_B`$ to fit the degree of splitting with that found in experiments. The amplitude of the pair potential is about 30-40 meV in typical high-$`T_c`$ materials. The degree of splitting is then estimated to be 1.2-1.6 meV for $`z_0=5`$, which is almost consistent with that found in experiments, for instance, 2 meV. covington  It is also necessary to consider electronic structures of high-$`T_c`$ materials for the quantitative agreement of the splitting width in theories with that in experiments.

In quasiclassical Green function theories, the conductance is proportional to the density of states at the surface of superconductors. We find that $`E`$ dependence of the density of states at the interface in Eq. (63) is apparently different from that of the conductance even in the clean junctions shown in Eq. (76). We also find that $`N_s(E,x=0)N_0N_s(E,x=\xi _0)`$ as shown in Fig. 2 (b). The density of states near the interface averaged over $`\xi _0`$ in this paper may corresponds to the surface density of states in the quasiclassical approximation, where the rapidly oscillating part of the wave function are neglected and the smallest length scales is given by $`\xi _0`$. It is impossible to directly compare the present theory with the quasiclassical Green function theories because the position of impurities is a key parameter for the split of the ZBCP in our theory.

The Abrikosov-Gor’kov (AG) theory agd  is a useful approach to discuss the impurity scattering in superconductors. The applicability of the AG theory is limited to the superconductivity in diffusive metals, where the mean free path is much smaller than the size of the disordered region in superconductors. Here we briefly discuss a relation between the AG theory and ours. Since the AG theory assumes the diffusive transport regime, we define the mean free path of a quasiparticle in the Born approximation,

$`\mathrm{}=`$ $`v_F\tau ,`$ (89)
$`{\displaystyle \frac{\mathrm{}}{\tau }}=`$ $`2\pi \overline{\rho }_iN_0V_i^2,`$ (90)
$`\overline{\rho }_i=`$ $`N_i/WL_s,`$ (91)
where we consider $`W\times L_s`$ disordered region on two-dimensional superconductor as in Fig. 6. The impurity density $`\overline{\rho }_i`$ can be replaced by the dimensionless impurity density $`\rho _i=(2\pi /k_F)^2\overline{\rho }_i`$. The ratio of the mean free path and the coherence length is given by

$$\frac{\mathrm{}}{\xi _0}=\frac{\mathrm{\Delta }_0}{\mu _F}\frac{\pi }{\rho _i(N_0V_i)^2}$$
(92)
In this paper, $`\mathrm{}/\xi _0`$ is larger than unity even in the limit of $`\rho _i=1`$ because we assume $`\mathrm{\Delta }_0/\mu _F=0.1`$ and $`N_0V_i<0.1`$. On the other hand, the dirty limit is defined by a relation $`\mathrm{}/\xi _01`$. When we choose $`L_s`$ being a few coherence length as shown in Fig. 6, where $`L_s3\xi _0`$, we find that $`\mathrm{}/L_s`$ is still larger than unity. The disordered region is not in the diffusive regime but in the quasi-ballistic regime because the diffusive regime is characterized by a relation $`\mathrm{}/L_s1`$. Thus it is basically impossible to apply the Abrikosov-Gor’kov theory to the model in this paper. The scattering theory asano5  used in this paper is a suitable analytic method to discuss the conductance of such NS junctions. Although the impurity scattering near the interface is very weak in normal states, it drastically affects the conductance below the critical temperature. Impurities located at the resonant states seriously suppress the degree of the resonance even if their potentials are weak. In the AG theory, the real part of the self-energy is usually neglected. In our theory, this approximation corresponds to an equation $`s_2=0`$. However $`s_2`$ plays an important role in the peak splitting. When we omit $`s_2`$, $`I_2`$ and $`I_3`$ must be also zero, which leads to positive $`\mathrm{\Gamma }_j`$ and no splitting irrespective of $`V_iN_0`$ and $`x_jk_F`$.

We show that the impurity scattering causes the splitting of the ZBCP. The conclusion, however, does not deny a possibility of the BTRSS. In a previous paper, kitaura  we assume the $`s+id`$ symmetry near the NS interface and numerically study the tunneling conductance. The results show that the split of the ZBCP is insensitive to the potential disorder. Thus peak splitting would be always expected in low transparent junctions Tanuma2001  if the BTRSS appears at the NS interface. At present, we have only limited information on the BTRSS within the mean-field theories. To understand the nature of the BTRSS further beyond the mean-field theory, we have to analyze electronic structure of high-$`T_c`$ superconductors based on microscopic models and make clear effects of the surface, the electron correlation and the random potentials on the superconducting state. This is an important future problem.

Since the formation of the ZES is a universal phenomenon in superconductors with unconventional pairing symmetries, the ZES is also expected at a surface of spin-triplet superconductors. Buch  It is interesting to study effects of impurities on transport properties in spin-triplet superconductor junctions. Yama1 ; Yama2 ; Yama3 ; asano2003-3 ; Y1 ; H1 ; H2 ; Kuroki1 ; Kuroki2 ; Kusakabe ; Honer ; St ; Sen ; asano6 ; asano8 ; asano9

## V conclusion

We have discussed effects of impurity scatterings on the conductance in normal-metal/$`d`$ wave superconductor junctions. The conductance is calculated from the Andreev and the normal reflection coefficients which are estimated by using the single-site approximation. We consider impurities near the junction interface on the superconductor side. The strength of the impurity scattering strongly depends on the transparency of the junction, the position of impurities and the energy of a quasiparticle because the ZES are formed at the NS interface. We conclude that the impurity scattering causes the split of the zero-bias conductance peak. The results are consistent with previous numerical simulations. We have also shown that characteristic behaviors of the conductance spectra at finite temperatures and under external magnetic fields qualitatively agree with those reported in experiments.

## Appendix A Transmission and reflection coefficients

In the clean NS junctions, the transmission and the reflection coefficients can be calculated from the appropriate boundary condition of the wave function. The calculated results are shown below.

$`r_{NN}^{he}(l)=`$ $`{\displaystyle \frac{\overline{k}_l^2}{\mathrm{\Xi }_l}}{\displaystyle \frac{\mathrm{\Delta }_l}{2}}e^{i\phi },`$ (93)
$`r_{NN}^{ee}(l)=`$ $`{\displaystyle \frac{iz_0(\overline{k}_liz_0)}{\mathrm{\Xi }_l}}E,`$ (94)
$`t_{SN}^{ee}(l)=`$ $`{\displaystyle \frac{\overline{k}_l(\overline{k}_liz_0)}{\mathrm{\Xi }_l}}Eu_le^{i\phi /2},`$ (95)
$`t_{SN}^{he}(l)=`$ $`{\displaystyle \frac{iz_0\overline{k}_l}{\mathrm{\Xi }_l}}Ev_le^{i\phi /2},`$ (96)
$`r_{NN}^{eh}(l)=`$ $`{\displaystyle \frac{\overline{k}_l^2}{\mathrm{\Xi }_l}}{\displaystyle \frac{\mathrm{\Delta }_l}{2}}e^{i\phi },`$ (97)
$`r_{NN}^{hh}(l)=`$ $`{\displaystyle \frac{iz_0(\overline{k}_l+iz_0)}{\mathrm{\Xi }_l}}E,`$ (98)
$`t_{SN}^{hh}(l)=`$ $`{\displaystyle \frac{\overline{k}_l(\overline{k}_l+iz_0)}{\mathrm{\Xi }_l}}Eu_le^{i\phi /2},`$ (99)
$`t_{SN}^{eh}(l)=`$ $`{\displaystyle \frac{iz_0\overline{k}_l}{\mathrm{\Xi }_l}}Ev_le^{i\phi /2},`$ (100)
$`r_{SS}^{he}(l)=`$ $`{\displaystyle \frac{\overline{k}_l^2+2z_0^2}{\mathrm{\Xi }_l}}{\displaystyle \frac{\mathrm{\Delta }_l}{2}},`$ (101)
$`r_{SS}^{ee}(l)=`$ $`{\displaystyle \frac{iz_0(\overline{k}_liz_0)}{\mathrm{\Xi }_l}}\mathrm{\Omega }_l,`$ (102)
$`t_{NS}^{ee}(l)=`$ $`{\displaystyle \frac{\overline{k}_l(\overline{k}_liz_0)}{\mathrm{\Xi }_l}}\mathrm{\Omega }_lu_le^{i\phi /2},`$ (103)
$`t_{NS}^{he}(l)=`$ $`{\displaystyle \frac{iz_0\overline{k}_l}{\mathrm{\Xi }_l}}\mathrm{\Omega }_lv_le^{i\phi /2},`$ (104)
$`r_{SS}^{eh}(l)=`$ $`{\displaystyle \frac{(\overline{k}_l^2+2z_0^2)}{\mathrm{\Xi }_l}}{\displaystyle \frac{\mathrm{\Delta }_l}{2}},`$ (106)
$`r_{SS}^{hh}(l)=`$ $`{\displaystyle \frac{iz_0(\overline{k}_l+iz_0)}{\mathrm{\Xi }_l}}\mathrm{\Omega }_l,`$ (107)
$`t_{NS}^{hh}(l)=`$ $`{\displaystyle \frac{\overline{k}_l(\overline{k}_l+iz_0)}{\mathrm{\Xi }_l}}\mathrm{\Omega }_lu_le^{i\phi /2},`$ (108)
$`t_{NS}^{eh}(l)=`$ $`{\displaystyle \frac{iz_0\overline{k}_l}{\mathrm{\Xi }_l}}\mathrm{\Omega }_lv_le^{i\phi /2},`$ (109)
$`\mathrm{\Xi }_l=`$ $`Ez_0^2+\overline{k}_l^2\left({\displaystyle \frac{E+\mathrm{\Omega }_l}{2}}\right).`$ (110)
For instance, $`t_{NS}^{he}(l)`$ is the transmission coefficients from the electron branch in a superconductor to the hole branch in a normal metal. In above coefficients, we use a relation $`q_l^\pm =k_l^\pm k_l`$ for simplicity. The conclusions in this paper remain unchanged in this approximation.

## Appendix B Green functions

The real space retarded Green function in clean junctions can be calculated by using the transmission and the reflection coefficients in Appendix A. For $`x<x^{}<0`$, the Green function from a normal metal to a normal metal is

$`\widehat{G}_0^{NN}(𝒓,𝒓^{})=i{\displaystyle \frac{\pi N_0}{W}}{\displaystyle \underset{k_y^{(l)}}{}}e^{ik_y^{(l)}(yy^{})}`$
$`\times `$ $`[{\displaystyle \frac{1}{q^+}}\{\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)e^{iq_l^+|xx^{}|}+\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)e^{iq_l^{}x}e^{iq_l^+x^{}}r_{NN}^{he}(l)+\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)e^{iq_l^+(x+x^{})}r_{NN}^{ee}(l)\}`$ (117)
$`+{\displaystyle \frac{1}{q^{}}}\{\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)e^{iq_l^{}|xx^{}|}+\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)e^{iq_l^+x}e^{iq_l^{}x^{}}r_{NN}^{eh}(l)+\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)e^{iq_l^{}(x+x^{})}r_{NN}^{hh}(l)\}],`$ (124)
$`N_0={\displaystyle \frac{m}{\pi \mathrm{}^2}}.`$ (125)
For $`x>x^{}>0`$, the Green function from a superconductor to a superconductor is

$`\widehat{G}_0^{SS}(𝒓,𝒓^{})=i{\displaystyle \frac{\pi N_0}{W}}{\displaystyle \underset{k_y^{(l)}}{}}e^{ik_y^{(l)}(yy^{})}{\displaystyle \frac{E}{\mathrm{\Omega }_l}}`$
$`\times `$ $`\widehat{\mathrm{\Phi }}[{\displaystyle \frac{1}{k^+}}\{\left(\begin{array}{cc}u_l^2& u_lv_l\\ u_lv_l& v_l^2\end{array}\right)e^{ik_l^+|xx^{}|}+\left(\begin{array}{cc}u_lv_l& v_l^2\\ u_l^2& u_lv_l\end{array}\right)e^{ik_l^{}x+ik_l^+x^{}}r_{SS}^{he}(l)+\left(\begin{array}{cc}u_l^2& u_lv_l\\ u_lv_l& v_l^2\end{array}\right)e^{ik_l^+(x+x^{})}r_{SS}^{ee}(l)\}`$ (132)
$`+{\displaystyle \frac{1}{k^{}}}\{\left(\begin{array}{cc}v_l^2& u_lv_l\\ u_lv_l& u_l^2\end{array}\right)e^{ik_l^{}|xx^{}|}+\left(\begin{array}{cc}u_lv_l& u_l^2\\ v_l^2& u_lv_l\end{array}\right)e^{ik_l^+xik_l^{}x^{}}r_{SS}^{eh}(l)+\left(\begin{array}{cc}v_l^2& u_lv_l\\ u_lv_l& u_l^2\end{array}\right)e^{ik_l^{}(x+x^{})}r_{SS}^{hh}(l)\}]\widehat{\mathrm{\Phi }}^{},`$ (139)
$`\widehat{\mathrm{\Phi }}=\left(\begin{array}{cc}e^{i\frac{\phi }{2}}& 0\\ 0& e^{i\frac{\phi }{2}}\end{array}\right),`$ (142)
where $`\phi `$ is the phase of a superconductor. For $`x>0>x^{}`$, the Green function from a normal metal to a superconductor is

$`\widehat{G}_0^{SN}(𝒓,𝒓^{})=`$ $`i{\displaystyle \frac{\pi N_0}{W}}{\displaystyle \underset{k_y^{(l)}}{}}e^{ik_y^{(l)}(yy^{})}`$
$`\times `$ $`\widehat{\mathrm{\Phi }}[{\displaystyle \frac{1}{q^+}}\{\left(\begin{array}{cc}u_l& 0\\ v_l& 0\end{array}\right)e^{ik_l^+xiq_l^+x^{}}t_{SN}^{ee}(l)+\left(\begin{array}{cc}v_l& 0\\ u_l& 0\end{array}\right)e^{ik_l^{}xiq_l^+x^{}}t_{SN}^{he}(l)\}`$ (147)
$`+{\displaystyle \frac{1}{q^{}}}\{\left(\begin{array}{cc}0& v_l\\ 0& u_l\end{array}\right)e^{ik_l^{}x+iq^{}x^{}}t_{SN}^{hh}(l)+\left(\begin{array}{cc}0& u_l\\ 0& v_l\end{array}\right)e^{ik_l^+x+iq_l^{}x^{}}t_{SN}^{eh}(l)\}].`$ (152)
For $`x<0<x^{}`$, the Green function from a superconductor to a normal metal is

$`\widehat{G}_0^{NS}(𝒓,𝒓^{})=`$ $`i{\displaystyle \frac{\pi N_0}{W}}{\displaystyle \underset{k_y^{(l)}}{}}e^{ik_y^{(l)}(yy^{})}{\displaystyle \frac{E}{\mathrm{\Omega }_l}}`$
$`\times `$ $`[{\displaystyle \frac{1}{k^+}}\{\left(\begin{array}{cc}u_l& v_l\\ 0& 0\end{array}\right)e^{iq_l^+x+ik_l^+x^{}}t_{NS}^{ee}(l)+\left(\begin{array}{cc}0& 0\\ u_l& v_l\end{array}\right)e^{iq_lx+ik_l^+x^{}}t_{NS}^{he}(l)\}`$ (157)
$`+{\displaystyle \frac{1}{k^{}}}\{\left(\begin{array}{cc}0& 0\\ v_l& u_l\end{array}\right)e^{iq_l^{}xik_l^{}x^{}}t_{NS}^{hh}(l)+\left(\begin{array}{cc}v_l& u_l\\ 0& 0\end{array}\right)e^{iq_l^+xik_l^{}x^{}}t_{NS}^{eh}(l)\}]\widehat{\mathrm{\Phi }}^{}.`$ (162)


format: text

text:
# Contents

## 1 Introduction

The physics of massive and mixed neutrinos is one of the hot topics in today’s research in high energy physics, following the discovery of oscillations of atmospheric and solar neutrinos (see the reviews in Refs.  and references therein, as well as the references in ).

Neutrino oscillations have been proposed in the late 50’s by Pontecorvo . The oscillations are generated by the interference of different massive neutrinos, which are produced and detected coherently because of their very small mass difference.

In 1962 Maki, Nakagawa and Sakata  considered for the first time a model with mixing of different neutrino flavors. In 1967 Pontecorvo proposed the possibility of solar neutrino oscillations , before the discovery in 1968 of the solar neutrino problem in the Homestake experiment . In 1969 Pontecorvo and Gribov considered $`\nu _e\nu _\mu `$ oscillations as a possible explanation of the solar neutrino problem .

The theory of neutrino oscillations in the plane-wave approximation was developed in the middle 70’s by Eliezer and Swift , Fritzsch and Minkowski , Bilenky and Pontecorvo , and beautifully reviewed by Bilenky and Pontecorvo in Ref. .

In 1976 Nussinov  for the first time considered the wave packet nature of propagating neutrinos and inferred the existence of a coherence length, beyond which the interference of different massive neutrinos is not observable. This is due to the different group velocities of different massive neutrinos, that causes a separation of their wave packets. In 1996 Kiers, Nussinov and Weiss  first pointed out the importance of the detection process for the coherence of neutrino oscillations and discussed some implications for the wave packet approach.

In 1981 Kayser  presented the first detailed discussion of the quantum mechanical problems of neutrino oscillations, pointing out the necessity of a wave packet treatment.

Wave packet models of neutrino oscillations have been later developed in the framework of quantum mechanics  and in the framework of quantum field theory  (see also Ref.  and the reviews in Refs. ).

In spite of the well-known fact that in quantum theory localized particles are described by wave packets and in spite of the success of the wave packet treatment of neutrino oscillations, some authors have presented arguments against such approach. In Section 2 we briefly review well-known evidences of the wave packet nature of light, which, by analogy, imply a wave packet nature of neutrinos. In Sections 3 and 4 we present objections to arguments against a wave packet description of neutrinos. In Section 4 we present also a wave packet derivation of neutrino oscillations in the density matrix formalism and we discuss some implications for the case of stationary beams.

Another debated problem in the theory of neutrino oscillations is the determination of the energies of massive neutrinos. In Section 5 we present objections to recently proposed arguments in favor of the equality of the energies of different massive neutrinos.

## 2 No source of waves vibrates indefinitely

The title of this section is taken from the beginning sentence of Section 11.11 of Ref. , which concludes with the following paragraph:

> In light sources, the radiating atoms emit wave trains of finite length. Usually, because of collisions or damping arising from other causes, these packets are very short. According to the theorem mentioned above111 “The largest the number $`N`$ of waves in the group, the smaller the spread $`\mathrm{\Delta }\lambda `$, and in fact theory shows that $`\mathrm{\Delta }\lambda /\lambda _0`$ is approximately equal to $`1/N`$.” (Section 11.11 of Ref. .) , the consequence is that the spectrum lines will not be very narrow but will have an appreciable width $`\mathrm{\Delta }\lambda `$. A measurement of this width will yield the effective “lifetime” of the electromagnetic oscillators in the atoms and the average length of the wave packets. A low-pressure discharge through the vapor of mercury containing the single isotope $`{}_{}{}^{198}\mathrm{Hg}`$ yields very sharp spectral lines, of width about $`0.005\text{ }\mathrm{\AA }`$. Taking the wavelength of one of the brightest lines, $`5461\text{ }\mathrm{\AA }`$, we may estimate that there are roughly $`10^6`$ waves in a packet and that the packets themselves are some 50 cm long.

The broadening of optical lines due to the finite lifetime of atomic transitions is known as *natural linewidth*.

It is interesting to note that the broadening of optical lines was well known to experimental physicists in the nineteen century and explained by classical models before the advent of quantum theory (see Ref.  and references therein). In 1895 Michelson  listed among the hypotheses formulated before that time to account for line broadening

> 3. The exponential diminution in amplitude of the vibrations due to communication of energy to the surrounding medium or to other causes.

As explained in ,

> We consider an emitting atom which we shall proceed to remove to infinity and reduce the “temperature” to the point where, classically at least, no translational motion exists. Now from the classical picture of a vibrating electron or the simple picture of a pair of energy levels between which our radiation transition takes place, we should expect these conditions to yield a spectral line of a single frequency. We, of course, do not obtain this result, but, rather, we obtain the familiar natural line shape which is attributable to Michelson’s Cause 3.

A classical derivation of the natural linewidth can be found in Section 17.7 of Ref. . Let us emphasize that the natural linewidth of atomic lines has been *observed experimentally* (see Section 21.4 of Ref. )!

Another important cause of line broadening known in 1895 was

> 4. The change in wavelength due to the Doppler effect of the component of the velocity of the vibrating atom in the line of sight.

The Doppler line broadening due to the thermal motion of atoms in a medium was calculated by Rayleigh in 1889 . This important effect, that must be aways taken into account in calculating the spectral shape of monochromatic beams, does not concern us here, because it does not generate a coherent broadening. It is simply due to the different motion of different atoms, whose radiation is incoherent.

To these (and other) causes of line broadening Michelson added in 1895

> 5. The limitation of the number of regular vibrations by more or less abrupt changes of phase amplitude or plane of vibration caused by collisions.

This important coherent effect has been called with several names, among which *collision broadening*, *pressure broadening* and *interruption broadening*, it has been studied in depth by many authors (see, for example, Refs.) and it has been *observed experimentally* (see Section 21.5 of Ref. )!

In quantum theory, the fact that “no source of waves vibrates indefinitely” implies that all particles are produced as wave packets, whose size is determined by the finite lifetime of the parent particle, or by its finite mean free path if the production process occurs in a medium. Since also no wave detector vibrates indefinitely, it is clear that all particles are also detected as wave packets.

Usually, at least in high energy physics, the wave packet character of particles is not important and they can be well approximated by plane waves. The phenomenon of neutrino oscillations is an exception, as shown in 1981 by Kayser , because the localization of the source and detector requires a wave packet treatment.

The wave packet approach to neutrino oscillations implemented in Refs.  in the framework of quantum mechanics and in Refs.  in the framework of quantum field theory allows to derive the neutrino oscillation probability from first principles in a consistent framework.

In the light of these considerations, it seems rather surprising that some physicists do not agree with a wave packet description of neutrinos and put forward arguments in favor of a unique value of neutrino energy. A tentative explanation of this approach stems from the common education and practice in physics of working with energy eigenstates, forgetting that they are only approximations of the states of physical systems in the real world222 I would like to thank I. Pesando for an illuminating discussion about this point. . In the following Sections we present a critical discussion of recently advocated approaches to neutrino oscillations in which neutrinos are not described by wave packets with some spread in momentum and energy.

## 3 Bound states

Bound states of atomic electrons are usually calculated with different degrees of approximation taking into account the static potential generated by the nucleus and other electrons and neglecting the coupling of the electron with the electromagnetic radiation field. In this approximation, bound states are stationary and have definite energy. In spite of the obvious fact that in such approximation electrons cannot jump from one bound state to another, the energy difference between two stationary bound states is sometimes associated without any further explanation with the energy of the photon emitted or absorbed in the transition of an electron from one bound state to the other.

However, in order to describe the transitions between bound states it is necessary to take into account the coupling of the electron with the time-dependent electromagnetic radiation field. Since the bound states are not eigenstates of the full Hamiltonian, in reality they are not stationary and they do not have definite energy. The lack of a definite energy is a necessary requirement for the existence of transitions: since the Schrödinger equation implies that the time evolution of energy eigenstates evolve by a phase, orthogonality between different energy eigenstates is conserved in time, excluding transitions.

The fact that unstable systems do not have a definite energy is well-known to high energy physicists from their experience with resonances (hadronic excited states) and highly unstable particles ($`Z`$ and $`W`$ bosons), which do not have a definite mass (energy in the rest frame). Their mass width is inversely proportional to the decay rate.

From these considerations it follows that the quantum field theoretical model of neutrino oscillations in the process

$$np+e^{}+\mathrm{bar}\nu _e\stackrel{\mathrm{bar}\nu _e\mathrm{bar}\nu _e}{}\mathrm{bar}\nu _e+e^{}\mathrm{bar}\nu _e+e^{}$$
(3.1)
presented in Ref. , in which a neutrino is assumed to be produced by a nuclear $`\beta `$-decay between two stationary bound states with definite energy and detected through scattering with an electron in a stationary atomic state with definite energy, is unrealistic, not to mention the description of the final-state electrons by plane waves, which is in contradiction with the fact that the electrons interact with the surrounding medium and with the necessity to observe the final lepton in the detection process, as already noted in . Of course, once the matters of principle are clear, one can consider the assumptions in Ref.  as approximations acceptable in some cases. One must also notice that Ref.  is very interesting for the useful theorem proved in the appendix, which has been used in the quantum field theoretical calculation of neutrino oscillations in the wave packet approach in Refs. .

Developing the technique proposed in Ref. , the authors of Refs.  considered neutrino oscillations in the process

$$\mu ^+e^++\nu _e+\mathrm{bar}\nu _\mu \stackrel{\mathrm{bar}\nu _\mu \mathrm{bar}\nu _e}{}\mathrm{bar}\nu _e+pn+e^+$$
(3.2)
with the muon correctly described by a wave packet, but the proton was still described by an unrealistic stationary bound state with definite energy. Furthermore, all final particles were described by plane waves, in contradiction with the fact that the positrons and neutron interact with the surrounding medium and the necessity to observe the final positron in the detection process in order to measure oscillations. Localized interactions of a particle with a detector obviously imply that the particle cannot be described by an unlocalized plane wave and a wave packet description is necessary.

## 4 Stationary beams

The author of the very interesting Ref.  wrote:

> …many wavepacket discussions for the coherence properties of particle beams are unnecessary since they deal with stationary sources; and when the problem is stationary, essentially all information is in the energy spectrum.

Indeed, as proved in Ref. ,

> …under very general assumptions it is not possible to distinguish experimentally neutrinos produced in some region of space as wave packets from those produced in the same region of space as plane waves with the same energy distribution.

It is not clear if the purpose of Ref.  was to show that neutrinos are not described by wave packets, but it is certain that it has been interpreted in this way by some physicists. However, Ref.  certainly does not provide any help for the derivation of neutrino oscillations in the framework of quantum field theory. Even in the framework of quantum mechanics Ref.  does not go beyond the well-known discussion of neutrino oscillations assuming equal energies for the different massive neutrino components , which has been shown to be unrealistic in Refs.  and is further criticized in the Section 5.

An interpretation of Ref.  as a proof that neutrinos are not described by wave packets seems to stem from a confusion between microscopic and macroscopic stationarity, as already remarked in . Since most neutrino sources emit neutrinos at a constant rate over macroscopic intervals of time, they can be considered macroscopically stationary. However, the microscopic processes of neutrino production (nuclear or particle decay, etc.) are certainly not stationary and the single neutrino does not know if it will be part of a stationary beam or not. Because of the localization in space and time of the production and detection processes, a description of neutrino oscillations in terms of wave packets is inescapable.

Of course, after the matters of principle are settled, neutrino oscillations are derived from first principles in a consistent theoretical framework, and an estimation of the coherence length has shown that it is much longer than the source-detector distance in a neutrino oscillation experiment, the analysis of the experimental data can be performed using the standard oscillation probability obtained in the plane wave approximation. This turns out to be the case in all present-day neutrino oscillation experiments.

In the following part of this section we present a wave packet derivation of neutrino oscillations in the density matrix formalism, which is suitable for the description of a stationary beam, as emphasized in Ref. . Our calculation follows a method similar to the one presented in Ref. .

In a quantum-mechanical wave packet treatment of neutrino oscillations a neutrino produced at the origin of the space-time coordinates by a weak interaction process with definite flavor $`\alpha `$ ($`\alpha =e,\mu ,\tau `$) and propagating along the $`x`$ axis is described by the state

$$|\nu _\alpha (x,t)=\underset{k}{}U_{\alpha k}^{}\psi _k(x,t)|\nu _k,$$
(4.1)
where $`U`$ is the mixing matrix, $`|\nu _k`$ is the state of a neutrino with mass $`m_k`$ and $`\psi _k(x,t)`$ is its wave function. Approximating the momentum distribution of the massive neutrino $`\nu _k`$ with a gaussian,

$$\psi _k(p)=\left(2\pi \sigma _{p}^{\mathrm{P}}{}_{}{}^{2}\right)^{1/4}\mathrm{exp}\left[\frac{\left(pp_k\right)^2}{4\sigma _{p}^{\mathrm{P}}{}_{}{}^{2}}\right],$$
(4.2)
where $`p_k`$ is the average momentum and $`\sigma _p^\mathrm{P}`$ is the momentum uncertainty determined by the production process, the wave function is

$$\psi _k(x,t)=\frac{\mathrm{d}p}{\sqrt{2\pi }}\psi _k(p)e^{ipxiE_k(p)t},$$
(4.3)
with the energy $`E_k(p)`$ given by333 We do not make unnecessary assumptions about the values of the energies and momenta of different massive neutrinos. The claim that different massive neutrinos must have the same energy has been confuted in Refs.  and in further discussed in Section 5.

$$E_k(p)=\sqrt{p^2+m_k^2}.$$
(4.4)
If the gaussian momentum distribution (4.2) is sharply peaked around the mean momentum $`p_k`$ (i.e. if $`\sigma _p^\mathrm{P}E_k^2(p_k)/m_k`$), the energy $`E_k(p)`$ can be approximated by

$$E_k(p)E_k+v_k\left(pp_k\right),$$
(4.5)
where

$$E_kE_k(p_k)=\sqrt{p_k^2+m_k^2}$$
(4.6)
is the average energy and

$$v_k\frac{E_k(p)}{p}|_{p=p_k}=\frac{p_k}{E_k}.$$
(4.7)
Under this approximation, the integration over $`p`$ in Eq. (4.3) is gaussian and leads to the wave packet in coordinate space

$$\psi _k(x,t)=\left(2\pi \sigma _{x}^{\mathrm{P}}{}_{}{}^{2}\right)^{1/4}\mathrm{exp}\left[iE_kt+ip_kx\frac{\left(xv_kt\right)^2}{4\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}}\right],$$
(4.8)
where

$$\sigma _x^\mathrm{P}=\frac{1}{2\sigma _p^\mathrm{P}}$$
(4.9)
is the spatial width of the wave packet. From Eq. (4.8) one can see that $`v_k`$ is the group velocity of the wave packet of the massive neutrino $`\nu _k`$.

The pure state in Eq. (4.1) can be written in the form of a density matrix operator as

$$\widehat{\rho }_\alpha (x,t)=|\nu _\alpha (x,t)\nu _\alpha (x,t)|.$$
(4.10)
Using Eq. (4.8), we have

$`\widehat{\rho }_\alpha (x,t)={\displaystyle \frac{1}{\sqrt{2\pi \sigma _{x}^{\mathrm{P}}{}_{}{}^{2}}}}{\displaystyle \underset{k,j}{}}U_{\alpha k}^{}U_{\alpha j}\mathrm{exp}[i(E_kE_j)t+i(p_kp_j)x`$
$`{\displaystyle \frac{\left(xv_kt\right)^2}{4\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}}}{\displaystyle \frac{\left(xv_jt\right)^2}{4\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}}}]|\nu _k\nu _j|.`$ (4.11)
This space and time dependent density matrix operator describes neutrino oscillations in space and time. Although in laboratory experiments it is possible to measure neutrino oscillations in time through the measurement of both the production and detection processes (see Ref. )444 Let us also mention the possibility of time-dependent flavor oscillations of neutrinos in the primordial cosmological plasma (see Ref.  and references therein). , in all existing neutrino oscillation experiments only the source-detector distance is known. In this case the relevant density matrix operator $`\widehat{\rho }_\alpha (x)`$ is given by the time average of $`\widehat{\rho }_\alpha (x,t)`$. Since the integral over time is gaussian, one easily obtains555 The density matrix operator in Eq. (4.12) is normalized by $`_\alpha \widehat{\rho }_\alpha (x)=\widehat{I}`$, where $`\widehat{I}`$ is the identity operator.

$`\widehat{\rho }_\alpha (x)={\displaystyle \underset{k,j}{}}U_{\alpha k}^{}U_{\alpha j}\mathrm{exp}\{i[{\displaystyle \frac{v_k+v_j}{v_k^2+v_j^2}}(E_kE_j)(p_kp_j)]x`$
$`{\displaystyle \frac{\left(v_kv_j\right)^2x^2}{4\left(v_k^2+v_j^2\right)\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}}}{\displaystyle \frac{\left(E_kE_j\right)^2}{4\left(v_k^2+v_j^2\right)\sigma _{p}^{\mathrm{P}}{}_{}{}^{2}}}\}|\nu _k\nu _j|.`$ (4.12)
Since this density matrix operator is time independent, it is suitable for the description of a stationary beam in neutrino oscillation experiments, as noticed in Ref. .

In order to obtain the oscillation probability, it is convenient to simplify Eq. (4.12) for the realistic case of extremely relativistic neutrinos (see the discussion in Section 2 of Ref. ). In general, the average massive neutrino energies $`E_k`$ can be written as

$$E_kE+\xi _\mathrm{P}\frac{m_k^2}{2E},$$
(4.13)
where $`E`$ is the neutrino energy in the limit of zero mass and $`\xi _\mathrm{P}`$ is a dimensionless quantity that depends from the characteristics of the production process666 The quantity $`\xi _\mathrm{P}`$ cannot be calculated in a quantum mechanical framework, but its value can be estimated from energy-momentum conservation in the production process. The calculation of the value of $`\xi _\mathrm{P}`$ requires a quantum field theoretical treatment (see Ref. ). . From Eq. (4.4), the corresponding momentum in the relativistic approximation is

$$p_kE\left(1\xi _\mathrm{P}\right)\frac{m_k^2}{2E}.$$
(4.14)
The generality of these relations can be understood by noting that, since the neutrino mass $`m_k`$ enters quadratically in the energy-momentum dispersion relation (4.4), its the first order contribution to the energy and momentum must be proportional to $`m_k^2`$ and inversely proportional to $`E`$, which is the only available quantity with dimension of energy.

Using Eqs. (4.13), (4.14) and the relativistic approximation

$$v_k1\frac{m_k^2}{2E_k^2},$$
(4.15)
we obtain

$$\widehat{\rho }_\alpha (x)=\underset{k,j}{}U_{\alpha k}^{}U_{\alpha j}\mathrm{exp}\left[i\frac{\mathrm{\Delta }m_{kj}^2x}{2E}\left(\frac{\mathrm{\Delta }m_{kj}^2x}{4\sqrt{2}E^2\sigma _x^\mathrm{P}}\right)^2\left(\xi _\mathrm{P}\frac{\mathrm{\Delta }m_{kj}^2}{4\sqrt{2}E\sigma _p^\mathrm{P}}\right)^2\right]|\nu _k\nu _j|,$$
(4.16)
with $`\mathrm{\Delta }m_{kj}^2m_k^2m_j^2`$.

In analogy with the production process, we describe the process of detection of a neutrino with flavor $`\beta `$ at the coordinate $`x=L`$ with the operator

$`\widehat{𝒪}_\beta `$ $`(xL)={\displaystyle \underset{k,j}{}}U_{\beta k}^{}U_{\beta j}`$
$`\times \mathrm{exp}\left[i{\displaystyle \frac{\mathrm{\Delta }m_{kj}^2\left(xL\right)}{2E}}\left({\displaystyle \frac{\mathrm{\Delta }m_{kj}^2\left(xL\right)}{4\sqrt{2}E^2\sigma _x^\mathrm{D}}}\right)^2\left(\xi _\mathrm{D}{\displaystyle \frac{\mathrm{\Delta }m_{kj}^2}{4\sqrt{2}E\sigma _p^\mathrm{D}}}\right)^2\right]|\nu _k\nu _j|,`$ (4.17)
where $`\xi _\mathrm{D}`$ is a dimensionless quantity that depends from the characteristics of the detection process (see footnote 6), $`\sigma _p^\mathrm{D}`$ is the momentum uncertainty of the detection process and $`\sigma _x^\mathrm{D}=1/2\sigma _p^\mathrm{D}`$.

The probability of $`\nu _\alpha \nu _\beta `$ transitions is given by

$$P_{\nu _\alpha \nu _\beta }(L)=\mathrm{Tr}\left(\widehat{\rho }_\alpha (x)\widehat{𝒪}_\beta (xL)\right)=dx\underset{k}{}\nu _k|\widehat{\rho }_\alpha (x)\widehat{𝒪}_\beta (xL)|\nu _k,$$
(4.18)
which leads to

$$P_{\nu _\alpha \nu _\beta }(L)=\underset{k,j}{}U_{\alpha k}^{}U_{\alpha j}U_{\beta k}U_{\beta j}^{}\mathrm{exp}\left[2\pi i\frac{L}{L_{kj}^{\mathrm{osc}}}\left(\frac{L}{L_{kj}^{\mathrm{coh}}}\right)^22\pi ^2\xi ^2\left(\frac{\sigma _x}{L_{kj}^{\mathrm{osc}}}\right)^2\right],$$
(4.19)
with the oscillation and coherence lengths

$`L_{kj}^{\mathrm{osc}}={\displaystyle \frac{4\pi E}{\mathrm{\Delta }m_{kj}^2}},`$ (4.20)
$`L_{kj}^{\mathrm{coh}}={\displaystyle \frac{4\sqrt{2}E^2}{|\mathrm{\Delta }m_{kj}^2|}}\sigma _x,`$ (4.21)
and

$`\sigma _x^2=\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}+\sigma _{x}^{\mathrm{D}}{}_{}{}^{2},`$ (4.22)
$`\xi ^2\sigma _x^2=\xi _\mathrm{P}^2\sigma _{x}^{\mathrm{P}}{}_{}{}^{2}+\xi _\mathrm{D}^2\sigma _{x}^{\mathrm{D}}{}_{}{}^{2}.`$ (4.23)
This result was already obtained777 The only improvement in the expressions above with respect to the analogous ones in Ref.  is the introduction of $`\xi _\mathrm{D}`$, that takes into account the properties of the detection process.  in a quantum mechanical framework in Ref. , where the wave packet treatment with pure states presented several years before in Ref.  was extended in order to take into account the coherence properties of the detection process, whose importance was first recognized in Ref. . Expressions for the $`\nu _\alpha \nu _\beta `$ transition probability similar to Eq. (4.19) have been also obtained with wave packet treatments of neutrino oscillations in a quantum field theoretical framework in Refs. .

It is important to remark that the wave packet treatment of neutrino oscillations confirms the standard value in Eq. (4.20) for the oscillation length $`L_{kj}^{\mathrm{osc}}`$. The coherence length $`L_{kj}^{\mathrm{coh}}`$ in Eq. (4.21) is the distance beyond which the interference of the massive neutrinos $`\nu _k`$ and $`\nu _j`$ is suppressed because the separation of their wave packets when they arrive at the detector is so large that they cannot be absorbed coherently. As shown by Eqs. (4.21) and (4.22) the coherence lengths $`L_{kj}^{\mathrm{coh}}`$ are proportional to the total coherence size, that is dominated by the largest between the coherence sizes $`\sigma _x^\mathrm{P}`$ and $`\sigma _x^\mathrm{D}`$ of the production and detection process. The last term in the exponential of Eq. (4.19) implies that the interference of the massive neutrinos $`\nu _k`$ and $`\nu _j`$ is observable only if the localization of the production and detection processes is smaller than the oscillation length.

The localization term is important for the distinction of neutrino oscillation experiments from experiments on the measurement of neutrino masses. As first shown by Kayser in Ref. , neutrino oscillations are suppressed in experiments able to measure the value of a neutrino mass, because the measurement of a neutrino mass implies that only the corresponding massive neutrino is produced or detected.

Kayser’s  argument goes as follows. Since a neutrino mass is measured from energy-momentum conservation in a process in which a neutrino is produced or detected, from the energy-momentum dispersion relation (4.4), the uncertainty of the mass determination is

$$\delta m_{k}^{}{}_{}{}^{2}=\sqrt{\left(2E_k\delta E_k\right)^2+\left(2p_k\delta p_k\right)^2}2\sqrt{2}E\sigma _p,$$
(4.24)
where the approximation holds for realistic extremely relativistic neutrinos and $`\sigma _p`$ is the momentum uncertainty. If the mass of $`\nu _k`$ is measured with an accuracy better than $`\mathrm{\Delta }m_{kj}^2`$, i.e.

$$\delta m_{k}^{}{}_{}{}^{2}<|\mathrm{\Delta }m_{kj}^2|\frac{|\mathrm{\Delta }m_{kj}^2|}{2\sqrt{2}E\sigma _p}>1,$$
(4.25)
the neutrino $`\nu _j`$ is not produced or detected888 The smallness of $`\sigma _p`$ required to satisfy the condition (4.25) implies that $`\nu _k`$ and $`\nu _j`$ cannot be produced or detected in the same process, because their momentum difference is larger than the momentum uncertainty, or their energy difference is larger than the energy uncertainty (for extremely relativistic neutrinos the energy uncertainty is practically equal to the momentum uncertainty $`\sigma _p`$). In this case only one of the two massive neutrinos is produced or detected.  and the interference of $`\nu _k`$ and $`\nu _j`$ is not observed.

The localization term in the oscillation probability (4.19) automatically implements Kayser’s mechanism. Indeed, it can be written as

$$2\pi ^2\xi ^2\left(\frac{\sigma _x}{L_{kj}^{\mathrm{osc}}}\right)^2=\xi ^2\left(\frac{\mathrm{\Delta }m_{kj}^2}{4\sqrt{2}E\sigma _p}\right)^2,$$
(4.26)
with the momentum uncertainty

$$\frac{1}{\sigma _p^2}=4\sigma _x^2=\frac{1}{\sigma _{p}^{\mathrm{P}}{}_{}{}^{2}}+\frac{1}{\sigma _{p}^{\mathrm{D}}{}_{}{}^{2}}.$$
(4.27)
If the condition (4.25) for neutrino mass measurement is satisfied, the localization term (4.26) suppresses999 This argument assumes that $`\xi `$ is not too small. As discussed in Ref. , in realistic cases $`\xi `$ is a number of order one, but, at least in principle, there could be some cases in which $`\xi `$ is very small, or even zero. In these cases Kayser’s mechanism cannot be implemented in the simplified quantum mechanical framework presented here. However, it has been shown in Ref.  that a wave packet quantum field theoretical derivation of the neutrino oscillation probability produces additional terms, independent from $`\xi `$ or similar quantities, which suppress oscillations when the condition (4.25) is satisfied. I would like to thank M. Beuthe for a illuminating discussion about this problem.  the interference of $`\nu _k`$ and $`\nu _j`$.

In the above derivation of neutrino oscillations the wave-packet description of massive neutrinos is crucial in order to allow the integration over time of the time-dependent density matrix operator $`\widehat{\rho }_\alpha (x,t)`$ in Eq. (4.11), which leads to the time-independent density matrix operator $`\widehat{\rho }_\alpha (x)`$ in Eq. (4.12). Of course, as emphasized in Ref. ,

> A single, given density matrix can arise in different ways, especially when incoherence is involved.

As already noticed in Ref.  the density matrix operator $`\widehat{\rho }_\alpha (x)`$ can also be generated through an appropriate incoherent average over the energy spectrum. The same is obviously true for the oscillation probability (4.19).

In order to illustrate this point, let us consider the simplest case of two-neutrino mixing. The wave packet oscillation probability (4.19) becomes

$$P_{\nu _\alpha \nu _\beta }(L)=\frac{1}{2}\mathrm{sin}^22\vartheta \left\{1\mathrm{cos}\left(\frac{\mathrm{\Delta }m^2L}{2E}\right)\mathrm{exp}\left[\left(\frac{L}{L^{\mathrm{coh}}}\right)^22\pi ^2\xi ^2\left(\frac{\sigma _x}{L^{\mathrm{osc}}}\right)^2\right]\right\},$$
(4.28)
with the oscillation and coherence lengths $`L^{\mathrm{osc}}`$ and $`L^{\mathrm{coh}}`$ given, respectively, by Eqs. (4.20) and (4.21) with $`\mathrm{\Delta }m_{kj}^2=\mathrm{\Delta }m^2`$ ($`\mathrm{\Delta }m^2`$ is the squared-mass difference and $`\vartheta `$ is the mixing angle; see the reviews in Refs. ). An incoherent average of the probability in the plane-wave approximation over a gaussian energy spectrum with width $`\sigma _E`$ reads

$$P_{\nu _\alpha \nu _\beta }^{\mathrm{incoh}}(L)=\frac{1}{2}\mathrm{sin}^22\vartheta \left\{1\frac{\mathrm{d}E^{}}{\sqrt{2\pi }\sigma _E}\mathrm{cos}\left(\frac{\mathrm{\Delta }m^2L}{2E^{}}\right)\mathrm{exp}\left[\left(\frac{\left(EE^{}\right)^2}{2\sigma _E^2}\right)^2\right]\right\}.$$
(4.29)
Figure 1 shows the values of the probabilities (4.29) and (4.28) as functions of the source-detector distance $`L`$ for $`\mathrm{\Delta }m^2=2.5\times 10^3\text{ }\mathrm{eV}^2`$, $`\mathrm{sin}^22\vartheta =1`$, $`E=10\text{ }\mathrm{GeV}`$, and $`\sigma _E=\sigma _p=1\text{ }\mathrm{GeV}`$. One can see that it is hard to distinguish the dashed line obtained with Eq. (4.29) from the solid line obtained with Eq. (4.28).

This result does not mean that wave packet effects are irrelevant101010 Let us clarify here that with “wave packet effects” we mean effects due to the momentum and energy uncertainty of a single process associated with the interruption of the emitted wave train caused by decay (natural linewidth) or collisions (collision broadening) discussed in Section 2. These effects are usually not taken into account in the calculation of neutrino energy spectra (see, for example, Ref. ). . It means that in practice one can calculate with reasonable approximation the decoherence of oscillations due to wave packet effects either taking into account the momentum spread in the calculation of the amplitude, as done in the derivation of Eq. (4.28), or averaging the probability over the same energy spread, as done in the derivation of Eq. (4.29)111111 One must only be careful to notice that the energy distribution in the incoherent average of the probability must be normalized to one, whereas the squared modulus of the momentum distribution of the wave packets must be normalized to one. . As emphasized in Ref. , it also means that if one does not have control on both the production and detection processes, one cannot know if oscillations are suppressed because of incoherent averaging over energy of different microscopic processes or because of decoherence due to the separation of wave packets. However, having good control of both the production and detection processes, it may be possible in the far future to reduce the causes of incoherent broadening of the energy spectrum and prove experimentally that oscillations can be suppressed because of wave packet effects (in practice using an approach similar to the one adopted for the measurement of the natural linewidth and collision broadening of atomic lines discussed in Section 2).

Let us finally mention that the fact that coherent and incoherent stationary beams are indistinguishable without a theoretical analysis is well known in optics (see Section 7.5.8 of Ref. ), in neutron interferometry  and in general stationary particle beams .

## 5 Energy of massive neutrinos

In the recent Ref.  the claim that different massive neutrinos must have the same energy , which was refuted in Refs. , has been renewed. If such claim were correct, it would mean that, at least in a quantum mechanical framework, a wave packet description of neutrinos is not necessary. Indeed, in this case one can construct a time-independent density matrix operator from the plane wave state

$$|\stackrel{~}{\nu }_\alpha (x,t)=\underset{k}{}U_{\alpha k}^{}e^{ip_kxiEt}|\nu _k,$$
(5.1)
with $`p_k=\sqrt{E^2m_k^2}`$:

$$\widehat{\stackrel{~}{\rho }}_\alpha (x)=|\stackrel{~}{\nu }_\alpha (x,t)\stackrel{~}{\nu }_\alpha (x,t)|=\underset{k,j}{}U_{\alpha k}^{}U_{\alpha j}\mathrm{exp}\left[i\left(p_kp_j\right)x\right].$$
(5.2)
Using this density matrix operator one can easily derive the standard oscillation probability without need of neutrino wave packets. Therefore, in this section we consider recent claims in favor of an equal energy of different massive neutrinos and show that they are faulty.

The author of Ref.  wrote:

> …states with different ENERGIES ARE NEVER COHERENT in any realistic experiment. States of the same energy and different momenta can be coherent, but may not be.

The reason offered by the author of Ref.  is

> The usual detector is a nucleon, which changes its state after absorbing a neutrino and emitting a charged lepton, and is initially either in an energy eigenstate or in a statistical mixture in thermal equilibrium with its surroundings. No neutrino detector has ever been prepared in a coherent mixture of energy eigenstates and no such detector has been proposed for future experiments.

First, one must note that in terminology of Ref.  “detector” is not a macroscopic device, but a nucleon hit by a neutrino. Then, the claim in Ref.  is easily proved to be wrong using the arguments presented in Section 2, where we have shown that all kinds of waves, classical or quantistic, can be emitted or absorbed only as wave packets, i.e. as coherent superpositions of plane waves with different frequency (energy). If the claim in Ref.  were correct, photons, as well as other particles, could never be emitted or absorbed. The claim that a nucleon should be initially in an energy eigenstate has been shown to be incorrect in Section 3. Furthermore, it is a mystery why all the claims on energy in Ref. , if they were true, do not apply also to momentum.

Also the author of Ref.  claimed that different massive neutrinos have the same energy:

> *Energy or momentum?* In mixing problems, where we have to deal with linear combinations of particles of different mass, the question comes up as to whether one should deal with states of the same energy or the same momentum. Since as stated above, for stationary conditions we are to perform the calculation as an incoherent sum over energies, we have given the answer “energy”. Evidently, for stationary problems it is most natural to use stationary wave functions $`e^{iEt}`$.

The meaning of this paragraph seems to be that different massive neutrinos are supposed to manage to have equal energy in order to satisfy the calculational needs of some theoretical physicists. Therefore, we cannot consider it a proof of anything.

A simple way to show the absurdity of the claims that massive neutrinos must have equal energies is to consider the transformation properties of energy and momentum from one inertial frame to another, following the discussion presented in Ref. .

If the arguments presented in Ref.  were correct, in order to be produced and detected coherently massive neutrinos should necessarily have exactly equal energies in the inertial reference frames of both the production and detection processes (which presumably coincide with the rest frames of the initial nucleons in the production and detection processes, in the approach of Ref. ; let us skip the useless discussion of the determination of the relevant inertial reference frame in neutrino-electron scattering). Such a requirement would imply that the production and detection processes have to be at rest in exactly the same inertial system, because, as shown in Ref.  the equal-energy requirement is not Lorentz invariant.

Therefore if the arguments presented in Ref.  were correct it would mean that in practice neutrino oscillations are not observable, since it is practically impossible to suppress random thermal motion in the source and detector. Certainly, oscillations of solar and atmospheric neutrinos would be impossible. Let me remind that in this discussion we are concerned with the oscillatory terms in the flavor transition probability (4.19), due to the interference of different massive neutrinos, which require coherence. If these terms are suppressed, it is still possible to measure the distance-independent and energy-independent flavor-changing probability

$$P_{\nu _\alpha \nu _\beta }=\underset{k}{}|U_{\alpha k}|^2|U_{\beta k}|^2,$$
(5.3)
or the distance-independent flavor-changing probability calculated in Ref.  if matter effects are important . According to the results of solar neutrino experiments , the solar neutrino problem is very likely due to neutrino oscillations in matter, in which the distance-independent flavor-changing probability calculated in Ref.  is relevant. However, recently the KamLAND experiment  observed indications in favor of distance and energy dependent oscillations in vacuum of reactor antineutrinos due to the same mass-squared difference responsible of solar neutrino oscillations. Distance and energy dependent flavor-changing transitions have been also observed in the Super-Kamiokande , Soudan 2  and MACRO  atmospheric neutrino experiments and in the long-baseline experiment K2K . These experimental evidences show without any doubt that the arguments presented in Ref.  in favor of an equal energy for different massive neutrinos are wrong.

Indeed, if we consider for example atmospheric neutrinos produced by decays in the atmosphere of highly energetic pions and muons, the rest frames of the source and detector are very different and different massive neutrinos cannot have the same energy in both frames. If the arguments presented in Ref.  were correct different massive neutrinos could not be either produced or detected coherently, in contradiction with the observed oscillations. The same reasoning applies to reactor and accelerator neutrinos, remembering that according to the arguments presented in Ref.  coherence is possible only if the source and detector rest frames coincide *exactly* (which, by the way, is an absurd concept, if velocity is a real-valued quantity).

The only rigorous and correct calculation which could be erroneously interpreted121212 Indeed, I made such mistake in the first version of this paper appeared in the hep-ph electronic archive. I am deeply indebted to M. Beuthe for his comments that helped to correct this mistake and my wrong criticisms of Refs. .  as a proof in favor of the equal energy assumption has been presented in Ref.  (see also Section 2.5 of Ref. ). The argument can be illustrated with the density matrix model of neutrino oscillations presented in Section 4. Instead of writing the density matrix operator (4.10) in the form (4.11) using the wave packets (4.8) in which the integration over the momentum $`p`$ in Eq. (4.3) has been already performed, one could write the density matrix operator (4.10) using the wave functions (4.3) as

$$\widehat{\rho }_\alpha (x,t)=\underset{k,j}{}U_{\alpha k}^{}U_{\alpha j}\frac{\mathrm{d}p\mathrm{d}p^{}}{2\pi }\psi _k(p)\psi _j(p^{})e^{i(pp^{})xi(E_k(p)E_j(p^{}))t}|\nu _k\nu _j|.$$
(5.4)
As in Section 4, the time-independent density matrix operator that describes a stationary beam is given by the average over time of $`\widehat{\rho }_\alpha (x,t)`$. The integral over $`t`$ of $`\mathrm{exp}[i(E_k(p)E_j(p^{}))t]`$ yields a $`\delta (E_k(p)E_j(p^{}))`$, which means that “interference occurs only between wave packet components with the same energy”  (see Refs. for similar results in electron interferometry and kaon oscillations). This means that the time average of $`\widehat{\rho }_\alpha (x,t)`$ is equivalent to an appropriate incoherent average over the contributions of the single-energy components of the wave packets. This is the reason why coherent and incoherent stationary beams cannot be distinguished without a theoretical analysis, as discussed in Section 4.

Let us emphasize, however, that this argument does not mean that different massive neutrinos have the same energy, or that wave packets are unnecessary. The wave packet nature of massive neutrinos is necessary for the existence of components with the same energy that produce the observable interference, and the average energies of the wave packets of different massive neutrino are in general different.

## 6 Conclusions

Starting from the well-known fact that “no source of waves vibrates indefinitely”, we have argued that neutrinos, as all other particles, are naturally described by wave packets. This is in agreement with the well known fact that in quantum theory localized particles are described by wave packets. Even when neutrinos are produced or detected through interactions with particles in bound states, they are described by wave packets, because bound states in interaction with an external field do not have a definite energy.

Since the production and detection processes in neutrino oscillation experiments are localized in space-time, a consistent description of neutrino oscillations requires a wave packet treatment. In particular, the wave packet character of massive neutrinos is crucial for the derivation of neutrino oscillations in space, because the group velocity establishes a connection between space and time and allows the time average of the space and time dependent probability (or density matrix).

Finally, we have shown that the claimed arguments in favor of an equal energy of massive neutrinos (which could allow a quantum mechanical derivation of neutrino oscillations without wave packets) are in contradiction with well known physical laws and phenomena.

## Acknowledgements

I would like to thank M. Beuthe and C. Cardall for very interesting and useful comments on the first version of this paper appeared in the hep-ph electronic archive.