format: text

text:
# Quantum Disorder and Quantum Chaos in Andreev Billiards

## I Introduction

Condensed matter physics usually deals with quantum disorder. Quantum disordered systems contain random impurities with the characteristic size smaller than the electron Fermi wavelength $`\lambda _\mathrm{F}`$. In these systems, electrons travel between scattering off impurities along classical trajectories, but the scattering is described by quantum mechanics. Averaging over impurity configurations leads to quantum effects in disordered systems.WL

Recently, new quantum systems - quantum dots - were created.qd  In quantum dots the size of irregularities is much larger than the electron wavelength, therefore electron motion may be determined by a semiclassical theory.

In the semiclassical limit certain thermodynamic and kinetic quantities can be calculated by expansion in powers of the Plank’s constant $`\mathrm{}`$.LL  This expansion corresponds to a series in powers of $`\lambda _\mathrm{F}/L`$, where $`L`$ is the characteristic length scale of the system and $`\lambda _\mathrm{F}`$ is the electron Fermi wave length. For other quantities such semiclassical description is justified only for times shorter than the Ehrenfest time $`t_\mathrm{E}`$. At longer times a crossover from semiclassical to quantum description of these quantities takes place. \[EhrenfestE  demonstrated that the semiclassical description of motion of a quantum particle is justified for short times.\]

What is the Ehrenfest time in classically chaotic systems, where the classical motion is governed by exponential divergence of trajectories? It was shown in Refs. \[LO, ; Z, ; AL1, \], that the Ehrenfest time logarithmically depends on $`\mathrm{}`$. Indeed, in semiclassical theory, a particle is represented by a wave packet - a superposition of particle wave functions. The size of the wave packet cannot be smaller than the electron wavelength $`\lambda _\mathrm{F}`$ and cannot be larger than the characteristic scale $`L`$ of variations of the potential energy. In chaotic systems the size of any wave packet exponentially increases in time. Thus, the wave packet increases from its initial size $`\lambda _\mathrm{F}`$ to size $`d(t)=\lambda _\mathrm{F}\mathrm{exp}(\lambda t)`$ as time $`t`$ increases and expands to scale of the system $`L`$ after the Ehrenfest time $`t_\mathrm{E}=\lambda ^1\mathrm{ln}L/\lambda _\mathrm{F}`$, where $`\lambda `$ is a typical value of the Lyapunov exponent.

We refer to a particle motion at time longer than the Ehrenfest time as a regime of quantum chaos. Quantum chaos affects different phenomena of condensed matter physics. In a metal with random long-range potential the weak localization correction to the a.c.-conductance at frequency $`\omega `$ is absent at high frequencies, $`\omega t_\mathrm{E}1`$.AL1  Nonetheless, at lower frequencies, $`\omega t_\mathrm{E}1`$, the weak localization correction to the conductance reaches its universal quantum limit.

The statistics of electron levels is correctly described by the Gutzwiller trace formula at large energies, $`\epsilon t_\mathrm{E}1`$. At small energies, $`\epsilon t_\mathrm{E}1`$, energy levels obey the Wigner–Dyson statistics. At intermediate energies, $`\epsilon t_\mathrm{E}1`$, the weak localization effects are important for understanding of the level statistics.AL2

The shot noise through a chaotic quantum dot is absent, if the electron escape time $`\tau _{\mathrm{esc}}`$ from the dot is shorter than the Ehrenfest time, $`\tau _{\mathrm{esc}}t_\mathrm{E}`$. In the opposite limit of long $`\tau _{\mathrm{esc}}t_\mathrm{E}`$, the shot noise reaches its universal value.AAL ; SS

The semiclassical method in theory of superconductivity was developed by de Gennes and Tinkham dGT  and Shapoval Shap  for the description of the diffusive electron scattering at surfaces. The applicationsdGT ; Shap ; dGbook  of the semiclassical method gave a basis for a beliefGL ; LOF  that in order to obtain electron properties of a superconductor, first, electron states in the superconductor can be found as solutions of semiclassical equations and, then, these solutions can be averaged over different classical trajectories. However, a more detailed analysisLO  demonstrated that in certain problems \[*e.g.* in calculations of the magnetic penetration depth as a function of an applied magnetic field\], the semiclassical method does not work, and the averaging over classical trajectories gives a different result than that obtained by the standard technique of averaging over impurities.AGD

According to ref. \[LO, \], even for a system with smooth long-range disorder, $`L\lambda _\mathrm{F}`$, the semiclassical method is applicable only at time smaller than the Ehrenfest time $`t_\mathrm{E}`$. Paper \[LO, \] explains the discrepancy between semiclassical and quantum methods in theory of superconductivity: within the semiclassical approach electrons and holes move in opposite directions along identical trajectories inside a superconductor. Because quantum diffraction switches electrons and holes between different trajectories, such semiclassical description breaks down at times larger than the Ehrenfest time $`t_\mathrm{E}`$. The Ehrenfest time $`t_\mathrm{E}`$ and the Lyapunov exponent $`\lambda `$ were calculated for electron systems with smooth impurities of the size $`L\lambda _\mathrm{F}`$ in ref. \[LO, \].

In this paper we study the density of states in an Andreev billiard. An Andreev billiard is a small chaotic metal grain or semiconductor quantum dotqd  connected to a superconductor through a (SN) contact. Electron reflection at the contact is described by the Andreev reflection mechanism,Arefl  when an electron is reflected as a hole, moving in the opposite direction. The density of states in Andreev billiards was studied both by semiclassicalLN ; ScB ; MBFB  and quantum methods.MBFB  The semiclassical method gives an exponentially small but non-zero tail of the density of states for any finite energy, see refs. \[LN, ; ScB, \]. On the contrary, quantum calculationsMBFB  within random matrix theory predict a gap in the energy spectrum.fn1  \[Both curves of the density of states are presented below in Sec. III.2.2.\]

The difference between semiclassical and random matrix results for the density of states has attracted theoretical interest to this problem.TSA ; AB ; Bnew  Based on the approach of refs. \[AL1, ; AL3, \], the authors of ref. \[TSA, \] gave a qualitative explanation of the difference between the semiclassical and the random matrix methods. They argued that the crossover happens at energy of the order of the inverse Ehrenfest time, $`1/\tau _\mathrm{E}`$. Papers \[AB, ; Bnew, \] presented quantitative results for the energy gap and we discuss these results in Sec. V.

We show that in both disordered and chaotic Andreev billiards the density of states has a gap and a square root singularity above the gap. We study a disordered Andreev billiard, containing isotropic point scatters, characterized by arbitrary mean free path $`l_0`$. We investigate the crossover from the semiclassicalLN ; ScB  to quantum regimeMBFB  as parameter $`l_0/l_{\mathrm{esc}}`$ decreases (here $`l_{\mathrm{esc}}=\tau _{\mathrm{esc}}v_\mathrm{F}`$ is the average length of trajectories in the dot and $`v_\mathrm{F}`$ is the Fermi velocity). We further demonstrate that the random matrix result of ref. \[MBFB, \] for reflectionless contact is applicable only at intermediate impurity concentration, when the mean free path $`l_0`$ is smaller than a typical length of trajectories in the dot, $`v_\mathrm{F}\tau _{\mathrm{esc}}`$, and longer than size $`b`$ of the SN contact, $`bl_0v_\mathrm{F}\tau _{\mathrm{esc}}`$. As the mean free path becomes shorter than the SN contact size, $`l_0b`$, we find the density of states, described by a solution of the Usadel equation.fn2

For the clean chaotic billiards, we apply general methods developed in refs. \[AL1, ; AL2, ; AL3, ; AAL, \] to study the density of states. We investigate the crossover from the semiclassical limitLN ; ScB , $`t_\mathrm{E}/\tau _{\mathrm{esc}}\mathrm{}`$, to a random matrix limit,MBFB  $`t_\mathrm{E}/\tau _{\mathrm{esc}}0`$, and calculate the energy gap $`E_\mathrm{g}`$ as a function of $`t_\mathrm{E}/\tau _{\mathrm{esc}}`$. We find that the energy dependence of the density of states has an oscillating component if the Ehrenfest time is finite, see also ref. \[AB, \].

We demonstrate that the density of states has different shapes for systems with different values of the strength of disorder or chaos. In all cases, the density of states has a finite energy gap. The corresponding curves of the density of states are different and cannot be transformed into each other by change of energy scales. Particularly, we consider the high energy asymptote of the density of states, $`\epsilon \tau _{\mathrm{esc}}1`$, and show that the crossover between the semiclassical and quantum regimes can be understood as a suppression of repetitive trajectories.

This paper is organized as follows. In the next section we present basic parameters of the model and introduce the Eilenberger equation. In section III we consider systems with quantum disorder. In section IV we study the effect of quantum chaos on the density of states. Section V contains discussion and conclusions.

## II Model

We consider a small $`d`$-dimensional metal grain connected to a superconductor, see Fig. 1. The characteristic size of the grain is $`L`$ and its volume is $`𝒱_\mathrm{g}L^d`$, with $`d`$ being dimension of a system (in two dimensional system $`𝒱_\mathrm{g}L^2`$ stands for the dot’s area). We assume that the size $`b`$ of the contact between the metal grain and the superconductor is small, so that $`Lb`$, and denote the $`(d1)`$ dimension area of the contact by $`S_\mathrm{c}b^{d1}`$. We also include a random elastic potential into our model, which determines the effect of impurities or irregularities inside the grain and at its boundaries. The length scale between consequent scatterings off this potential is referred to as the mean free path $`l_0=v_\mathrm{F}\tau _0`$ and $`\tau _0`$ is the mean free time. According to the relation between mean free path $`l_0`$ and the size of the system $`L`$ we can distinguish ballistic, $`l_0L`$, or diffusive, $`l_0L`$, limits.

An important energy parameter for a metal grain is the Thouless energy, which is inversely proportional to the ergodic time $`\tau _{\mathrm{erg}}`$. The ergodic time is the characteristic time scale which determines how long it takes for an electron to reach boundaries of the system from arbitrary point in the grain, i.e. to explore the whole phase space of the system. The ergodic time is different in ballistic and diffusive systems. In ballistic systems, the ergodic time is determined by time between collisions with the system boundaries, $`\tau _{\mathrm{erg}}L/v_\mathrm{F}`$, on the other hand in the diffusive regime ergodic time is defined by the diffusion time between boundaries, $`\tau _{\mathrm{erg}}L^2/D`$, where $`D=v_\mathrm{F}l_0/d`$ is the diffusion constant.

Another characteristic time of the system is the escape time $`\tau _{\mathrm{esc}}`$, which determines how long electrons spend in the grain before colliding with the contact. This time can be estimated as the ratio of the total boundary area and area of the contact, multiplied by the ergodic time. Indeed, the ergodic time determines time between collisions with boundaries, and the areas ratio determines the average number of collisions with closed boundaries before a collision with the contact takes place. We have $`\tau _{\mathrm{esc}}\tau _{\mathrm{erg}}(L/b)^{d1}`$. We quantify the escape time from the dot for a two dimensional system by considering the volume $`\mathrm{\Gamma }`$ of the phase space of the grain. The integral over the phase space of the dot can be calculated using two different approaches. First, directly integrating over the direction of the momentum $`𝐧`$ and the volume of the grain, we obtain $`\mathrm{\Gamma }=2\pi 𝒱_\mathrm{g}`$. On the other hand, every point in the phase space is uniquely determined by a trajectory, starting at the superconductor normal metal contact (SN contact). The integral over the phase space in this case can be represented as

$$\mathrm{\Gamma }=_{\pi /2}^{\pi /2}𝑑\theta _{S_\mathrm{c}}\mathrm{cos}\theta l(\theta ,b)𝑑b=2l_{\mathrm{esc}}S_\mathrm{c},$$
(1)
and we obtain the escape time $`\tau _{\mathrm{esc}}=l_{\mathrm{esc}}/v_\mathrm{F}`$:

$$\gamma _{\mathrm{esc}}=\frac{1}{\tau _{\mathrm{esc}}}=\frac{v_\mathrm{F}S_\mathrm{c}}{\pi 𝒱_\mathrm{g}},$$
(2)
where $`𝒱_\mathrm{g}`$ is the area of two dimensional metal. Below we will usually use the escape rate $`\gamma _{\mathrm{esc}}=1/\tau _{\mathrm{esc}}`$.

We emphasize, that the escape time $`\tau _{\mathrm{esc}}`$, Eq. (2), based on the trajectory calculation coincides with the escape time, defined in the random matrix theory as $`\mathrm{}/\tau _{\mathrm{esc}}=N\delta _1/2\pi `$, where $`N`$ is the number of open channels in the contact, $`N=2b/\lambda _\mathrm{F}`$, $`\delta _1=2\pi \mathrm{}^2/m_\mathrm{e}S_\mathrm{g}`$ is the mean level spacing, $`\lambda _\mathrm{F}`$ is the Fermi wave length, $`m_\mathrm{e}`$ is the electron effective mass and $`S_\mathrm{g}`$ is the area of the dot. In this paper we assume $`\delta _1`$ to be negligibly small and $`N`$ to be large.

The superconductivity order parameter $`\mathrm{\Delta }`$ is constant inside the superconducting lead and vanishes fast at the contact with the metal grain. We consider the limit of strong superconductivity, $`\mathrm{\Delta }\gamma _{\mathrm{esc}}`$, and study the density of states at energy $`\epsilon \gamma _{\mathrm{esc}}`$ much smaller than the energy gap $`\mathrm{\Delta }`$ in the superconductor, $`\epsilon \mathrm{\Delta }`$.

The electron Hamiltonian is the sum of kinetic and potential terms:

$$H=\frac{𝐩^2}{2m}\epsilon _\mathrm{F}+V_\mathrm{s}(𝐫)+V_\mathrm{q}(𝐫),$$
(3)
where the potential term $`V_\mathrm{s}(𝐫)+V_\mathrm{q}(𝐫)`$ is divided into semiclassical smooth $`V_\mathrm{s}`$ part, which describes the effect of the grain boundaries, and the quantum part, $`V_\mathrm{q}`$, which represents impurity scattering: $`V_q(𝐫)=_iV_i\delta (𝐫𝐑_i).`$ We write the full Hamiltonian in the Gorkov-Nambu space as

$$\widehat{H}_{\mathrm{GN}}=H\widehat{\tau }_z+\widehat{\mathrm{\Delta }}(𝐫),$$
(4)
where the term $`\widehat{\mathrm{\Delta }}(𝐫)=\widehat{\tau }_x\mathrm{\Delta }\stackrel{~}{\theta }(𝐫)`$ takes into account superconductor pairing, $`\widehat{\tau }_{x,z}`$ are the Pauli matrices and $`\stackrel{~}{\theta }(𝐫)=1`$ inside the superconductor and vanishes otherwise. This approximation for the coordinate dependence of the order parameter is justified in the limit of strong superconductivity, $`\mathrm{\Delta }\epsilon `$.

We introduce the Green’s function as an inverse function of the Hamiltonian:

$$\left(\epsilon \widehat{H}_{\mathrm{GN}}\right)\widehat{G}(\epsilon ,𝐫,𝐫^{})=\delta (𝐫𝐫^{})\widehat{1}$$
(5)
and perform the Wigner transformation

$$\widehat{G}(\epsilon ,𝐑,𝐩)=e^{i\mathrm{𝐩𝐫}}\widehat{G}(\epsilon ,𝐑+𝐫/2,𝐑𝐫/2)𝑑𝐫.$$
(6)
Next, we change coordinates in the phase space from $`𝐫`$ and $`𝐩`$ to $`\xi `$ and $`𝐱`$, where $`\xi =H(𝐫,𝐩)\epsilon _\mathrm{F}`$ and $`𝐱=(𝐫,𝐧)`$ with $`𝐧=𝐩/|𝐩|`$. The reduced Green’s function is obtained by integration over $`\xi `$:

$$\widehat{𝒢}(\epsilon ,𝐱)=\left(\begin{array}{cc}g(\epsilon ,𝐱)& f_+(\epsilon ,𝐱)\\ f_{}(\epsilon ,𝐱)& g(\epsilon ,𝐱)\end{array}\right)=\frac{i}{\pi \nu }\widehat{G}(\epsilon ,\xi ,𝐱)𝑑\xi ,$$
(7)
where $`\nu `$ is the density of states in the grain per spin degree of freedom in the absence of superconductor. The semiclassical Green’s function allows us to calculate electron density of states in the system as an integral over the phase space of electrons in the grain:

$$𝒩(\epsilon )=\nu \mathrm{Re}\mathrm{Tr}\left\{[\widehat{\tau }_z,\widehat{𝒢}(\epsilon ,𝐱)]\right\}\frac{d𝐱}{\mathrm{\Omega }_d}$$
(8)
Here $`[,]`$ stands for a commutator. An element of the phase space $`d𝐱`$ is $`d𝐱=d𝐫d\mathrm{\Omega }_𝐧`$; $`d\mathrm{\Omega }_𝐧`$ denotes a measure on a sphere of a unit radius, representing electron momentum direction $`𝐧`$, and $`\mathrm{\Omega }_d`$ is the total area of a sphere of a unit radius in $`d`$dimensional space, $`S_2=2\pi `$ and $`S_3=4\pi `$.

The reduced Green’s function $`\widehat{𝒢}(\epsilon ,𝐱)`$ satisfies the Eilenberger equation: Eilen ; LO

$`[i\epsilon \widehat{\tau }_z+\widehat{\mathrm{\Delta }}(𝐫),\widehat{𝒢}(\epsilon ,𝐱)]\widehat{𝒢}(\epsilon ,𝐱)=I[\widehat{𝒢}(\epsilon ,𝐱)],`$ (9)
where

$$=\frac{𝐩}{m}\frac{}{𝐫}\frac{V_\mathrm{s}(𝐫)}{𝐫}\frac{}{𝐩}$$
(10)
is the Liouville operator and $`I[\widehat{𝒢}(\epsilon ,𝐱)]`$ is the scattering term. The scattering term takes into account quantum part of the potential $`V_\mathrm{q}`$, introduced in Eq. (3). A general form for the scattering term in the Bohr approximation is:

$$I[\widehat{𝒢}(𝐫,𝐧)]=[\widehat{\mathrm{\Sigma }}(𝐫,𝐧),\widehat{𝒢}(𝐫,𝐧)],$$
(11a)
where $`\widehat{\mathrm{\Sigma }}(𝐫)`$ stands for the integral
$`\widehat{\mathrm{\Sigma }}(\epsilon ,𝐫,𝐧)`$ $`=`$ $`\left(\begin{array}{cc}\sigma _0(\epsilon ,𝐫,𝐧)& \sigma _+(\epsilon ,𝐫,𝐧)\\ \sigma _{}(\epsilon ,𝐫,𝐧)& \sigma _0(\epsilon ,𝐫,𝐧)\end{array}\right)`$
$`=`$ $`\pi \nu n_\mathrm{i}{\displaystyle \widehat{𝒢}(𝐫,𝐧^{})|V_\mathrm{q}(𝐧^{}𝐧)|^2\frac{d\mathrm{\Omega }_𝐧^{}}{\mathrm{\Omega }_d}}`$
over direction of momentum $`𝐧`$ at the Fermi surface, $`V_\mathrm{i}(𝐧𝐧^{})`$ represents the scattering amplitude off a single impurity, the impurity concentration is $`n_\mathrm{i}`$.

A solution of the Eilenberger equation, Eq. (9), is restricted by the constraints:

$$\widehat{𝒢}^2(\epsilon ,𝐱)=\widehat{1},\mathrm{tr}\widehat{𝒢}(\epsilon ,𝐱)=0,$$
(12)
which represent the normalization conditions for $`\widehat{𝒢}(\epsilon ,𝐱)`$.

Inside the superconductor at distance $`\xi `$ away from the SN contact the Green function is given by

$$\widehat{𝒢}(\epsilon ,\xi )=\widehat{\tau }_x+\alpha _+e^{2\mathrm{\Delta }\xi }\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right)+\alpha _{}e^{2\mathrm{\Delta }\xi }\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right),$$
(13)
where $`\xi <0`$ for electrons moving in the superconductor towards the grain and $`\xi >0`$ for electrons moving in the superconductor away from the grain. The requirement that the Green function inside the superconductor relaxes to $`\widehat{𝒢}(\epsilon ,\xi )=\widehat{\tau }_x`$, gives the following form of the Green function $`\widehat{𝒢}_\mathrm{s}(\epsilon ,𝐱)`$ at the SN contact inside the superconductor:LN

$$\widehat{𝒢}_\mathrm{s}(\epsilon ,𝐱)=\widehat{\tau }_x+\alpha \left(\begin{array}{cc}1& \varsigma \\ \varsigma & 1\end{array}\right),$$
(14)
where $`\varsigma =+1`$ for electrons incoming to the grain, *i.e.* $`\mathrm{𝐧𝐧}_\mathrm{s}>0`$, and $`\varsigma =1`$ for outgoing electrons, $`\mathrm{𝐧𝐧}_\mathrm{s}<0`$. Here $`𝐧_\mathrm{s}`$ is a unit vector perpendicular to the superconductor-metal interface and directed to the grain. The second term of Eq. (14) has a specific form, so that it corresponds to a vanishing solution of the Eilenberger equation inside the superconductor.LN  Coefficient $`\alpha `$ is a free parameter. For a reflectionless SN contact the Green function is continuous and Eq. (14) represents the boundary conditions for electron Green functions inside the grain at the SN contact.

We represent a point of the phase space by its coordinate $`𝐫`$ and the direction of the electron momentum $`𝐧`$. Alternatively, we can parameterize the phase space by trajectories, which begin and end at the SN interface. Since only one trajectory goes through point $`𝐫`$ in direction $`𝐧`$, the trajectory and the distance $`\zeta `$ along this trajectory, measured from the starting point, uniquely represent point $`𝐱=(𝐫,𝐧)`$ of the phase space. Below we often use the trajectory parameterization of the phase. A point at the SN contact with $`\mathrm{𝐧𝐧}_s>0`$ has always $`\zeta =0`$, and a point at the contact with $`\mathrm{𝐧𝐧}_s<0`$ has $`\zeta =l`$, where $`l`$ is the length of the corresponding trajectory. We omit the Fermi velocity and imply that the distance $`\zeta `$ is measured in time units. In this case, $`\zeta `$ corresponds to time, which electron travels from the beginning of a trajectory to a given point, and $`l`$ is the total time of motion along this trajectory.

The normalization conditions, Eqs. (12), for the semiclassical Green’s function $`\widehat{𝒢}(\epsilon ,𝐱)`$ leaves only two independent functions, which parameterize the matrix $`\widehat{𝒢}(\epsilon ,𝐱)`$. For certain purposes it is convenient to choose the Riccati parametrization in terms of functions $`a_\pm =a_\pm (\epsilon ,𝐱)`$:Maki

$$\widehat{𝒢}=\frac{1}{1+a_+a_{}}\left(\begin{array}{cc}1a_+a_{}& 2a_+\\ 2a_{}& a_+a_{}1\end{array}\right)$$
(15)
The Eilenberger equations for functions $`a_\pm (\epsilon ,𝐱)`$ acquire the following form:

$`\left(2i\epsilon \right)a_+`$ $`=`$ $`\sigma _{}a_+^2+2\sigma _0a_+\sigma _+,`$ (16a)
$`\left(2i\epsilon +\right)a_{}`$ $`=`$ $`\sigma _+a_{}^2+2\sigma _0a_{}\sigma _{}.`$ (16b)
Here the self energy elements $`\sigma _\pm =\sigma _\pm (\epsilon ,𝐱)`$ and $`\sigma _0=\sigma _0(\epsilon ,𝐱)`$ are defined by Eq. (11).

The boundary conditions, Eq. (14), acquire a simple form in the Riccati parameterization. Using the continuity condition for the semiclassical Green’s function $`\widehat{𝒢}(\epsilon ,𝐱)`$ at the SN interface, and comparing Eqs. (15) and (14) we find that $`a_+(\epsilon ,l)=1`$ at the end of a trajectory of length $`l`$, \[$`\varsigma <0`$ in Eq. (14)\] and $`a_{}(\epsilon ,\zeta =l)`$ is not restrained and reflects a freedom of factor $`\alpha `$ in Eq. (14). Similarly, $`a_{}(\epsilon ,0)=1`$ for incoming electrons, $`\varsigma >0`$, while $`a_+(\epsilon ,\zeta )`$ is a free parameter, see Table 1.

In the limit of strong superconductivity, $`\mathrm{\Delta }\epsilon `$, an electron moves along some trajectory and when it is Andreev reflected as a hole, the hole moves in the opposite direction along the same trajectory. Functions $`a_\pm (\epsilon ,𝐫,𝐧)`$, which describe electron Green function in the Riccati parametrization, Eq. (15), have the meaning of the phase factor, which electron-hole pairs acquire as a result of a single Andreev reflection at the SN contact while they travel from the contact to point $`\{𝐫,𝐧\}`$ (note the boundary conditions in Table 1). We refer to $`a_\pm (\epsilon ,𝐫,𝐧)`$ as wave functions of an electron-hole pairs at point $`\{𝐫,𝐧\}`$. The off-diagonal elements of the Green function $`\widehat{𝒢}(\epsilon ,𝐱)`$ may be represented as a series in $`a_\pm (\epsilon ,𝐫,𝐧)`$ and correspond to the full amplitude, which takes into account multiple Andreev reflections, see Fig. 2. The diagonal element of the Green function may also be written in the form of a series, which gives us the following expression for the electron density of states in the grain:

$`𝒩(\epsilon )`$ $`=`$ $`{\displaystyle \frac{2}{\delta _1}}{\displaystyle \frac{d𝐱}{𝒱_\mathrm{g}\mathrm{\Omega }_d}}`$ (17)
$`\times `$ $`\left(1+2{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}(1)^k\mathrm{Re}\left\{a_+^k(\epsilon ,𝐱)a_{}^k(\epsilon ,𝐱)\right\}\right).`$
Here $`\delta _1=\nu 𝒱_\mathrm{g}`$ is the mean level spacing of electron states in the grain of the volume $`𝒱_\mathrm{g}`$. The $`k`$th term of this series represents the effect of $`k`$ Andreev reflections. We notice, that the first term of this expansion (unity) correspond to electron Green function in the grain without superconductivity at the contact.

In conclusion of this Section, we make the following observation. The scattering term vanishes after integration over the direction of momentum at each point of coordinate space. If we integrate the Eilenberger equation in the form of Eq. (9) over the phase space, we obtain

$$_{S_\mathrm{c}}\mathrm{𝐧𝐧}_s\widehat{𝒢}(\epsilon ,𝐫,𝐧)dS_\mathrm{c}\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}=i\epsilon _{𝒱_\mathrm{g}}[\widehat{\tau }_z,\widehat{𝒢}(\epsilon ,𝐫,𝐧)]𝑑𝐫\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d},$$
(18)
where $`dS_\mathrm{c}`$ is an element of the SN interface and $`𝐧_\mathrm{s}`$ is a unit vector, perpendicular to $`dS_\mathrm{c}`$. We notice that the diagonal part of Eq. (18) gives

$$_{S_\mathrm{c}}\mathrm{𝐧𝐧}_sg(\epsilon ,𝐫,𝐧)dS_\mathrm{c}d\mathrm{\Omega }_𝐧=0.$$
(19a)
The latter equation represents conservation of the electric charge in the grain. The left hand side part of the off-diagonal part of Eq. (18) contains the following difference $`f_+(\epsilon ,l)f_+(\epsilon ,0)`$ of the off-diagonal elements of the Green function at the terminals of trajectory of length $`l`$. From the boundary conditions, defined by Eq. (14), we know that $`f_+(\epsilon ,\zeta )`$ and $`g(\epsilon ,\zeta )`$ are related at the trajectory terminals:
$`f_+(\epsilon ,0)`$ $`=`$ $`1g(\epsilon ,0),`$
$`f_+(\epsilon ,l)`$ $`=`$ $`1+g(\epsilon ,l).`$
These relations allow us to write the following equation for the off-diagonal part of Eq. (18):
$$_{S_\mathrm{c}}|\mathrm{𝐧𝐧}_s|g(\epsilon ,𝐫,𝐧)dS_\mathrm{c}\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}=2i\epsilon _{𝒱_\mathrm{g}}f(\epsilon ,𝐫,𝐧)𝑑𝐫\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}.$$
(19b)
This equation makes a connection between the boundary values of the Green’s function and the integral of the Green’s function over the phase space in the metal grain.

## III Quantum Disorder

In this section we study density of states in the dirty system, when the quantum part of the potential is produced by small (point) impurities. The collision term can be written in the form of Eqs. (11).

The mean free time is defined as the integral over different directions of momentum:

$$\frac{1}{\tau _0}=2\pi \nu n_\mathrm{i}|V_\mathrm{q}(𝐧𝐧^{})|^2\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}.$$
(20)
Note, that the angular dependence of the Green’s function introduces scattering time, different from the mean free path $`\tau _0`$. Particularly, in the Usadel equation angular dependence of the Green’s function on the momentum direction leads to the appearance of the transport mean free time $`\tau _{\mathrm{tr}}`$, see Eq. (22) below.

In this section we assume that the scattering potential has no sharp angle dependence, and we use the mean free time as a characteristic value of the impurity scattering strength.

### III.1 Diffusive contact

In the limit when the scattering is strong, Eilenberger equation reduces to the Usadel equationU  \[we will discuss conditions for applicability of the Usadel equation below\]. In this limit, the Eilenberger Green’s function weakly depends on the direction of the momentum $`𝐧`$ and we represent Green’s function in the form: $`\widehat{𝒢}(\epsilon ,𝐫,𝐧)=\widehat{𝒢}_0(\epsilon ,𝐫)+𝐧\widehat{𝐠}_1(\epsilon ,𝐫)`$, where we assume $`\widehat{𝒢}_0(\epsilon ,𝐫)\mathrm{𝐧𝐠}_1(\epsilon ,𝐫)`$. Taking into account the normalization condition, Eq. (12), we find

$$𝐠_1(\epsilon ,𝐫)=v_\mathrm{F}\tau _{\mathrm{tr}}\widehat{𝒢}(\epsilon ,𝐫)\frac{}{𝐫}\widehat{𝒢}(\epsilon ,𝐫),$$
(21)
where we introduced the transport mean free time $`\tau _{\mathrm{tr}}`$:

$$\frac{1}{\tau _{\mathrm{tr}}}=2\nu n_\mathrm{i}(1\mathrm{𝐧𝐧}^{})|V_\mathrm{q}(𝐧^{}𝐧)|^2\frac{d^d𝐧}{\mathrm{\Omega }_d}.$$
(22)
Substituting $`𝐠_1(\epsilon ,𝐫)`$ into the Eilenberger equation, Eq. (9), we obtain the Usadel equation for a system with time-reversal symmetry:

$$D^2\mathrm{\Theta }=2i\epsilon \mathrm{sin}\mathrm{\Theta }2\mathrm{\Delta }(𝐫)\mathrm{cos}\mathrm{\Theta },$$
(23)
where we parameterize the Green function by $`\mathrm{\Theta }=\mathrm{\Theta }(\epsilon ,𝐫)`$

$$\widehat{𝒢}=\left(\begin{array}{cc}\mathrm{cos}\mathrm{\Theta }& \mathrm{sin}\mathrm{\Theta }\\ \mathrm{sin}\mathrm{\Theta }& \mathrm{cos}\mathrm{\Theta }\end{array}\right)$$
(24)
and $`D=v_\mathrm{F}^2\tau _{\mathrm{tr}}/d`$ is the diffusion coefficient.

We solve the Usadel equation for a $`d`$-dimensional metal grain of volume $`𝒱_\mathrm{g}L^d`$, connected to a superconductor by a small contact of the area $`S_\mathrm{c}`$, see Fig. 1. In the normal region we substitute $`\mathrm{\Delta }=0`$ into Eq. (23) and integrate the resulting equation over the volume $`𝒱_\mathrm{g}`$ of the grain:

$$D_{S_\mathrm{c}}𝐧_\mathrm{s}\mathrm{\Theta }(\epsilon ,𝐬)𝑑S_\mathrm{c}=2i\epsilon _{𝒱_\mathrm{g}}\mathrm{sin}\mathrm{\Theta }(\epsilon ,𝐫)𝑑𝐫,$$
(25)
where the integral in the left hand side is taken over the area of the contact $`S_\mathrm{c}`$, and $`𝐧_\mathrm{s}`$ denotes the normal vector at the SN contact. Since $`\mathrm{\Theta }(𝐫)`$ varies only close to the contact, we neglect variation of $`\mathrm{\Theta }`$ and replace the right hand side of Eq. (25) by $`2\epsilon 𝒱_\mathrm{g}\mathrm{sin}\mathrm{\Theta }_0`$, where $`\mathrm{\Theta }_0`$ is an asymptotic value of $`\mathrm{\Theta }(𝐫)`$ away from the contact.

To calculate the left hand side of Eq. (25) we solve the following equation

$$D^2\mathrm{\Theta }=0,$$
(26)
neglecting the term with energy $`\epsilon `$ in the Usadel equation. This approximation is justified for energy smaller than the characteristic value of the left hand side of Eq. (26), which we estimate as $`D/b^2`$.

We notice that Eq. (26) coincides with equation $`^2\phi (𝐫)=0`$ for an electrostatic potential $`\phi (𝐫)`$ in the problem of calculation of the current

$$I=_{S_\mathrm{c}}𝐧_\mathrm{s}𝐣𝑑S_\mathrm{c}$$

through the contact. We can write the Ohms law $`I=(\phi _\mathrm{c}\phi _0)/R_\mathrm{c}`$ for the current $`I`$, where $`\phi _\mathrm{c}`$ and $`\phi _0`$ are the values of the electric potential $`\phi `$ at the contact and inside the grain, respectively, and $`R_\mathrm{c}`$ is the resistance of the contact. The contact resistance $`R_\mathrm{c}`$ depends on the geometry of the contact. For a circular contact with diameter $`b`$ we have $`\sigma R_\mathrm{c}1/b`$ for a three dimensional system and $`\sigma R_\mathrm{c}\mathrm{log}L/b`$ for a two dimensional system.

Identifying the left hand side of Eq. (25) with the current $`I`$ and the Usadel parameter $`\mathrm{\Theta }(𝐫)`$ with the electric potential $`\phi (𝐫)`$, we reduce the solution of the Usadel equation to the electrostatic problem. Then, using the boundary condition for $`\mathrm{\Theta }=\pi /2`$ at the contact and introducing the value of the parameter $`\mathrm{\Theta }`$ in the grain as $`\mathrm{\Theta }=\pi /2\vartheta `$, we obtain the equation for $`\vartheta `$:

$$\vartheta =2i\frac{\epsilon }{\gamma _\mathrm{c}}\mathrm{cos}\vartheta ,\gamma _\mathrm{c}=\frac{D}{\sigma R_\mathrm{c}𝒱_\mathrm{g}}.$$
(27)
Here we introduced the quantity $`\gamma _\mathrm{c}`$, which has a meaning of the electron escape rate from the grain with a diffusive contact. Because $`\sigma e^2\nu D/\mathrm{}`$, we have $`\gamma _\mathrm{c}(\mathrm{}/e^2R_\mathrm{c})(1/\nu 𝒱_\mathrm{g})=N_{\mathrm{trans}}\delta _1`$, where $`N_{\mathrm{trans}}=\mathrm{}/e^2R_\mathrm{c}`$ is the number of transparent channels in the SN contact \[compare $`\gamma _\mathrm{c}`$ to $`\gamma _{\mathrm{esc}}=N\delta _1/2\pi `$ for a clean contact with $`N`$ channels\].

The density of states in the grain is given by the real part of the Eilenberger Green function

$$𝒩(\epsilon )=\frac{2}{\delta _1}\mathrm{Re}\mathrm{sin}\vartheta .$$
(28)
and vanishes, if $`\vartheta `$ is purely imaginary. For small energy $`\epsilon <E_\mathrm{g}^\mathrm{U}`$, Eq. (27) has onlu imaginary solutions, and the spectrum begins at threshold energy $`E_\mathrm{g}^\mathrm{U}`$, where $`\vartheta `$ acquires a real part. From Eq. (27) we find the gap energy $`E_\mathrm{g}`$

$$E_\mathrm{g}^\mathrm{U}0.331\gamma _\mathrm{c}.$$
(29a)
Above the gap energy $`E_\mathrm{g}^\mathrm{U}`$ the density of states exhibits a square root singularity:
$$𝒩(\epsilon )=\frac{c_\mathrm{U}}{\delta _1}\sqrt{\frac{\epsilon }{E_\mathrm{g}^\mathrm{U}}1}$$
(29b)
and the numerical constant $`c_\mathrm{U}3.62`$. At higher energy the density of states approaches its metal value $`\nu `$, according to the following expression:

$$𝒩(\epsilon )=\frac{2}{\delta _1}\left(1+\frac{\pi ^2}{32}\frac{\gamma _\mathrm{c}^2}{\epsilon ^2}\right).$$
(30)
We plot the density of states, described by Eq. (28) through a solution of Eq. (27), in Fig. 3 by a solid line. For comparison, we present the random matrix resultMBFB  for the case of reflectionless channels by a solid line.fn2

This result is applicable for $`\epsilon D/b^2`$. At larger energies, $`\epsilon D/b^2`$, Eq. (26) and, consequently, Eq. (27) for the parameter $`\vartheta `$ are no longer valid.

It is instructive to compare the result of Eqs. (29) with the solution of Usadel equation in the one dimensional case for wires of length $`L`$ with both ends connected to a superconductor, see refs. \[1dU, ; OSF, \]. The value of the energy gap and density of states above the gap are given by:

$$E_\mathrm{g}3.12\frac{D}{L^2},𝒩(\epsilon )\frac{2.30}{\delta _1}\sqrt{\frac{\epsilon }{E_\mathrm{g}}1}.$$
(31)
We observe that $`E_\mathrm{g}D/L^2`$. In this case the Laplasian term in Eq. (25) is comparable with the energy term and approximation of Eq. (26) is not applicable. The Green’s function is obtained as a result of exact integration of the one-dimensional sine-Gordon equation.1dU ; OSF

Deriving Eq. (21) we have to assume that $`𝐧\widehat{𝐠}_11`$, on the other hand, using Eq. (21) and taking into account that $`\widehat{𝒢}(𝐫)1/b`$, we obtain $`𝐧\widehat{𝐠}_1l_0/b`$. The result of Eqs. (29) and (30) is valid in the limit when $`l_0b`$. Beyond this limit the Usadel equation is not applicable.

For one dimensional systems, the density of states was studied beyond the Usadel equation in refs. \[PBB, \], where an arbitrary relation between the system size $`L`$ and the mean free path $`l_0`$ was considered. The conclusion of refs. \[PBB, \] is that there always exists a finite energy gap in the density of states.

### III.2 Ballistic contact

In this subsection we study the density of states in the limit, when the mean free path is larger than the size of the contact $`l_0b`$. In the intermediate regime, $`bl_0`$, the result depends on specific geometry of the contact, and the solution of the Eilenberger equation may be studied numerically in the spirit of refs. \[PBB, \].

#### III.2.1 Random Matrix Limit

Here we consider the limit of short mean free time $`\tau _0\tau _{\mathrm{esc}}`$. The Green’s function coincides with the Green function at the contact, see Eq. (14) only at small distances $`\zeta \tau _0`$ or $`l\zeta \tau _0`$ along trajectories from the terminals of these trajectories. After several scatterings, at distances larger than the mean free path, the Green’s function averages out and becomes independent from momentum direction and coordinate. We refer to this coordinate independent part of the Green’s function as a zero mode and define it by

$`\widehat{𝒢}_0`$ $`=`$ $`\left(\begin{array}{cc}g_0(\epsilon )& f_0(\epsilon )\\ f_0(\epsilon )& g_0(\epsilon )\end{array}\right)`$ (34)
$`=`$ $`{\displaystyle \frac{1}{1+a_0^2(\epsilon )}}\left(\begin{array}{cc}1a_0^2(\epsilon )& 2a_0(\epsilon )\\ 2a_0(\epsilon )& 1+a_0^2(\epsilon )\end{array}\right).`$ (37)
Here the second equation introduces the Riccati parameterization of the zero mode component in terms of a single function $`a_0(\epsilon )`$, see Eq. (15).

The Green function $`\widehat{𝒢}^\mathrm{T}(\epsilon ,\zeta )`$ near the boundaries ($`\zeta =0,l`$) is:

$`\widehat{𝒢}(\epsilon ,0)`$ $`=`$ $`\widehat{𝒢}^\mathrm{T}(\epsilon ,l)=\left(\begin{array}{cc}g(\epsilon ,0)& f_+(\epsilon ,0)\\ f_{}(\epsilon ,0)& g(\epsilon ,0)\end{array}\right)`$ (40)
$`=`$ $`{\displaystyle \frac{1}{1+a_0(\epsilon )}}(\begin{array}{cc}1a_0(\epsilon )& 2a_0(\epsilon )\\ 2& 1+a_0(\epsilon )\end{array}).`$ (43)
Here we took into account the boundary conditions $`a_{}(\epsilon )=1`$ for incoming electrons and $`a_+(\epsilon )=1`$ for outgoing electrons, see Table 1. The unrestrained values of $`a_+(\epsilon )`$ and $`a_{}(\epsilon )`$ at the SN contact coincide with the zero mode value $`a_0(\epsilon )`$ inside the grain.

At sufficiently small values of the mean free time $`\tau _0`$, we neglect the contribution to the average value of the Green’s function from the parts of trajectories close to the SN contact, where $`\widehat{𝒢}(\epsilon ,\zeta )`$ is given by Eq. (43), since this contribution is as small as $`\gamma _{\mathrm{esc}}\tau _0`$, where $`\gamma _{\mathrm{esc}}`$ is defined by Eq. (2) in terms of the contact area and the grain volume. Then, the integral equation (19b) reduces to

$$\gamma _{\mathrm{esc}}g(\epsilon ,0)=i\epsilon f_0(\epsilon ).$$
(44)
Substituting $`f_0(\epsilon )`$ from Eq. (37) and $`g(\epsilon ,0)`$ from Eq. (43) into Eq. (44), we rewrite Eq. (44) in the form:

$$\gamma _{\mathrm{esc}}\frac{1a_0(\epsilon )}{1+a_0(\epsilon )}=i\epsilon \frac{2a_0(\epsilon )}{1+a_0^2(\epsilon )}.$$
(45)
We substitute $`a_0(\epsilon )=\mathrm{exp}(2i\psi )`$ and reduce Eq. (45) to

$$\mathrm{tan}\psi \mathrm{cos}2\psi =\frac{\epsilon }{\gamma _{\mathrm{esc}}}.$$
(46)
The density of states

$$𝒩(\epsilon )=\frac{2}{\delta _1}\mathrm{Re}[g_0(\epsilon )]=\frac{2}{\delta _1}\mathrm{Im}\mathrm{tan}2\psi $$
(47)
is defined by Eq. (37) in terms of the solution $`a_0(\epsilon )`$ and is finite only if $`\psi `$ has an imaginary part. The function in the left hand side of Eq. (46) has a maximum at certain value $`\psi _\mathrm{m}=\mathrm{arccos}(\sqrt{1+\sqrt{5}}/2)0.45`$ and we define

$`E_\mathrm{g}^{\mathrm{rmt}}`$ $`=`$ $`\mathrm{tan}\psi _\mathrm{m}\mathrm{cos}2\psi _\mathrm{m}\gamma _{\mathrm{esc}}`$ (48a)
$`=`$ $`\sqrt{{\displaystyle \frac{5\sqrt{5}11}{2}}}\gamma _{\mathrm{esc}}0.30\gamma _{\mathrm{esc}}.`$
At small energy $`\epsilon <E_\mathrm{g}^{\mathrm{rmt}}`$, Eq. (46) has real solutions and the density of states vanishes. Above the gap we expand the left hand side of Eq. (46) in $`\psi \psi _\mathrm{m}`$ as $`\upsilon (\psi \psi _\mathrm{m})^2=E_\mathrm{g}^{\mathrm{rmt}}\epsilon `$, with $`\upsilon 2.17`$, and obtain that the density of states has a square root singularity above the gap
$$𝒩(\epsilon )=\frac{c_{\mathrm{rmt}}}{\delta _1}\sqrt{\frac{\epsilon }{E_\mathrm{g}^{\mathrm{rmt}}}1}$$
(48b)
with $`c_{\mathrm{rmt}}=2\sqrt{2+4/\sqrt{5}}3.89`$.

At energies larger than the escape rate $`\gamma _{\mathrm{esc}}`$, the density of states approaches the normal metal value $`2/\delta _1`$, according to

$$𝒩(\epsilon )=\frac{2}{\delta _1}\left(1+\frac{\gamma _{\mathrm{esc}}^2}{2\epsilon ^2}\right).$$
(49)
The plot of the density of states is presented in Fig. 3 by the dashed line. We emphasize that the above result for the density of states, Eqs. (46)-(49) coincides with the result of ref. \[MBFB, \] obtained within random matrix theory.fn1

#### III.2.2 Crossover regime

Here we consider arbitrary mean free time $`\tau _0\tau _{\mathrm{esc}}`$ due to scattering off isotropic impurities. In this case we have the following form of the self energy:

$`\widehat{\mathrm{\Sigma }}(\epsilon )={\displaystyle \frac{1}{2\tau _0}}{\displaystyle \widehat{𝒢}(\epsilon ,𝐫,𝐧)\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}}`$ (50)
i.e. the self energy is proportional to the angle average of the Green’s function at point $`𝐫`$ inside the grain. For small dots, when the Thouless energy $`E_{\mathrm{Th}}`$ is the largest energy scale of the grain, we will make a conjecture, that the averaging over direction of electron momentum is equivalent to averaging over the full phase space of the system:

$`\widehat{𝒢}(\epsilon ,𝐫,𝐧)`$ $`=`$ $`{\displaystyle \frac{1}{𝒱_\mathrm{g}}}{\displaystyle \widehat{𝒢}(\epsilon ,𝐫,𝐧)𝑑𝐫\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}}`$ (51)
$`=`$ $`{\displaystyle \widehat{𝒢}(\epsilon ,𝐫,𝐧)\frac{d\mathrm{\Omega }_𝐧}{\mathrm{\Omega }_d}}.`$
Indeed, different directions of the momentum represent contributions from different points of the grain, since we assume that the Green’s function does not change fast at short parts of trajectories, corresponding to electron motion during the ergodic time $`\tau _{\mathrm{erg}}=1/E_{\mathrm{Th}}`$.

The Green’s function $`\widehat{𝒢}(\epsilon ,𝐱)`$ satisfies the Eilenberger equation in the form:

$$\widehat{𝒢}(\epsilon ,𝐫,𝐧)=[\widehat{H}_{\mathrm{eff}}(\epsilon ),\widehat{𝒢}(\epsilon ,𝐫,𝐧)],$$
(52a)
where
$$\widehat{H}_{\mathrm{eff}}=\left(\begin{array}{cc}G_0& F_0\\ F_0& G_0\end{array}\right)=i\epsilon \widehat{\tau }_z+\frac{1}{2\tau _0}\widehat{𝒢}(\epsilon ,𝐫,𝐧).$$
(52b)
Assuming that $`\widehat{H}_{\mathrm{eff}}(\epsilon )`$ is fixed, Eq. (52a) is linear in $`\widehat{𝒢}(\epsilon ,𝐫,𝐧)`$. We choose the parametrization of the grain’s phase space in terms of the classical trajectories, which go through space point $`𝐫`$ and in the direction of $`𝐧`$. Then, the solution of Eq. (52a) can be found for each trajectory as a function of the coordinate along the trajectory $`\zeta `$, since the Liouville operator has the form $`=/\zeta `$. In this trajectory parametrization of the phase space, the smooth part of the potential $`V_\mathrm{s}`$ is taken into account as the refraction of trajectories.

The solution to Eq. (52a) can be written in terms of three matrices:

$`\widehat{\varrho }_z={\displaystyle \frac{1}{D_0}}\widehat{H}_{\mathrm{eff}};\widehat{\varrho }_\pm ={\displaystyle \frac{1}{2D_0}}\left(\begin{array}{cc}F_0& G_0\pm D_0\\ G_0D_0& F_0\end{array}\right).`$ (55)
Matrices $`\widehat{\varrho }_z`$ and $`\widehat{\varrho }_\pm `$ are related to the Pauli matrices $`\widehat{\tau }_z`$ and $`\widehat{\tau }_\pm =(\widehat{\tau }_x\pm i\widehat{\tau }_y)/2`$ through a rotation in the Gorkov-Nambu space, which diagonalizes the Hamiltonian $`\widehat{H}_{\mathrm{eff}}`$, see Eq. (52b). If $`\widehat{H}_{\mathrm{eff}}`$ is diagonal, matrices $`\widehat{\varrho }_i`$ coincide with matrices $`\widehat{\tau }_i`$. A general solution of Eq. (52a) is

$$\widehat{𝒢}(\epsilon ,\zeta )=\alpha _z\widehat{\varrho }_z+\alpha _+\widehat{\varrho }_+e^{2D_0(\zeta l/2)}+\alpha _{}\widehat{\varrho }_{}e^{2D_0(\zeta l/2)},$$
(56)
where $`D_0^2=G_0^2+F_0^2`$ and coordinate $`\zeta `$ is a length along a trajectory of length $`l`$. Coefficients $`\alpha _z`$ and $`\alpha _\pm `$ are uniquely determined by the boundary conditions, Eq. (14). We have

$`\alpha _z`$ $`=`$ $`{\displaystyle \frac{F_0\mathrm{cosh}D_0l+D_0\mathrm{sinh}D_0l}{D_0\mathrm{cosh}D_0l+F_0\mathrm{sinh}D_0l}},`$ (57a)
$`\alpha _\pm `$ $`=`$ $`{\displaystyle \frac{G_0}{D_0\mathrm{cosh}D_0l+F_0\mathrm{sinh}D_0l}}.`$ (57b)
The average of the Green’s function is defined as the integral over the phase space, which in the trajectory parameterization of the phase space is

$$\widehat{𝒢}(\epsilon ,𝐫,𝐧)=\gamma _{\mathrm{esc}}_0^{\mathrm{}}𝑑lP(l)_0^l\widehat{𝒢}_l(\epsilon ,x)𝑑\zeta ,$$
(58)
where we first integrate $`\widehat{𝒢}(\epsilon ,x)`$ along a trajectory of length $`l`$, and then average the result over trajectories of different length $`l`$ with the weight

$$P(l)=\gamma _{\mathrm{esc}}\mathrm{exp}(\gamma _{\mathrm{esc}}l).$$
(59)
Here $`\gamma _{\mathrm{esc}}`$ is escape rate, or inverse average escape time $`\tau _{\mathrm{esc}}`$, defined by Eq. (2).

Performing the integrations, we obtain:

$$\widehat{𝒢}(\epsilon ,𝐫,𝐧)=\frac{1}{\mathrm{cosh}\theta }\left(\begin{array}{cc}\eta _1\mathrm{sinh}\theta \eta _2& \eta _1+\eta _2\mathrm{sinh}\theta \\ \eta _1+\eta _2\mathrm{sinh}\theta & \eta _2\eta _1\mathrm{sinh}\theta \end{array}\right),$$
(60)
and coefficients $`\eta _{1,2}`$ are

$`\eta _1`$ $`=`$ $`1+2{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k\gamma _{\mathrm{esc}}^2}{(2F_0k\mathrm{cosh}\theta +\gamma _{\mathrm{esc}})^2}}\mathrm{tanh}^{2k}{\displaystyle \frac{\theta }{2}},`$
$`\eta _2`$ $`=`$ $`{\displaystyle \frac{\gamma _{\mathrm{esc}}}{F_0\mathrm{cosh}\theta }}\mathrm{tanh}{\displaystyle \frac{\theta }{2}}`$
$`+`$ $`{\displaystyle \frac{2\gamma _{\mathrm{esc}}}{F_0\mathrm{sinh}\theta }}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k\gamma _{\mathrm{esc}}}{2F_0k\mathrm{cosh}\theta +\gamma _{\mathrm{esc}}}}\mathrm{tanh}^{2k}{\displaystyle \frac{\theta }{2}}.`$
Here we introduced a shorthand notation $`\theta `$, defined through $`G_0=F_0\mathrm{sinh}\theta `$ and $`D_0=F_0\mathrm{cosh}\theta `$.

In these calculations, Eqs. (52b) and (60) represent a self consistent system of equation for $`F_0`$ and $`G_0`$ (or $`\theta `$). The solution of this self consistency equation gives the electron density of states in the grain:

$$𝒩(\epsilon )=\frac{2}{\delta _1}\mathrm{Re}\left\{\eta _1\mathrm{tanh}\theta \eta _2\frac{1}{\mathrm{cosh}\theta }\right\}.$$
(62)
The density of states, in the form of Eq. (62), supplied with Eqs. (52) and (60), are obtained here for the first time. These equations can be solved numerically for an arbitrary value of $`\tau _0\gamma _{\mathrm{esc}}`$. In Fig. 4 we plot $`𝒩(\epsilon )`$ for weak disorder, $`\tau _0\gamma _{\mathrm{esc}}=10`$, and for moderately strong disorder, $`\tau _0\gamma _{\mathrm{esc}}=0.7`$, by dash-dotted and dotted lines, respectively. For weak disorder, the density of states is similar to the semiclassical result \[dashed line\], although it exhibits the gap $`E_\mathrm{g}0.09\gamma _{\mathrm{esc}}`$. The gap is not noticeable in Fig. 4, since the overall number of states near the gap is already exponentially small. As the strength of disorder increases, $`\tau _0^1\gamma _{\mathrm{esc}}`$, the density of states acquire the form, similar to the random matrix result \[solid line\], see the dotted line in Fig. 4 for the $`\tau _0\gamma _{\mathrm{esc}}=0.7`$ case shown only for $`\epsilon <0.5\gamma _{\mathrm{esc}}`$. We also present the energy gap $`E_\mathrm{g}`$ as a function of the mean free time $`\tau _0`$. As $`\tau _0`$ increases, the gap size monotonically decreases, see Fig. 5.

Now we show, that both the semiclassical resultLN ; ScB  and the random matrix resultMBFB  may be obtained from Eqs. (52), (60) and (62).

In the absence of impurity scattering, $`\tau _0\mathrm{}`$, we obtain the ballistic semiclassical result. Indeed, we have $`G_0=D_0=i\epsilon `$, $`F_00`$ and $`\theta \mathrm{}`$. The density of states is determined by $`\eta _1`$:

$$𝒩(\epsilon )=\frac{2}{\delta _1}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\frac{(1)^k\gamma _{\mathrm{esc}}^2}{(2ik\epsilon +\gamma _{\mathrm{esc}})^2}.$$
(63)
This result is equivalent to the semiclassical expression for the density of states in the Andreev grain:

$$𝒩(\epsilon )=\frac{2}{\delta _1}\frac{\pi ^2\gamma _{\mathrm{esc}}^2}{4\epsilon ^2}\frac{\mathrm{cosh}\pi \gamma _{\mathrm{esc}}/2\epsilon }{\mathrm{sinh}^2\pi \gamma _{\mathrm{esc}}/2\epsilon },$$
(64)
which was derived in ref. \[ScB, \]. At large energy, the density of states in the clean semiclassical system approaches its value in the normal state according to

$$𝒩(\epsilon )=\frac{2}{\delta _1}\left(1+\frac{\pi ^2}{24}\frac{\gamma _{\mathrm{esc}}^2}{\epsilon ^2}\right).$$
(65)
In the opposite limit, $`\epsilon \tau _01`$ and $`\gamma _{\mathrm{esc}}\tau _01`$, we neglect series in Eq. (52b) for $`\eta _{1,2}`$, since the series are of the second order in $`\gamma _{\mathrm{esc}}\tau _0`$. To the zeroth order in $`\epsilon \tau _0`$ and $`\gamma _{\mathrm{esc}}\tau _0`$ the self-consistency equation, Eq. (52b), is satisfied by arbitrary $`\theta `$. Keeping the first order terms in $`\epsilon \tau _0`$ and $`\gamma _{\mathrm{esc}}\tau _0`$, we obtain the equation for $`\theta `$:

$$i\epsilon =\gamma _{\mathrm{esc}}\mathrm{cosh}\theta \mathrm{tanh}\frac{\theta }{2}$$
(66)
and the expression for the density of states

$$𝒩(\epsilon )=\frac{2}{\delta _1}\mathrm{Re}\left\{\mathrm{tanh}\theta \right\}.$$
(67)
Equations (66) and (67) are equivalent to Eqs. (46) and (47) after the substitution $`\theta =2i\psi `$. We conclude, that the density of states to the lowest order in $`\epsilon \tau _0`$ and $`\gamma _{\mathrm{esc}}\tau _0`$ corresponds to the random matrix theory result, see Eqs. (48) and (49) and ref. \[MBFB, \].

The most drastic difference between the density of states corresponding to different values of the mean free time $`\tau _0`$ appears at small energy, see Figs. 4 and 5. This difference was discussed in refs. \[ScB, ; LN, ; MBFB, ; TSA, ; AB, ; Bnew, \]. Here we emphasize that the difference exists even at high energy. Comparing Eqs. (49) and (65), we notice that already the lowest order in $`\gamma _{\mathrm{esc}}/\epsilon `$ term differs by a numerical factor of $`\pi ^2/12`$. In general, a high energy asymptote $`\epsilon \gamma _{\mathrm{esc}}`$ of the density of states can be written in the form:

$$𝒩(\epsilon )\frac{2}{\delta _1}\left(1+\phi (\epsilon \tau _0)\frac{\gamma _{\mathrm{esc}}^2}{2\epsilon ^2}\right),$$
(68a)
where $`\phi (\epsilon \tau _0)`$ is
$$\phi (y)=\mathrm{Re}\left\{\frac{1+2i2y\pi ^2y^2/12}{(1+iy)^2}\right\}.$$
(68b)
We observe, that transition from the dirty limit, $`\epsilon \tau _01`$, to the clean limit, $`\epsilon \tau _01`$, occurs at energy of the order of the scattering rate, $`\epsilon 1/\tau _0`$.

We discuss the difference between the high energy limits of the semiclassical and the random matrix theoryfn1  results in more details. Mathematically, the different coefficient in front of $`\epsilon ^2/\gamma _{\mathrm{esc}}^2`$ in Eqs. (49) and (65) may be attributed to $`(1)^k/k^2`$ over $`k`$, when the random matrix theory result is obtained only from the first term of the sum, while the semiclassical result corresponds to $`_{k=1}^{\mathrm{}}(1)^k/k^2=\pi ^2/12`$. We argue that this mathematical difference between Eqs. (49) and (65) has a physical meaning.

In the semiclassical limit, the energy levels in the grain are usually interpreted as standing waves along classical trajectories.LN ; ScB  Alternatively the contribution to the density of states may be explained as the effect of interference of repetitive motion of electron-hole pairs along a classical trajectory, see Eq. (63). A scattering off impurities switches electrons and holes between different trajectories and destroys the interference effect of such repetitive motion. Due to the scattering, an electron-hole pair travels along any trajectory only once, and the sum over $`k`$ in Eq. (63) – the number of repetitions – is limited to $`k=\pm 1`$.

## IV Quantum Chaos

In this Section we consider a clean metal grain with chaotic smooth potential and consider the effect of quantum diffraction on electron density of states.

### IV.1 Qualitative discussion

In quantum mechanics electron propagation is described by the evolution of a wave packet. The electron Green function, see Eq. (7), at some point $`𝐱`$ of the phase space is actually determined by averaging over a set of wave functions, which constitute a wave packet with center of mass at point $`𝐱`$. In this sense, the effect of quantum diffraction is similar to the effect of disorder. We emphasize, however, that the quantum diffraction mixes wave functions of electron-hole pairs at close points of the phase space within a minimal wave packet, see below. On the contrary, scattering off point impurities in disordered system mixes electron-hole wave functions propagating in completely different directions. Thus, the scattering off point impurities results in the mixing of distant points of the phase space.

First we estimate the size of a typical wave packet in the phase space and then discuss the effect of chaotic dynamics on the electron Green function.

We consider a wave packet, which consists of two electron-hole pairs, travelling from points $`C^{}`$ and $`A^{}`$ to points $`B^{}`$ and $`D^{}`$ in Fig. 6. The width of the wave packet is $`\delta \rho `$ and pairs’ momenta are misaligned within angle $`\delta \varphi `$, so that the momentum distribution of the wave packet has width $`p_{}=p_\mathrm{F}\delta \varphi `$, where $`p_\mathrm{F}=2\pi \mathrm{}/\lambda _\mathrm{F}`$ is the Fermi momentum and $`\lambda _\mathrm{F}`$ is the Fermi wavelength. According to the uncertainty principle, we have $`p_{}\delta \rho 2\pi \mathrm{}`$. Taking into account the relation $`\delta \rho L\delta \varphi `$, see Fig. 6, we obtain

$$\varphi _{\mathrm{min}}=\sqrt{\lambda _\mathrm{F}/L},\rho _{\mathrm{min}}=\sqrt{\lambda _\mathrm{F}L}$$
(69)
for the momentum deviation angle and for the coordinate width of a minimal wave packet, respectively.

The Eilenberger equation, Eq. (9), for a clean system does not describe quantum diffraction. The way to mimic quantum diffraction is to introduce small angle scattering potential, see refs. \[AL1, ; AAL, \]. This scattering introduces mixing, or “interaction”, between electron-hole wave functions, which approach each other within a small distance in the phase space, described by Eqs. (69). The small angle scattering suppresses the oscillating components of the Green functions and destroys the semiclassical description of electron-hole motion in the system.

We develop the application of the method of small angle potential to calculations of density of states in Andreev billiard in Appendix A, using the technique developed in ref. \[AL1, \]. Below in the text we utilize somewhat more intuitive approach. We represent the Green function in the Riccati parameterization, see Eq. (15) and rewrite it in terms of well defined series in $`a_\pm ^k(\epsilon ,𝐱)`$. To include the effect of quantum diffraction, we identify powers of $`a_\pm (\epsilon ,𝐱)`$ as weighted average over the minimal wave packet of the electron-hole pairs, see Fig. 6:

$$a_\pm ^k(\epsilon ,𝐱)a_\pm ^k(\epsilon ,𝐱)=\underset{i=1}{\overset{k}{}}P(𝐱𝐱_i)a_\pm (\epsilon ,𝐱_i)𝑑𝐱_i.$$
(70)
Function $`P(𝐱𝐱^{})`$ represents the coordinate and momentum distribution within a minimal wave packet and vanishes fast on coordinate scale of $`\rho _{\mathrm{min}}`$ or for momentum deviation on angle larger than $`\delta \varphi _{\mathrm{min}}`$, see Eqs. (69). A detailed form of the function $`P(𝐱𝐱^{})`$ is not important, since, as we show below, the density of states depends logarithmically on the width of the function $`P(𝐱𝐱^{})`$.

In chaotic system the distance between two trajectories grows exponentially in time $`t`$: $`\rho (t)=\rho _{\mathrm{min}}\mathrm{exp}(\lambda t)`$ and $`\varphi (t)=\varphi _{\mathrm{min}}\mathrm{exp}(\lambda t)`$, where $`\lambda >0`$ is the Lyapunov exponent of the classical chaotic motion. Consequently, even though the size of the wave packet is initially small, time evolution leads to the spreading of the wave packet over the phase space. Generally, the spreading of the wave packet is characterized by the Ehrenfest time $`t_\mathrm{E}`$. The Ehrenfest time shows how long it takes for a wave packet to spread over the whole phase space of the system and is equal to

$$t_\mathrm{E}=\frac{1}{\lambda }\mathrm{ln}\frac{L}{\lambda _\mathrm{F}}=\frac{2}{\lambda }\mathrm{ln}\frac{L}{\rho _{\mathrm{min}}}.$$
(71)
We argue, however, that the relevant time scale for the problem of proximity effect in Andreev billiards depends on the size of the contact $`b`$ and is given by

$$\tau _\mathrm{E}=\frac{1}{\lambda }\mathrm{ln}\frac{b}{\rho _{\mathrm{min}}}=\frac{t_\mathrm{E}}{2}\frac{1}{\lambda }\mathrm{ln}\frac{L}{b}.$$
(72)
The second term in the right hand side of Eq. (72) may be omitted, since it contains a logarithm of macroscopic parameters $`b`$ and $`L`$. Nevertheless we explain the appearance of the scale $`b`$ as follows. If electron propagates from points $`B^{}`$ or $`D^{}`$ along trajectories $`B^{}B`$ or $`D^{}D`$, see Fig. 6a, for time smaller than $`\tau _\mathrm{E}`$, the spreading of the initial wave packet is smaller than the size $`b`$ of the SN contact and the two trajectories ends nearly simultaneously at points $`B`$ and $`D`$. We conclude that the contribution of electron-hole wave functions to the electron Green function originates from coherent motion of electron-hole pairs along these two trajectories. This coherent contribution may be viewed as a repetitive motion along the same trajectory, see the text at the end of Sec. III.2.2.

On the other hand, if electron propagates along $`B^{}B`$ and $`D^{}D`$ for time longer than $`\tau _\mathrm{E}`$, the wave packet spreads on the scale larger than the size $`b`$ of the SN contact. As shown in Fig. 6b, although the electron-hole pair travelling along $`D^{}D`$ experiences Andreev reflection at point $`D`$, the other pair travels along $`B^{}B`$ and misses the SN contact. It continues its motion further from point $`B`$ and explores the phase space, different from the one reached by the $`D^{}D`$ trajectory. Thus, even though the wave packet originates from the contribution of electron-hole wave functions at close points in the phase space, the corresponding wave functions are quite different, since they represent motion in distant parts of the phase space. The averaging over the wave packet leads to appearance of a component of the Green function, produced from different parts of the system phase space, known as a “zero mode”.

We expect that if the Ehrenfest time $`\tau _\mathrm{E}`$ is much longer than the time of electron propagation along a typical trajectory, then the electron Green function may be determined by semiclassical methods, see refs. \[LN, ; ScB, \], while in the opposite limit of short $`\tau _\mathrm{E}`$, the quantum diffraction destroys the repetitive modes and leads to strong coupling between electron-hole wave functions from distant parts of the phase space. In this sense, motion of electron-hole pairs in systems with short $`\tau _\mathrm{E}\tau _{\mathrm{esc}}`$ actually coincides with motion in disordered systems, and the random matrix resultMBFB  describes the density of states, see Sec. III.2.1.

### IV.2 Quantitative analysis

In this Subsection we perform an analytical analysis of the density of states in clean chaotic Andreev billiards. First we consider the high energy limit $`\epsilon \gamma _{\mathrm{esc}}`$, when the perturbation theory is applicable. In the perturbation theory we identify the relevant features of the Eilenberger equation for a description of motion of electron-hole pairs in the regime of quantum chaos. Next, we construct a self-consistency equation, which determines the density of states at arbitrary energy $`\epsilon \gamma _{\mathrm{esc}}`$.

#### IV.2.1 High energy limit

As we showed in section III.2.2, the density of states of disordered systems approaches the normal density of states as $`\gamma _{\mathrm{esc}}^2/\epsilon ^2`$ at large energy $`\epsilon \gamma _{\mathrm{esc}}`$, but the coefficient of this asymptote depends on the strength of disorder, characterized by the mean free time $`\tau _0`$, see Eq. (68b). A similar asymptotic behavior exists in quantum chaotic systems and now the coefficient depends on the value of the Ehrenfest time $`\tau _\mathrm{E}`$. To prove this statement, we calculate the density of states at high energy ($`\epsilon \gamma _{\mathrm{esc}}`$).

We use the Riccati parameterization of the Green function and apply Eq. (70) for calculation of powers $`a_\pm (\epsilon ,𝐱)`$ at coincident points, which appear in the expression for the density of states, see Eq. (17). A more formal derivation with Perron-Frobenius differential operator is presented in Appendix A.

The density of states is determined by the sum of $`\overline{a_+^k(\epsilon ,𝐱)a_{}^k(\epsilon ,𝐱)}`$ over all positive integer $`k`$. Here $`\overline{(\mathrm{})}`$ stands for the averaging over the phase space of electrons in the grain, while $`\mathrm{}`$ denotes averaging within a minimal wave packet according to Eq. (70). We observe that $`a_+^k(\epsilon ,𝐱)`$ and $`a_{}^k(\epsilon ,𝐱)`$ are mutually independent, since the corresponding electron-hole pairs move in completely different parts of the phase space, approaching point $`𝐱`$ from two opposite directions (notice the different signs before the Liouville operators in Eqs. (16)). Due to independence of $`a_+^k(\epsilon ,𝐱)`$ and $`a_{}^k(\epsilon ,𝐱)`$, we write

$$\overline{a_+^k(\epsilon ,𝐱)a_{}^k(\epsilon ,𝐱)}=\overline{a_+^k(\epsilon ,𝐱)}\overline{a_{}^k(\epsilon ,𝐱)},$$
(73)
so that the density of states in the grain is determined by the product of $`\overline{a_+^k(\epsilon ,𝐱)}`$ and $`\overline{a_{}^k(\epsilon ,𝐱)}`$.

We start with calculation of the averages of electron-hole wave functions $`\overline{a_\pm (\epsilon ,𝐱)}`$. We integrate both sides of Eq. (16) over the phase space of electrons in the grain $`𝒱_\mathrm{g}\mathrm{\Omega }_d`$. As we already mentioned, the scattering term vanishes after integration, see Eq. (18). The integral of $`_{\mathrm{cm}}a_\pm (\epsilon ,𝐱)`$ reduces to the values of $`a_\pm (\epsilon ,𝐱)`$ at the SN contact. At $`\epsilon \gamma _{\mathrm{esc}}`$, only the contribution from the trajectory terminal with restrained value of $`a_\pm (\epsilon ,\zeta )`$ is important, since the unrestrained value of $`a_\pm (\epsilon ,\zeta )`$ at the SN contact oscillates fast and thus averages out. We obtain

$$\overline{a_\pm (\epsilon ,𝐱)}\frac{\gamma _{\mathrm{esc}}}{2i\epsilon }.$$
(74)
Next, we calculate the contribution to the density of states from the $`k=2`$ term in Eq. (17), $`\overline{a_+^2(\epsilon ,𝐱)}`$, to the second order in $`\gamma _{\mathrm{esc}}/\epsilon `$. We show that this term behaves differently in semiclassical and quantum mechanics. Particularly, the quantum diffraction destroys the semiclassical contribution to the density of states, originating from the repetitive motion.

A similar problem was solved in ref. \[AL1, \], where the weak localization correction to the conductance of quantum chaotic systems was studied. It was shown that the weak localization correction originates as a result of quantum diffraction which couples Cooperon and diffuson modes. Following the technique, developed in refs. \[AL1, ; AAL, \], we introduce new functions $`A_2^\pm (\epsilon ,𝐱_1,𝐱_2)`$ with arbitrary separated points $`𝐱_1=(𝐫_1,𝐧_1)`$ and $`𝐱_2=(𝐫_2,𝐧_2)`$:

$$A_2^\pm (\epsilon ,𝐱_1,𝐱_2)=a_\pm (\epsilon ,𝐱_1)a_\pm (\epsilon ,𝐱_2).$$
(75)
We are interested in the limit of small relative angle $`\varphi `$ between vectors $`𝐧_1`$ and $`𝐧_2`$ within a minimal wave packet, Eq. (69). We change variables from $`𝐱_1`$ and $`𝐱_2`$ to the coordinates of the center of mass $`𝐑=(𝐫_1+𝐫_2)/2,𝐍(𝐧_1+𝐧_2)/2`$ (for $`|\varphi |1`$) and to the relative coordinates $`𝝆=𝐫_1𝐫_2`$ and $`\varphi `$. Under the condition $`|\varphi |1`$, the vector of relative distance $`𝝆`$ is perpendicular to the direction of $`𝐍`$ and we refer to it through its projection $`\rho `$ on the axis perpendicular to $`𝐍`$. At $`|\varphi |1`$ the equation for $`A_2^\pm (\epsilon ,𝐱_1,𝐱_2)`$ in new variables is

$`\left(4i\epsilon _{\mathrm{cm}}(\rho ,\varphi )\right)A_2^\pm (\epsilon ,𝐑,𝐍,\rho ,\varphi )=0,`$ (76)
where

$$_{\mathrm{cm}}=v_\mathrm{F}𝐍\frac{}{𝐑}\frac{1}{p_\mathrm{F}}\frac{U_\mathrm{s}}{𝐑}\frac{}{𝐍}$$
(77a)
is the part of the Liouville operator describing the motion of the center of mass, and
$$(\rho ,\varphi )=v_\mathrm{F}\varphi \frac{}{\rho }\frac{1}{p_\mathrm{F}}\frac{^2U_\mathrm{s}}{𝐑_{}^2}\rho \frac{}{\varphi }$$
(77b)
is the Liouville operator of relative motion of two electron-hole pairs.

We average $`A_2^\pm (\epsilon ,𝐑,𝐍,\rho ,\varphi )`$ over the position of center of mass $`𝐍`$ and $`𝐑`$ in the phase space. Integration of the second term with the Liouville operator $`_{\mathrm{cm}}`$ by parts gives values of $`A_2^\pm `$ at the boundaries. Particularly, due to the boundary conditions, see Table 1, for outgoing electrons $`a_\pm (\epsilon ,\zeta =l)=1`$ and we have $`A_2^+(\epsilon ,𝐑,𝐍,\rho ,\varphi )=1`$ at the SN contact, provided $`\mathrm{𝐍𝐧}_\mathrm{s}<0`$ ($`𝐧_\mathrm{s}`$ is a normal vector to the SN inerface). In the opposite limit $`\mathrm{𝐍𝐧}_\mathrm{s}>0`$, the value of $`A_2^+(\epsilon ,𝐑,𝐍,\rho ,\varphi )`$ is unrestrained. Similar analysis provides values of $`A_2^{}(\epsilon ,𝐑,𝐍,\rho ,\varphi )`$ at the SN contact. We introduce notation $`s_2=A_2^\pm _{\mathrm{sn}}/\overline{A_2^\pm }`$ for the ratio of the average unrestricted value $`A_2^\pm `$ at the SN contact, $`A_2^\pm _{\mathrm{sn}}`$, to the average of $`A_2^\pm `$ over the whole phase space. We obtain

$$\left(4i\epsilon +\frac{s_2}{\tau _{\mathrm{esc}}}(\rho ,\varphi )\right)A_2^\pm (\epsilon ,\rho ,\varphi )=\frac{1}{\tau _{\mathrm{esc}}}.$$
(78)
At high energy $`s_21`$.

The Liouville operator $`(\rho ,\varphi )`$ in Eq. (78) can be rewritten in new variables $`z`$ and $`\chi `$, introduced as

$$z=\frac{1}{2}\mathrm{ln}\frac{\varphi ^2L^2+\rho ^2}{L^2},\chi =\mathrm{arctan}\frac{\varphi L}{\rho }.$$
(79)
As was shown in ref. \[AL1, \], the variable $`\chi `$ is irrelevant for quantity $`A_2^\pm (\epsilon ,\rho ,\varphi )`$, described by a solution of Eq. (78). Indeed, we represent evolution of the wave packet in time $`t`$, which is a conjugate variable to energy variable $`\epsilon `$ in the left hand side of Eq. (78), see Fig. 7. The characteristic scale of the wave packet, represented by the variable $`\mathrm{exp}(z)`$ grows in time as $`\mathrm{exp}(\lambda t)`$ with the Lyapunov exponent $`\lambda `$. A particular configuration of the wave packet is described by variable $`\chi `$, which significantly fluctuates as the wave packet explores the phase space and gets averaged out after integration of $`A_2^\pm `$ over a position of the center of mass. \[Actually this averaging occurs at ergodic time scale, which is much shorter than any other time scale in the problem.\] The fluctuations of the Lyapunov exponent $`\lambda `$ may also be taken into account. The resulting form of Eq. (78) isAL1 :

$$\left(4i\epsilon +\gamma _{\mathrm{esc}}\lambda \frac{}{z}\frac{\lambda _2}{2}\frac{^2}{z^2}\right)A_2^\pm (\epsilon ,z)=\gamma _{\mathrm{esc}},$$
(80)
where $`\lambda `$ is the average value of the Lyapunov exponent over the electron phase space in the grain and $`\lambda _2`$ represents the deviation of Lyapunov exponent for different points of the phase space, see ref. \[AL1, \] for more details.

As the size of the wave packet reaches the size of the SN contact $`b`$, and $`z\mathrm{ln}b/L`$, the $`a_\pm (\epsilon ,𝐱)`$ functions in $`\overline{A_2^\pm (\epsilon ,𝐱_1,𝐱_2)}`$ become independent. Indeed, these functions are solutions of the Eilenberger equation along different trajectories of significantly varying lengths. We may write the following boundary condition:

$$A_2^\pm (\epsilon ,z=\mathrm{ln}b/L)=\overline{a_\pm (\epsilon ,𝐱)}^2,$$
(81)
where $`\overline{a_\pm (\epsilon ,𝐱)}`$ is the average of $`a_\pm (\epsilon ,𝐱)`$ over the phase space and is given by Eq. (74). \[We should distinguish the average $`a_\pm (\epsilon ,𝐱)`$ and the zero mode components of the Green function. More precisely, the right hand side of Eq. (81) contains the zero mode value of $`a_\pm (\epsilon ,𝐱)`$, as we will discuss in the next subsection. At high energy $`\epsilon \gamma _{\mathrm{esc}}`$ both quantities are equal.\]

To solve Eqs. (78) and (81), we represent $`A_2^\pm (\epsilon ,\rho ,\varphi )`$ as a sum of two functions:

$$A_2^\pm (\epsilon ,\rho ,\varphi )=A_{2,\mathrm{b}}^\pm (\epsilon )+A_{2,\mathrm{q}}^\pm (\epsilon ,z),$$
(82)
where

$$A_{2,\mathrm{b}}^\pm (\epsilon )=\frac{\gamma _{\mathrm{esc}}}{4i\epsilon +\gamma _{\mathrm{esc}}}\frac{\gamma _{\mathrm{esc}}}{4i\epsilon }$$
(83)
is a solution of Eq. (80) without any boundary conditions imposed and $`\overline{A_{2,q}^\pm (\epsilon ,z)}`$ is a solution of homogeneous equation Eq. (80) with certain boundary conditions at $`z=\mathrm{ln}b/L`$, so that $`A_2^\pm (\epsilon ,z)`$ in the form of Eq. (82) satisfies Eq. (81).

A general solution of the homogeneous Eq. (80) was found in ref. \[AL1, \] in the limit $`\gamma _{\mathrm{esc}},\epsilon \lambda `$:

$`w(\omega ,z)=\mathrm{exp}\left[\left({\displaystyle \frac{i\omega }{\lambda }}+{\displaystyle \frac{\lambda _2\omega ^2}{2\lambda ^3}}\right)z\right],`$ (84)
with $`\omega =4\epsilon i\gamma _{\mathrm{esc}}`$. Taking into account the boundary conditions, we find

$`A_{2,\mathrm{q}}^\pm (\epsilon ,z)`$ $`=`$ $`\left(\overline{a_\pm (\epsilon ,𝐱)}^2\overline{A_{2,b}^\pm (\epsilon )}\right)`$ (85)
$`\times `$ $`w(4\epsilon i\gamma _{\mathrm{esc}},z\mathrm{ln}b/L).`$
This equation describes the Lyapunov divergence of classical trajectories, once a finite displacement between two trajectories was created by quantum diffraction. As was discussed in ref. \[AL1, \], Eq. (84) neglects the effect of quantum diffraction and a more accurate solution for the function $`w(\omega ,z)`$ is needed.AL1  Indeed, according to Eq. (17), the density of states contains term $`A_2^\pm (\epsilon ,\rho ,\varphi )`$, taken at coincident points, i.e. $`\rho ,\varphi 0`$. This limit corresponds to $`z\mathrm{}`$, when $`w(\epsilon ,z\mathrm{})=0`$. The latter equation means that the effect of quantum diffraction, which couples Green functions at different semiclassical trajectories is not described by the Eilenberger equation without scattering term ($`1/\tau _\mathrm{q}=0`$). The scattering term modifies the equation for function $`w(\epsilon ,z)`$, see Appendix A for details. We argue that instead of working in the limit $`z\mathrm{}`$, within a logarithmic accuracy we can average quantity $`A_2^\pm (\epsilon ,\rho ,\varphi )`$ over $`\rho `$ and $`\varphi `$ within a minimal wave packet, Eq. (69). We write

$`\mathrm{\Gamma }_2(\epsilon )`$ $`=`$ $`w(4\epsilon i\gamma _{\mathrm{esc}},z\mathrm{ln}{\displaystyle \frac{b}{L}}).`$ (86)
Performing averaging over the minimal wave packet

$`w(\omega ,z)`$ $``$ $`{\displaystyle w(\omega ,z)e^{(\varphi ^2L^2+\rho ^2)/2\rho _{\mathrm{min}}^2}𝑑\varphi 𝑑\rho }`$
with $`z=\mathrm{ln}\sqrt{(\varphi ^2L^2+\rho ^2)/b^2}`$ we obtain

$$\mathrm{\Gamma }_2(\epsilon )=\mathrm{exp}\left(\gamma _{\mathrm{esc}}\tau _\mathrm{E}4i\epsilon \tau _\mathrm{E}\frac{8\lambda _2\epsilon ^2}{\lambda ^2}\tau _\mathrm{E}\right)$$
(87)
and $`\tau _\mathrm{E}=(1/\lambda )\mathrm{ln}b/\rho _{\mathrm{min}}`$ is the Ehrenfest time, defined by Eq. (72). Combining Eqs. (74), (83) and (85) we find

$$A_2^\pm (\epsilon ,z)=\frac{\gamma _{\mathrm{esc}}}{4i\epsilon }\left\{1\mathrm{\Gamma }_2(\epsilon )\right\}$$
(88)
to the lowest order in $`\gamma _{\mathrm{esc}}/\epsilon `$.

Equation (88) has a simple interpretation. Function $`\mathrm{\Gamma }_2(\epsilon )`$ is the probability amplitude, that two electron-hole pairs, initially situated within the same minimal wave packet at point $`𝐱`$ of the phase space, escape the grain independently. In other words, $`\mathrm{\Gamma }_2(\epsilon )`$ is the probability amplitude that a minimal wave packet, constructed at point $`𝐱`$ expand to size $`b`$ of the contact before it leaves the grain. Consequently, $`1\mathrm{\Gamma }_2(\epsilon )`$ is the probability amplitude for the wave packet to leave the grain before it decays due to quantum diffraction.

We generalize the above calculations of $`\overline{A_{k=2}^\pm }`$ for arbitrary integer $`k>1`$. We defineGM

$$A_k^\pm (\epsilon ,\{𝐱_i\})=\underset{i=1}{\overset{k}{}}a_\pm (\epsilon ,𝐱_i)$$
(89)
and write the equation for $`\overline{A_k^\pm (\epsilon ,\{𝐱_i\})}`$, averaged over the center of mass of variables $`\{x_i\}`$:

$$\left(2ki\epsilon +\frac{1}{\tau _{\mathrm{esc}}}(𝝆,\mathit{\varphi })\right)A_k^\pm (\epsilon ,𝝆,\mathit{\varphi })=\frac{1}{\tau _{\mathrm{esc}}},$$
(90)
which is a generalization of Eq. (78) to arbitrary $`k`$. Here the variables $`𝝆=\{\rho _i\}`$ and $`\mathit{\varphi }=\{\varphi _i\}`$ constitute a set of $`2(k1)`$ variables. We argue that $`2k3`$ of these variables determine the relative position of $`k`$ electron-hole pairs within the corresponding wave packet and thus are irrelevant, while one variable determines the scale of the wave packet in the phase space $`z=\mathrm{ln}(L^2\varphi _i^2+\rho _i^2)/2`$, see Fig. 7. When the wave packet reaches the size of the SN contact, the wave functions of electron-hole pairs $`a_\pm (\epsilon ,𝐱_i)`$ become independent and we have the boundary condition $`z=\mathrm{ln}b/L`$:

$$A_k^\pm (\epsilon ,z=\mathrm{ln}b/L)=\overline{a_\pm (\epsilon ,𝐱)}^k\frac{\gamma _{\mathrm{esc}}^k}{(i\epsilon )^k}.$$
(91)
Following further along the steps of calculations of $`A_2^\pm (\epsilon ,z=\mathrm{ln}b/L)`$, we obtain to the lowest order in $`\gamma _{\mathrm{esc}}/\epsilon `$

$$A_k^\pm (\epsilon ,z)=\frac{\gamma _{\mathrm{esc}}}{2ik\epsilon }\left\{1\mathrm{\Gamma }_k(\epsilon )\right\}$$
(92)
with

$$\mathrm{\Gamma }_k(\epsilon )=\mathrm{exp}\left(\gamma _{\mathrm{esc}}\tau _\mathrm{E}2ik\epsilon \tau _\mathrm{E}\frac{2k^2\lambda _2\epsilon ^2}{\lambda ^2}\tau _\mathrm{E}\right).$$
(93)
Using Eq. (17) we obtain the density of states at large energy, $`\epsilon \gamma _{\mathrm{esc}}`$ in the form

$$𝒩(\epsilon )=\frac{2}{\delta _1}\left[1+\frac{\gamma _{\mathrm{esc}}^2}{2\epsilon ^2}\left(1\underset{k=2}{\overset{\mathrm{}}{}}\frac{(1)^k}{k^2}\mathrm{Re}\left\{1\mathrm{exp}\left(\left(2ik\epsilon +\gamma _{\mathrm{esc}}+\frac{2k^2\lambda _2\epsilon ^2}{\lambda ^2}\right)\tau _\mathrm{E}\right)\right\}^2\right)\right].$$
(94)
We found that the averages of functions $`\overline{a_+^k(\epsilon ,𝐫,𝐧)}`$ are expressed in terms of $`\mathrm{\Gamma }_k(\epsilon )`$. Particularly, function $`\mathrm{\Gamma }_2(\epsilon )`$ also appears in the expressions for *i)* the weak localization correction to the conductivityAL1 , *ii)* level statistics correlation functionAL2 , *iii)* noise of a current through a chaotic cavity.AAL

Functions $`\mathrm{\Gamma }_k(\epsilon )`$ oscillate with a period proportional to $`\tau _\mathrm{E}^1`$. Similar oscillations were found in the weak localization correction to the conductivityAL1  and in the energy level statistics.AL2  The origin of these oscillations has the following explanation. The effect of quantum diffraction can be neglected for motion at time shorter than the Ehrenfest time. Nevertheless the diffraction effects appear nearly instantaneously, with characteristic turn-on time equal to $`\lambda ^1`$, as time of motion approaches the Ehrenfest time. Thus we can separate two types of contributions to the Green function. One contribution represents semiclassical motion of electron-hole pairs and is responsible for oscillations in the density of states in the phase space. The other contribution is a constant in phase space component of the Green function, which corresponds to the zero mode.

To illustrate the oscillations of the density of states, we plot the tail of the density of states $`𝒩(\epsilon )`$ at large energies $`\epsilon \gamma _{\mathrm{esc}}`$ for $`\tau _\mathrm{E}=1.5\tau _{\mathrm{esc}}`$, shown by a solid line in Fig. 8. For comparison, we also plot the density of states for the random matrix limit (dashed line) and semiclassical limit (dash-dotted line).

We emphasize that it is the sharp turn-on of the effects of quantum diffraction at the Ehrenfest time leads to oscillations of the density of states. Fluctuations of the Lyapunov exponent suppress these oscillations. Indeed, the Lyapunov exponent varies for different points of the phase space and its fluctuations affect functions $`\mathrm{\Gamma }_k(\epsilon )`$ through the term $`\lambda _2\epsilon ^2\tau _\mathrm{E}/\lambda ^2`$. Thus for small energies, $`\epsilon \sqrt{\lambda ^2/\lambda _2\tau _\mathrm{E}}`$, we can disregard fluctuations of the Ehrenfest time. As energy increases, the oscillations of the density of states become suppressed, see Eq. (94). We also notice that in the quantum disorder, the process of scattering off impurities is Poissonian, and fluctuations of the free path of the order of the mean free path itself. Consequently, oscillations of the density of states in disordered systems are absent, see Eqs. (68).

#### IV.2.2 Arbitrary energy

Although the high energy expansion cannot be applied to the calculation of the density of states at energy comparable with the escape rate $`\gamma _{\mathrm{esc}}`$, exactly the same mechanism of quantum diffraction is responsible for the appearance of the gap in the density of states.

Equation for $`A_k^\pm (\rho ,\varphi )`$ becomes significantly simplified after averaging over the phase space of the grain. Nevertheless, this averaging happens on time scale of the order of ergodic time $`\tau _{\mathrm{erg}}`$. If we are interested in the properties of the spectrum at energy much smaller than the Thouless energy $`E_{\mathrm{Th}}`$ (the latter is equal to the inverse ergodic time, $`E_{\mathrm{Th}}=\tau _{\mathrm{erg}}^1`$), we can assume that the motion of electron-hole pairs constituting $`A_k^\pm (\rho ,\varphi )`$, is described by equation, averaged over irrelevant coordinates:

$$\left(2i\epsilon k\frac{}{\zeta }\lambda \frac{}{z}\frac{\lambda _2}{2}\frac{^2}{z^2}\right)A_k^+(\epsilon ,\zeta ,z)=0.$$
(95)
Here $`\zeta `$ stands for a coordinate of the center of mass of a wave packet. As we already discussed, the term $`^2/z^2`$ represents fluctuations of the Laypunov exponent and is important at large energy $`\epsilon \sqrt{\lambda /\tau _\mathrm{E}}`$. Below we neglect this term.

The function $`A_k^+(\epsilon ,\zeta ,z)`$ is defined at $`0<\zeta <l`$, where $`l`$ is the length of a trajectory, and has the boundary condition $`A_k^+(\epsilon ,l,z)=1`$ at $`\zeta =l`$, see Table 1. The other boundary condition exists at $`z=\mathrm{ln}b/L`$, which represents the statistical independence of electron-hole wave functions, separated by a distance $`b`$. We say about these functions that they correspond to a zero mode $`\overline{a}`$. Thus the boundary condition is $`A_k^\pm (\epsilon ,\zeta ,z=\mathrm{ln}b/L)=\overline{a}^k`$, see Fig. 9.

We write the solution to Eq. (95) as

$`A_k^+(\epsilon ,\zeta ,z)`$ $`=`$ $`\overline{a}^ke^{2ik\epsilon (z\mathrm{ln}b/L)/\lambda }\theta \left({\displaystyle \frac{z\mathrm{ln}b/L}{\lambda }}\zeta +l\right)`$ (96a)
$`+`$ $`e^{2ik\epsilon (\zeta l)}\theta \left(\zeta l{\displaystyle \frac{z\mathrm{ln}b/L}{\lambda }}\right).`$
Repeating the above analysis for $`A_k^{}(\epsilon ,x,z)`$, we obtain
$`A_k^{}(\epsilon ,\zeta ,z)`$ $`=`$ $`\overline{a}^ke^{2ik\epsilon (z\mathrm{ln}b/L)/\lambda }\theta \left({\displaystyle \frac{z\mathrm{ln}b/L}{\lambda }}\zeta \right)`$ (96b)
$`+`$ $`e^{2ik\epsilon \zeta }\theta \left(\zeta {\displaystyle \frac{z\mathrm{ln}b/L}{\lambda }}\right).`$
The density of states is described by a solution of Eilenberger equation, represented in the form of series in $`a_\pm ^k(\epsilon ,𝐱)`$. According to Eq. (70), the quantum diffraction leads to a finite separation between points $`𝐱_i`$ in the expression for $`a_\pm ^k(\epsilon ,𝐱)`$. The separation is determined by $`\rho _{\mathrm{min}}`$ and $`\varphi _{\mathrm{min}}`$ defined by Eq. (69). We conclude, that the Green function may be represented as a sum of functions $`A_k^\pm (\epsilon ,\zeta ,\mathrm{ln}\rho _{\mathrm{min}}/L)`$ over all positive integers $`k`$, see Eqs (70) and (89).

The functions

$`a_+(\epsilon ,\zeta )`$ $`=`$ $`\{\begin{array}{cc}\mathrm{exp}(2i\epsilon (\zeta l)),\hfill & l\tau _\mathrm{E}<\zeta <l\text{;}\hfill \\ a_0,\hfill & \text{ }0<x<l\tau _\mathrm{E}\text{,}\hfill \end{array}`$ (97c)
$`a_{}(\epsilon ,\zeta )`$ $`=`$ $`\{\begin{array}{cc}\mathrm{exp}(2i\epsilon \zeta ),\hfill & 0<\zeta <\tau _\mathrm{E}\text{;}\hfill \\ a_0,\hfill & \tau _\mathrm{E}<\zeta <l\text{.}\hfill \end{array}`$ (97f)
allow us to rewrite the sum for $`\widehat{𝒢}(\epsilon ,\zeta )`$ in terms of $`A_k^\pm (\epsilon ,\zeta ,z)`$ in a more compact form of Eq. (15). We introduced a new notation $`a_0`$ for the zero mode, related to $`\overline{a}`$ by the following equation $`a_0=e^{2i\epsilon \tau _\mathrm{E}}\overline{a}`$. The functions $`a_\pm (\epsilon ,𝐱)`$ have all properties we discussed above. They satisfy the boundary conditions at one terminal of a classical trajectory see Table 1, and coincide with semiclassical solutions of the Eilenberger equation at small distance $`\zeta \tau _\mathrm{E}`$ or $`l\zeta \tau _\mathrm{E}`$from this terminal, when the quantum diffraction effects may be disregarded. As an electron travels for time interval equal to the Ehrenfest time, $`\tau _\mathrm{E}`$, the quantum scattering suppresses the oscillations and $`a_\pm (\epsilon ,\zeta )`$ jumps to a zero mode value $`a_0`$.

Using Eqs. (97) to calculate the density of states at large energy, $`\epsilon \gamma _{\mathrm{esc}}`$, we obtain the result of Eq. (94). Indeed, by direct averaging of $`a_\pm (\epsilon ,𝐱)`$ in the form given by Eqs. (97), and comparing the result with Eq. (74), we find $`a_0=e^{2i\epsilon \tau _\mathrm{E}}\gamma _{\mathrm{esc}}/(2i\epsilon )`$, and then again integrating $`a_\pm ^k(\epsilon ,𝐱)`$ over the phase space we reproduce Eq. (92), apart of the term $`\lambda _2\epsilon ^2\tau _\mathrm{E}/\lambda ^2`$, which is important only at the extremely large energy $`\epsilon \sqrt{\lambda /\tau _\mathrm{E}}`$.

Below we perform the self-consistent calculations of the Green function in the Riccati parameterization, Eq. (15), with functions $`a_\pm (\epsilon ,𝐱)`$ given by Eqs. (97). Within the anzats of Eqs. (97), there exists only one parameter $`a_0`$, which completely determines the Green function in the grain. This single parameter can be found with the help of Eq. (18). Indeed, we substitute $`g(\epsilon ,𝐱)`$ and $`f(\epsilon ,𝐱)`$ in terms of the Riccati variables $`a_\pm (\epsilon ,𝐱)`$, write $`a_0=\mathrm{exp}(2i\psi )`$ and obtain the following equation for $`\psi `$:

$`\mathrm{tan}\psi `$ $`=`$ $`{\displaystyle \frac{\epsilon }{\gamma _{\mathrm{esc}}}}\left({\displaystyle \frac{e^{\gamma _{\mathrm{esc}}\tau _\mathrm{E}}}{\mathrm{cos}2\psi }}+2(1e^{\gamma _{\mathrm{esc}}\tau _\mathrm{E}})\mathrm{cos}2\psi \right)+2e^{\gamma _{\mathrm{esc}}\tau _\mathrm{E}}\mathrm{sin}2\psi \mathrm{ln}{\displaystyle \frac{\mathrm{cos}\psi }{\mathrm{cos}(\psi \epsilon \tau _\mathrm{E})}}`$ (98)
$`+`$ $`{\displaystyle _{\tau _\mathrm{E}}^{2\tau _\mathrm{E}}}\gamma _{\mathrm{esc}}e^{\gamma _{\mathrm{esc}}(l\tau _\mathrm{E})}\left\{{\displaystyle \frac{\mathrm{sin}\epsilon (2\tau _\mathrm{E}l)}{\mathrm{cos}\epsilon l}}+2\mathrm{sin}2\psi \mathrm{ln}{\displaystyle \frac{\mathrm{cos}\psi }{\mathrm{cos}(\psi +\epsilon (l\tau _\mathrm{E}))}}\right\}𝑑l`$
Representing $`\widehat{𝒢}(\epsilon ,𝐱)`$ in terms of functions $`a_\pm (\epsilon ,\zeta )`$, we find the density of states from Eq. (8):

$`𝒩(\epsilon )`$ $`=`$ $`{\displaystyle \frac{2}{\delta _1}}[{\displaystyle \frac{\pi ^2\gamma _{\mathrm{esc}}^2}{\epsilon ^2}}({\displaystyle \underset{n=0}{\overset{N_1^{}}{}}}(n+1/2)e^{\pi \gamma _{\mathrm{esc}}(n+1/2)/\epsilon }+{\displaystyle \underset{n=N_1^{}+1}{\overset{N_2^{}}{}}}(2\tau _\mathrm{E}\epsilon /\pi +(n+1/2))e^{\pi \gamma _{\mathrm{esc}}(n+1/2)/\epsilon })`$ (99)
$`+`$ $`e^{2\gamma _{\mathrm{esc}}\tau _\mathrm{E}}\mathrm{Im}\{\mathrm{tan}2\psi +\mathrm{ln}{\displaystyle \frac{\mathrm{cos}(\psi \epsilon \tau _\mathrm{E})}{\mathrm{cos}(\psi +\epsilon \tau _\mathrm{E})}}\}+{\displaystyle _{\tau _\mathrm{E}}^{2\tau _\mathrm{E}}}\gamma _{\mathrm{esc}}e^{\gamma _{\mathrm{esc}}l}\mathrm{Im}\left\{\mathrm{ln}{\displaystyle \frac{\mathrm{cos}(\psi \epsilon (l\tau _\mathrm{E}))}{\mathrm{cos}(\psi +\epsilon (l\tau _\mathrm{E}))}}\right\}dl].`$
Here $`N_j^{}=\mathrm{int}(j\epsilon \tau _\mathrm{E}/\pi 1/2)`$ is the integer part of $`j\epsilon \tau _\mathrm{E}/\pi 1/2`$. The terms in the first line of Eq. (99) originate from purely ballistic contributions of the electron Green function, while the term, containing $`\mathrm{tanh}2\psi `$, represents the contribution of the zero mode. The remaining terms describe the contribution from the parts of the phase space, where both the zero mode and ballistic solutions for the electron-hole wave pairs are intermixed. In the semiclassical limit only the first term in the first line of Eq. (99) survives with $`N_1^{}\mathrm{}`$, and thus coincides with Eq. (64). In the quantum limit, $`\tau _\mathrm{E}0`$, only the term containing $`\mathrm{tanh}2\psi `$ remains in Eq. (99), and corresponds to the random matrix result of Eq. (47).

According to Eq. (99), the density of states is always finite at $`\epsilon >\pi /4\tau _\mathrm{E}`$, and is formed by at least the second term of the first line of Eq. (82). We investigate if the gap energy $`E_\mathrm{g}`$ is smaller than $`\pi /4\tau _\mathrm{E}`$, so that the density of states originates from electron states described by the zero mode.

From Eq. (99) we conclude that the density of states also exists below $`\pi /4\tau _\mathrm{E}`$ as long as the solution to Eq. (98) has an imaginary part. For small energies $`\epsilon `$, Eq. (98) has real solutions and $`𝒩(\epsilon )=0`$. As energy $`\epsilon `$ increases, two solutions approach each other and at certain value of energy $`\epsilon =E_\mathrm{g}`$, they are both equal to some value $`\psi _\mathrm{c}`$. The energy $`E_\mathrm{g}`$ is the gap energy and the density of states has a square root singularity just above the gap. Indeed, we can expand both sides of Eq. (98) to the lowest order in $`\psi \psi _\mathrm{c}`$ and $`\epsilon E_\mathrm{g}`$. The expansion in $`\epsilon E_\mathrm{g}`$ linear, but expansion in $`\psi \psi _\mathrm{c}`$ is quadratic, since the first order derivative of Eq. (98) with respect to $`\psi `$ vanishes at $`\psi _\mathrm{c}`$, see Sec. III.2.1 for a similar analysis.

Figure 10 demonstrates dependence of the energy gap $`E_\mathrm{g}`$ on the Ehrenfest time $`\tau _\mathrm{E}`$, shown by a solid line. To the lowest order in $`\gamma _{\mathrm{esc}}\tau _\mathrm{E}`$ the gap energy decreases linearly from the random matrix result,MBFB  see also Eqs. (48) and (49):

$$E_\mathrm{g}E_\mathrm{g}^{\mathrm{rmt}}(10.23\gamma _{\mathrm{esc}}\tau _\mathrm{E}).$$
(100)
The density of states above $`E_\mathrm{g}`$ is given by

$$𝒩(\epsilon )=\frac{c}{\delta _1}\sqrt{\frac{\epsilon }{E_\mathrm{g}}1},$$
(101)
where $`cc_{\mathrm{rmt}}(10.42\gamma _{\mathrm{esc}}\tau _\mathrm{E})`$. As the Ehrenfest time increases, the energy gap decreases and is given by the following asymptote at $`\tau _\mathrm{E}\gamma _{\mathrm{esc}}1`$:

$$E_\mathrm{g}\frac{\pi }{4\tau _\mathrm{E}}\left(1\frac{1}{\tau _\mathrm{E}\gamma _{\mathrm{esc}}}\right).$$
(102)
We plot the asymptote of Eq. (102) by a dotted line in Fig. 10.

## V Discussion and Conclusions

In this paper we considered the effect of quantum disorder and quantum chaos on electron density of states in an Andreev billiard - a metal grain or a semiconductor dot of irregular shape, connected to a superconductor.

First we restate our assumptions about the considered system and briefly review prior publications, discussing related systems. We study the density of states at energies much smaller than the gap energy in the superconductor and the Thouless energy $`E_{\mathrm{Th}}`$ (or the inverse ergodic time $`\tau _{\mathrm{erg}}`$). The opposite limit, when the escape time is equal or less than the ergodic time, was considered in refs. \[not0d, \] for Andreev billiard, in refs. \[PBB, ; OSF, ; 1dU, \] for one dimensional systems (wires) and in refs. \[SLL, ; LO, \] for bulk superconductors.

In our calculations we use the exponential distribution $`P(l)`$ of trajectory lengths $`l`$, described by Eq. (59). An analysis of motion in classically chaotic systems shows deviation of the distribution function from the exponential law due to anomalously long trajectories.ZE  The effect of these trajectories on the density of states of an Andreev billiard was studied in ref. \[IR, \] within semiclassical theory, when the contribution of long trajectories is especially important. We assume that the deviation from the exponential distribution occurs at lengths $`l`$ longer than either the Ehrenfest time or the mean free time, and we neglected the contribution of anomalously long trajectories to the electron density of states.

We study only non-integrable systems. The difference between integrable and chaotic non-integrable systems was examined for a quantum dot in refs. \[MBFB, ; CBA, \].

We do not consider the effect of finite mean level spacing in the grain. The finite mean level spacing leads to mesoscopic fluctuations of the energy levels, and particularly, to the fluctuations of the gap.VBAB  The fluctuations of the energy gap smear the square root singularity of the density of states, Eq. (48b), and produce an exponentially small tail of the ensemble average density of states below the gap energy, see refs. \[BNA, ; VBAB, ; OSF, ; OSF1, \]. Nevertheless, the corresponding scale for the gap fluctuations and the mean level spacing near the gap are determined by the following small energy scale: $`(\delta _1^2\gamma _{\mathrm{esc}})^{1/3}\gamma _{\mathrm{esc}}`$. Significant fluctuations of the first energy level (energy gap) occur in exponentially rare cases.BNA ; VBAB ; OSF ; OSF1

We assume that the reflection at the contact with the superconductor is described by the Andreev mechanism, i.e. the contact has no tunnel barrier between the grain and the superconductor. The effect of finite barrier was studied in ref. \[MBFB, \] in the random matrix limit, where it was shown that a weak (normal) reflection at the barrier does not qualitatively change properties of the curve of the density of states.

We do not take into account the Coulomb interaction between electrons in the grain responsible for the effects of zero bias anomaly and the Coulomb blockade. In our system, the conductance of the SN contact is large (the number of channels $`N1`$), and the Coulomb interaction may be disregarded. The influence of the zero bias anomaly on proximity effect was studied in ref. \[OBSA, \], and the Coulomb blockade in a small metal grain connected to a superconductor was studied in ref. \[MatGlaz, \].

we also assume that the repulsion in the Cooper channel is sufficiently small and we disregard it. The effect of repulsion was studied, for example in refs. \[Rep1, ; Rep2, ; Rep3, \].

The model of clean chaotic Andreev billiards, considered in Sec. IV of our paper, was also recently addressed in paper \[Bnew, \]. According to ref. \[Bnew, \], the energy gap and the density of states above the gap are determined by the random matrix theory resultMBFB  with the appropriately modified number of channels and the mean level spacing. Although the resultsBnew  qualitatively coincide with ours, they are quantitatively differ from our results, presented in Eqs. (100)-(102) and in Fig. 10. We believe that the calculationsMBFB  within random matrix theory are justified only as $`\tau _\mathrm{E}0`$, since these calculations take into account only the zero mode. In general, as we demonstrated in Sec. IV, the ballistic parts of trajectories at length smaller than $`\tau _\mathrm{E}`$, also known as Hikami boxes, give rise to a substantial contribution to both the density of states and the self-consistency equation. Consequently, Eqs. (98) and (99) cannot be reduced by a change of mean level spacing and the number of channels to the random matrix theory equations of ref. \[MBFB, \], see also Eqs. (46) and (47) in Sec. III.2.1, and thus the quantitative result of ref. \[Bnew, \] is questionable.

The most interesting qualitative result of our work is the oscillations of the density of states as a function of $`\epsilon \tau _\mathrm{E}`$. \[An oscillating character of the density of states near the gap was also found in ref. \[AB, \].\] These oscillations are related to the contribution to $`𝒩(\epsilon )`$ from the ballistic part of the electron Green function. Although the full length of different trajectories significantly varies, the Ehrenfest time is nearly the same for all long trajectories with length $`l>\tau _\mathrm{E}`$. Indeed, even strong fluctuations of the Lyapunov exponent are averaged out on time scale of the ergodic time $`\tau _{\mathrm{erg}}`$, resulting in weak fluctuations of the Ehrenfest time. The fluctuations of the Ehrenfest time are important only at parametrically large energies $`\epsilon \sqrt{\lambda /\tau _\mathrm{E}}`$ and lead to suppression of the oscillations. Similar oscillations also exist in frequency dependence of the weak localization correction to the conductivityAL1  and the correlation function of electron energy levels.AL2

We notice that the oscillations discussed above should not be associated with the oscillations of the energy levels in the Wigner-Dyson statistics. The latter oscillations have a much smaller period (equal to the mean level spacing) than $`\tau _\mathrm{E}^1`$. The energy level oscillations were studied in systems with quantum disorderEfetov  and with quantum chaos.Alt ; BK  Similar oscillations are also present in the density of states near the gap with the period $`\delta _1^{2/3}\gamma _{\mathrm{esc}}^{1/3}`$.VBAB ; OSF  In the present paper we do not take into account these oscillations, as well as other mesoscopic effects.

The observation of oscillations in the density of states with period $`\tau _\mathrm{E}^1`$ is the major distinguishing characteristic of quantum chaotic systems from disordered systems, since in quantum disorder such oscillations are absent. The period of the oscillations may be used to estimate the Ehrenfest time and the Lyapunov exponent in each particular metal grain.

In conclusion, our results allow us to describe the classical-to-quantum crossover in the density of states of Andreev billiards. Methods, developed in this paper should also be applicable to study this crossover in many other problems of solid state physics and optics.

## Acknowledgements

We are grateful to I. Aleiner for many fruitful discussions. We also thank P. Brouwer for important suggestions and I. Adagideli and C. Beenakker for their comments. The work was supported by NSF grant # 0120702.

## Appendix A

In this appendix we show how the small angle scattering term, which mimics the effect of quantum diffraction, may be used to obtain the results of Sec. IV.2.1.

A small angle scattering potential, with a typical scattering angle $`\varphi _{\mathrm{min}}\sqrt{\lambda _\mathrm{F}/L}`$, leads to the following form of the self energy part of the Eilenberger equation:

$`\widehat{\mathrm{\Sigma }}(\epsilon ,𝐫,𝐧)`$ $`=`$ $`{\displaystyle \frac{1}{\tau _0}}\widehat{𝒢}(\epsilon ,𝐫,𝐧)+{\displaystyle \frac{1}{2\tau _\mathrm{q}}}_𝐧^2\widehat{𝒢}(\epsilon ,𝐫,𝐧),`$ (103)
where we introduced the full mean free path $`\tau _0`$,

$`{\displaystyle \frac{1}{\tau _0}}`$ $`=`$ $`2\pi \nu n_\mathrm{i}{\displaystyle |V(𝐧𝐧^{})|^2\frac{d𝐧}{\mathrm{\Omega }_d}},`$ (104a)
and the transport time $`\tau _\mathrm{q}`$,
$`{\displaystyle \frac{1}{\tau _\mathrm{q}}}`$ $`=`$ $`2\pi \nu n_\mathrm{i}{\displaystyle |V(𝐧𝐧^{})|^2|𝐧\times 𝐧^{}|^2\frac{d𝐧}{\mathrm{\Omega }_d}}.`$ (104b)
The first term in the right hand side of Eq. (103) does not contribute to the scattering term of the Eilenberger equation. The second term represents small angle diffraction and the Eilenberger equation acquires the form with the Perron-Frobenius differential operator. For the two dimensional electron system $`_𝐧=/\varphi `$, where $`\varphi `$ is the angle in the direction of $`𝐧`$. we estimate $`\tau _\mathrm{q}`$ and $`\tau _0`$ below.

We use the Ricatti parameterization of the Green function, see Eq. (15). For the scattering term of the Eilenberger equation in the form of Eq. (103), we rewrite Eqs. (16) as:

$`\left(2i\epsilon \right)a_+`$ $`=`$ $`{\displaystyle \frac{1}{\tau _\mathrm{q}}}\left(_𝐧^2a_+{\displaystyle \frac{2a_{}[_𝐧a_+]^2}{1+a_+a_{}}}\right),`$ (105a)
$`\left(2i\epsilon +\right)a_{}`$ $`=`$ $`{\displaystyle \frac{1}{\tau _\mathrm{q}}}\left(_𝐧^2a_{}{\displaystyle \frac{2a_+[_𝐧a_{}]^2}{1+a_+a_{}}}\right)`$ (105b)
Here again we omitted variables $`\{\epsilon ,𝐱\}`$ of the functions $`a_\pm =a_\pm (\epsilon ,𝐱)`$.

We notice that the average of the electron-hole wave function over the phase space is small at large energy $`\epsilon \gamma _{\mathrm{esc}}`$, see Eq. (74). This observation allows us to neglect non-linear terms in Eqs. (105), since these terms would contribute to the higher than the second order in $`\gamma _{\mathrm{esc}}/\epsilon `$. We write down equation for $`A_2^\pm (\epsilon ,𝐱_1,𝐱_2)`$ using the linearized equations of motion for $`a_\pm (\epsilon ,𝐱)`$, Eqs. (105).

$`[4i\epsilon _{\mathrm{cm}}(\rho ,\varphi )+{\displaystyle \frac{1}{2\tau _\mathrm{q}}}_\varphi ^2]A_2^\pm (\epsilon ,𝐑,𝐍,\rho ,\varphi )=0,`$ (106)
where $`_{\mathrm{cm}}`$ and $`(\rho ,\varphi )`$ are defined by Eqs. (77). As we already mentioned, the non-linear terms in the right hand side of Eqs. (105) lead to higher terms in $`\gamma _{\mathrm{esc}}/\epsilon `$.

After averaging over the position of the center of mass in the phase space, we obtain a counterpart of Eq. (78):

$$\left(4i\epsilon +s_2\gamma _{\mathrm{esc}}_\rho +\frac{1}{2\tau _\mathrm{q}}_\varphi ^2\right)A_2^\pm (\epsilon ,\rho ,\varphi )=\gamma _{\mathrm{esc}}.$$
(107)
Below we assume that $`s_21`$. Introducing new variables $`z`$ and $`\chi `$, according to Eq. (79), and then averaging over $`\chi `$, we arriveAL1  to the following equation for $`A_2^\pm (\omega ,z)`$ with $`\omega =4\epsilon i\gamma _{\mathrm{esc}}`$:

$$\left[i\omega +\lambda _z+\frac{\lambda _2}{2}_z^2+\frac{e^{2z}}{2\tau _\mathrm{q}}\left(\frac{1\gamma }{2}_z^2+\gamma _z\right)\right]A_2^\pm =\gamma _{\mathrm{esc}}$$
(108)
where $`_z=/z`$, and $`\gamma `$ is a numerical coefficient of order of unity.AL1  Equation (108) should be supplied with the boundary conditions at $`z=\mathrm{ln}b/L`$, determined by Eq. (81). The solution can be found as a superposition of two terms, see Eq. (82). One term, $`A_{2,\mathrm{b}}^\pm (\epsilon ,z)`$, is given by Eq. (83) and the other term is given by

$`A_{2,\mathrm{q}}^\pm (\epsilon ,z)`$ $`=`$ $`\left[\overline{a_\pm (\epsilon ,𝐱)}^2A_{2,b}^\pm (\epsilon )\right]{\displaystyle \frac{w(\omega ,z)}{w(\omega ,\mathrm{ln}b/L)}}.`$ (109)
Function $`w(\omega ,z)`$ is a solutionAL1  of homogeneous Eq. (108):

$$w(\omega ,z)=\mathrm{exp}\left[\left(\frac{i\omega }{2\lambda }+\frac{\lambda _2\omega ^2}{4\lambda ^3}\right)\mathrm{ln}\frac{\lambda \tau _\mathrm{q}}{\lambda \tau _\mathrm{q}e^{2z}+\gamma /2}\right].$$
(110)
Now solutions for $`A_2^\pm (\epsilon ,z)`$ are regular functions as $`z\mathrm{}`$, and correspondingly the product of two electron-hole pairs, see Eq. (75), can be found at coincident points.

An exact form and strength of small angle scattering potential is not important, since as we will show later, the density of states depends on the strength of the potential only as a logarithm of the ratio of the transport time of this potential and the size of the SN contact. We will choose the strength of the scattering potential so that at small time this potential scattering would describe the quantum mechanical diffraction effect. The scattering time $`\tau _\mathrm{q}`$ may be estimated as $`\tau _\mathrm{q}=\tau _0\varphi _{\mathrm{min}}^2`$, compare Eqs. (104a) and (104b). We assume that $`\lambda \tau _01`$ (since both $`\lambda `$ and $`\tau _0`$ are characterized by the same potential energy of electrons in the dot) and then we estimate $`\lambda \tau _\mathrm{q}\varphi _{\mathrm{min}}^2=L/\lambda _\mathrm{F}`$, see Eqs. (104a) and (104b). Combining Eqs. (81) and (109) and we obtain Eq. (88).

In conclusion, we discuss the meaning of Eq. (110). Within logarithmic accuracy, the term $`\lambda \tau _\mathrm{q}e^{2z}+\gamma /2`$ can be replaced by $`\mathrm{max}\{\lambda \tau _\mathrm{q}e^{2z},\gamma /2\}`$. At large scales $`\rho `$ and $`\varphi `$ of a wave packet, $`z>z_\mathrm{c}=\mathrm{ln}\sqrt{2\lambda \tau _\mathrm{q}/\gamma }`$, we can disregard the quantum diffraction, described by the last term in Eq. (106). In this “classical” approximation we find

$$w(\omega ,z_\mathrm{c})=\mathrm{exp}\left[\left(\frac{i\omega }{2\lambda }+\frac{\lambda _2\omega ^2}{4\lambda ^3}\right)\mathrm{ln}\frac{2\lambda \tau _\mathrm{q}}{\gamma }\right].$$
(111)
At smaller scales of the wave packet, $`z<z_\mathrm{c}`$ or $`\rho ^2+L^2\varphi ^2\rho _{\mathrm{min}}^2`$, the quantum diffraction overpowers the Lyapunov classical divergence of trajectories. Consequently, the function $`w(\epsilon ,z)`$ weakly depends on the scale $`z`$, see Eq. (110), Eq. (110), and we use $`w(\epsilon ,z\mathrm{})w(\epsilon ,z_\mathrm{c})`$. It is exactly this property of the function $`w(\epsilon ,z)`$ was used in Sec. IV.2.1 of the present paper, when we used $`w(\epsilon ,z_\mathrm{c})`$

In this Appendix we used linearized Eqs. (105), which are justified at high energy $`\epsilon \gamma _{\mathrm{esc}}`$. In general, functions $`A_2^\pm (\epsilon ,z)`$ are described by non-linear Eqs. (105). Nevertheless, the right hand side of Eqs. (105), which contains non-linear terms, is necessary only to determine the function $`w(\epsilon ,z)`$ at small scale of the wave packet, i.e. for $`z<z_\mathrm{c}`$, and at already at scale $`zz_\mathrm{c}`$ the right hand side of Eqs. (105) does not affect the solution for $`w(\epsilon ,z)`$ within logarithmic accuracy. That is why the non-linear term is not important within logarithmic accuracy.

Within the same logarithmic accuracy, we use approximate solutions for $`A_k^\pm (\epsilon ,z)`$ at $`zz_\mathrm{c}`$ for $`k2`$ in Sec. IV.2.