{'fraction_non_alphanumeric': 0.05083101026853453, 'fraction_numerical': 0.013620805250714563, 'mean_word_length': 3.9553243574051407, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weak Charge Quantization on Superconducting Islands

We consider the Coulomb blockade on a superconductive quantum dot strongly coupled to a lead through a tunnelling barrier and/or normal diffusive metal. Andreev transport of the correlated pairs leads to quantum fluctuations of the charge on the dot. These fluctuations result in exponential renormalization of the effective charging energy. We employ two complimentary ways to approach the problem, leading to the coinciding results: the instanton and the functional RG treatment of the non-linear sigma model. We also derive the charging energy renormalization in terms of arbitrary transmission matrix of the multi-channel interface.


I. INTRODUCTION
Physics of interacting electronic systems in presence of disorder has been a subject of an intense study already for a few decades [1]. Various theoretical approaches have been developed for the description of both metallic and insulator phases. The non-linear σ-model (NLσM) in the replica [2,3] or dynamic (Keldysh) [4,5] formulation has proven to be the most powerful tool to deal with the weakly disordered (metallic) phase. It was shown that both versions of the σ-model may effectively treat the perturbation theory as well as the renormalization group (RG) formalism. Unlike the non-interacting models, where a whole specter of non-perturbative results is available [6], it was relatively little progress in the development of non-perturbative solutions of the interacting NLσM. Our goal is to make a step in this direction, using the Coulomb blockade (CB) on a superconductive (SC) quantum dot as a prototypical example.
The Coulomb blockade on a quantum dot coupled to a certain number of leads proved to be an extremely rich and fascinating phenomena both theoretically and experimentally (see Refs. [7,8] for review). From the theoretical point of view it provides a model, where the Coulomb interactions, being spatially localized, may be treated in a non-perturbative way. The interactions strongly affect charge (and spin) fluctuations between the dot and the leads, that manifests in the peculiar transport and thermodynamic behavior of the coupled dot-lead system. In the case of practically isolated dot, such that the dimensionless conductance (measured in units of G Q ≡ e 2 /(2πh)) of the dot-lead interface is small (G ≪ 1), the fluctuations usually may be taken into account perturbatively. (The notable exceptions are provided by the vicinity of the charge degeneracy point and by the Kondo effect on the dot [8]). Here we concentrate on the opposite scenario of the dot strongly connected to the leads (G > 1). In this case the Coulomb blockade is expected to be suppressed by the charge fluctuations and the overall effect of interactions on a single dot to be weak. It is, however, a challenging theoretical problem, lacking an obvious small parameter, to understand a remnants of the CB on a strongly connected dot. More importantly the weak CB may prove to be a strong phenomena in the granulated systems, with many (superconductive) dots incorporated in the conducting matrix.
For the normal dot the weak CB effect was relatively well understood from the various points of view. Matveev [9] gave a complete picture of the phenomena for the case of one or two conducting channels connecting the dot and the leads. For the multi-channel system the CB suppression was calculated with the exponential accuracy employing RG technique [10,11], the instanton calculus [12][13][14][15] and bosonization [8]. The factor, suppressing the thermodynamic CB oscillations, was formulated by Nazarov [14] in terms of transmission coefficients, T α , of the dot-lead interface as α (1 − T α ) 1/2 , where α = 1 . . . N and N is the number of spinless channels. The remarkable feature of this result is that the essentially many-body phenomena may be described via the single-particle (non-interacting) scattering matrix only. Yet the knowledge of the interface conductance G = α T α alone is not sufficient to describe the CB oscillations. For example, in a dot coupled to a lead through the tunnelling barrier with T α ≪ 1, one obtains for the CB suppression factor exp{−G T /2}. Another important case of a dot-lead interface is a coupling via a piece of diffusive metal. Employing Dorokhov statistical distribution [16] of the transmission coefficients P (T ) = G D /(2T √ 1 − T ), where G D is the conductance of the diffusive area, one obtains for the typical CB suppression factor exp{−π 2 G D /8}. Since the result is exponentially sensitive -the difference between e.g. G/2 and π 2 G/8 may be actually enormous.
The superconductive dot in contact with normal leads is even more challenging system. We consider the low temperatures T ≪ ∆, where ∆ is the SC gap, and therefore only very few (if at all) quasi-particle excitations are allowed to leave or enter the superconductor. The dominant mechanism of the charge transfer through the interface is thus the Andreev transport of the correlated pairs. The Andreev transmission of a given channel is known to be T α = T 2 α /(2 − T α ) 2 [17,18]. For the tunnelling barrier setup (T α ≪ 1) this leads to an overall Andreev conductance that scales like G A ∼ G 2 T /N (pair tunnelling probability). In most cases this is a very small number. The presence of the diffusive normal metal adjacent to the tunnelling barrier increases the Andreev conductance to G A ∼ G 2 T /G D (or, in the case G T > G D , to G A ∼ G D ). The physical reason of this phenomena is multiple attempts of Andreev transmission due to the coherent back-scattering on the normal impurities. The natural question is whether it is G or G A (or may be none of them) that determines the amplitude of the CB oscillations. In case G A is the relevant quantity, one may wonder whether the coherent back-scattering enhancement should be taken into account. The answers are not immediately obvious, since while the SC strongly prefers the pair transport -the Coulomb energy of the dot makes the entrance of two charges at once energetically costly. The interactions of the dot may also provide a dephasing mechanism which ruins coherent back-scattering. We show that parametrically it is indeed G A (including diffusive enhancement) which determines the CB amplitude. The coefficient in the exponent, however, is setup dependent, making the final answer sensitive to the ratio G T /G D and not only to G A . The effective charging energy turns out to be proportional to α (1 − T α ) 1/2 , where T α are the Andreev transmissions defined above. This expression provides a remarkable analogy between the normal and Andreev results.
The other interesting phenomena brought by superconductivity is the parity effect [19,20]. At very small temperature the pair tunnelling is the dominant mechanism and the period (in the gate voltage) of the CB oscillations is twice larger than in the normal dot. For a closed dot the parity phenomena is destroyed by the entropic effects at moderately small temperature T * = ∆/ ln(∆/δ) ≪ ∆, where δ is the mean single-particle level spacing on the dot [20]. The physical reason for this temperature to be much less than ∆ is that it is enough to have a single excited quasiparticle to destroy the parity effect. At larger temperature the system exhibits the "normal" oscillation period. We show that for an open dot the transition temperature between the normal and doubled period is somewhat larger than for a closed dot and is given by Technically we treat the problem from the two seemingly distinct perspectives. First we look for the spatiallydependent instanton solution of interacting NLσM on SC dot in contact with the normal diffusive region. The finite action of such instanton configuration results in the exponential suppression of CB amplitude. Alternatively we can treat the problem employing the functional RG technique. The latter approach encodes the entire dot-lead interface in an (infinite) set of coefficients, whose values are subsequently renormalized by the quantum fluctuations of the phase on the dot. Renormalization is terminated at the cut-off energy scale where the conductance reaches unity. The renormalized (exponentially small) cut-off energy is the effective charging energyẼ C , this dictates the amplitude of the CB oscillations. One of the messages of the present paper is that these two approaches lead to the identical results.
The paper is organized as follows: in section II we describe the setup and formulate the appropriate action in the interacting NLσM language. Section III is devoted to the instanton treatment of the problem in the real space. In section IV we derive the proximity action functional and obtain CB amplitude for an arbitrary set of transmission eigenvalues (using instanton approach) at moderately low temperaturesẼ C ≪ T ≪ E C . Then in section V the RG approach is employed to treat the CB in the same temperature range. It is demonstrated that the results coincide with those obtained by instanton techniques. Finally in section VI we discuss the physical results and their possible experimental signatures. Appendix contains derivation of the instanton action starting from the real-time Keldysh functional technique.

II. PROBLEM SETUP AND ACTION
We consider a large diffusive (or chaotic) SC dot. The mean single-particle level spacing of the dot, δ, is supposed to be the smallest energy scale in the problem. The SC gap, ∆, on the other hand, is the largest scale. The electrons on the dot interact via the capacitive interaction of the form whereN is the electron number operator and q is the rescaled dimensionless gate voltage potential. The charging energy, E C = e 2 /(2C), is assumed to satisfy the inequalities δ < E C < ∆. The dot is separated from the normal diffusive metal by the tunnelling barrier with the conductance G T , see Fig. 1. The piece of the quasi 1D or 2D metal having size L and conductance G D is in turn connected to a clean bulk 3D lead. By the reasons which are explained later we shall assume that the Thouless energy of the diffusive region, E Th = D/L 2 , where D is the diffusion coefficient, is larger than the charging energy of the dot, E Th > E C . The opposite limiting case requires a separate treatment and will be presented elsewhere.
We shall be interested in those thermodynamic characteristics of the dot, which depend in an oscillatory way on the gate voltage potential, q. Specifically, we look for the free energy F (q) = −T ln Z(q), where Z(q) is the partition function. Both of these quantities depend on a particular realization of disorder in the diffusive region and inside the dot. We shall therefore look for the disorder averaged quantity like F (q) . For an open dot F (q) is an oscillatory function, with a phase sensitive to a disorder realization. It is therefore convenient sometimes to calculate also the correlation function F (q)F (q ′ ) [15]. The latter carries information about the typical F (q), rather than the average one (that may spuriously vanish due to phase randomness).
The disorder averaging may be performed in two ways either by introducing p replica of the system [2,3] and sending p → 0 in the end, or by dealing with the dynamical Keldysh formulation [4,5]. We shall employ both of these approaches to demonstrates what they are consistent and interchangeable. In the replica formalism, the NLσM is formulated in terms of the matrix field theory of the matrix field Q ab ij (r, τ, τ ′ ), where a, b = 1, 2 . . . 2p are the replica indexes (one needs 2p replica to describe the correlation function Z p (q)Z p (q ′ ) ) and i, j = 1, 2 are Gorkov-Nambu indexes. The correlation function of the free energies is given by F (q)F (q ′ ) = lim p→0 p −2 ( Z p (q)Z p (q ′ ) − 1). The Q-matrix obeys the constraint Q 2 = 1 and the fermionic anti-periodic boundary condition in both of its imaginary time arguments τ, τ ′ ∈ [0, β]. The coordinate r runs over the volume of the SC dot and the normal diffusive region. The matrix Q-field describes dynamics of the electrons; it is coupled to a scalar bosonic vector field Φ a (τ ) which originates from the Coulomb interactions on the dot. Since we restrict ourselves by the simplest capacitive interaction Eq. (1), the field Φ a (τ ) is space independent throughout the dot and vanishes outside the dot. As a result, only spatially independent (q = 0) component of the superconductive Q-matrix of the dot appears to be coupled to the Coulomb field Φ. Therefore we assume that the Q-matrix is spatially constant inside the dot, while it may have a non-trivial spatial structure inside the normal diffusive region. The effective action for our geometry contains three terms where S dot and S D are the effective actions of the isolated SC dot and normal diffusive region correspondingly, while S T is the tunnelling action which provides the coupling between the two. We shall examine separately these three contributions.
A. Action of the dot According to our model -the dot is the only region, where electrons interact via the Coulomb interaction Eq. (1). As a result, the dot's action contains the two coupled fields: Φ a (τ ) and Q ab S (τ, τ ′ ). The first one originates from the Hubbard-Stratonovich decoupling of the Coulomb term [21], while the latter is the result of disorder averaging and integration of the fermionic degrees of freedom on the dot [2]. Both of these fields are spatially independent, reflecting the fact that the Thouless energy of the dot is assumed to be large. The action takes the standard by now form: Here the replica vector of gate voltages, q a , is defined in a way that q a = q for a ∈ [1, p] while q a = q ′ for a ∈ [p+ 1, 2p]. The SC order parameter∆ = (∆ a σ + − ∆ * a σ − )δ ab , is the same in all replica and is written as a matrix in the Gorkov-Nambu space, where σ ± = (σ 1 ± iσ 2 )/2. Notice, that the gate voltage dependence of the partition function originates entirely from exp{iπ a q a W a }, where πW a ≡ dτΦ a (τ ) = Φ a (β) − Φ a (0) is the zeroth Matsubara component of thė Φ a (τ ) field.
Since all the energy scales we consider are larger than δ, one may evaluate the second term in the dot's action (3) in the saddle point approximation over the field Q S . In so doing, one disregards the mesoscopic conductance fluctuations of the dot-lead interface [8]. The saddle point value of the Q S field is given by the Gorkov Green function gauged by the phase Φ(τ ): The Gorkov Green function, Λ S , has the standard form, which is more familiar in the Matsubara basis (we assume the phase of the SC dot without the Coulomb interactions to be zero) Λ ab S (n, m) = δ ab δ nm cos θ n sin θ n sin θ n − cos θ n , where cos θ n ≡ ǫ n and ǫ n = πT (2n + 1) is the fermionic Matsubara frequency. The phase rotation (4) preserves the fermionic antiperiodic boundary conditions if all W a are even integers. At small temperatures, T ≪ ∆, however, one has cos θ n ≈ 0 and sin θ n ≈ 1 and therefore the Gorkov matrix, Λ S is (almost) off-diagonal in the Nambu space and local in time. As a result, only 2Φ(τ ) participate in the phase rotation, Eq. (4). Therefore the odd-integer W a preserve the fermionic boundary conditions, as well. The odd-integers W a result in the doubling of the period of the F (q) function with respect to the normal case. This reflects the fact that the pair transfer is the dominant mechanism of the charge exchange.
In fact, one has to be more careful and recall that, according to Eq. (3), one has to perform integrations over all Matsubara components ofΦ a (τ ) fields [21]. All non-zero Matsubara components may be eliminated by the gauge transformation, Eq. (4). As a result, they have no effect on the thermodynamic of an isolated dot at all (they are of major importance, of course, once the dot is coupled to the leads). The remaining (usual) integral over the zeroth Matsubara component, W a , must be performed explicitly. To this end one notices that πT W a enters the action as an imaginary chemical potential in the replica a. The (replicated) free energy of an isolated dot is thus a periodic function of each of W a with the period 2 (indeed the chemical potential always enters as exp{µN /T }, since the number operator,N , has only integer eigenvalues -the periodicity is apparent). The free energy possesses deep minima at even integer values of W a with the quadratic behavior in their vicinity F dot (W ) ≈ π 2 T 2 W 2 /(2δ), where δ = (∂ 2 F dot /∂µ 2 ) −1 is the mean level spacing of the dot. For a sufficiently large dot, where δ < T , the integrals over W a may be performed in the saddle point approximation, which results in the even-integer quantization of W a . In a SC dot there are additional minima at odd-integer values of W a . Consequently the integration over W a is substituted by the summation over all integers. At T > 0 the additional minima of the free energy at odd integer W a are not as deep as at the even integers, reflecting the fact that the addition of an odd number of electrons is possible by creating a quasi-particle. We shall evaluate now the action cost of the odd minima with respect to the even ones.
For the even values, say W a = 0, one substitute the saddle point solution, Eq. (4), into the action, Eq. (3), and obtains (we disregard for a moment the first term in Eq. (3)) This sum is divergent. However, it is only the difference of the action between different replica that has a physical significance. The later quantity is convergent as we shall see momentarily. For the odd integers, say W b = 1, there is an imaginary component iπT of the chemical potential in Eq. (3). It may be eliminated by the gauge transformation which converts the antiperiodic boundary conditions for the fermions into the periodic one. As a result, one arrives to the same expression as Eq. (7) with the fermionic Matsubara frequency ǫ n = πT (2n + 1) substituted by the bosonic one ω n = 2πT n: Employing the Poisson summation formula, one obtains for the difference where the last equality assumes T ≪ ∆. In the opposite limit of the normal dot, T > ∆, one finds A(T ) ≈ π 2 T /δ, reflecting the fact that in a normal dot the odd-integer minima are absent. Although these minima persist up to T ≈ ∆, their contribution is exponentially suppressed at T > T * , where T * is determined from the condition As a result, there is an important temperature dependence associated with the odd-integer winding numbers at the scale T ≈ T * ≪ |∆|. We shall see below that for the dot strongly coupled to the leads, the corresponding temperature scale is slightly different.
We summarize now our discussion of the action of a large SC dot with the Coulomb interaction at T ≪ ∆. The scalar potential in each replica obeys the boundary condition where W a is an integer winding number. The corresponding saddle point value of the Q S -matrix field is given by Eq. (4). The action of the dot takes the form: where A(T ) is given by Eq. (9).

B. Tunnelling barrier and diffusive region action
The tunnelling action couples the Q S -field on the dot with the Q(r = 0)-field at the point r = 0 adjacent to the tunnelling barrier from the normal metal side. It has the standard form [6] where G T is the tunnelling conductance. The action of the normal diffusive region also has the standard form [2,6] where D is the diffusion constant and ν the density of states (per spin) of the normal region. The total conductance of the normal region is given by G D = 4πνD/L for the quasi 1D geometry of the normal region and G D = 8π 2 νD/ ln(L/d) for the 2D geometry. Here L is the length (radius) of the 1D (2D) region and d is the radius of the SC dot. We have omitted the frequency term ν Tr{ǫQ} on the r.h.s. of Eq. (14), because of the assumption that E Th > E C > T . At the point where the diffusive region is attached to the normal bulk lead one has to impose the boundary condition where Λ ab N (n, m) = σ 3 δ ab δ nm sign(ǫ n ) is the appropriate Q-matrix of the normal bulk lead. Alternatively, one may imagine integrating out the Q(r)-field of the normal region, subject to the boundary condition Eq. (15) and weighted by the action S T + S D . This procedure (we shall describe it in details in section IV) leads to the effective action S T D of the interface plus diffusive region written in terms of Q S and Λ N only. If all the relevant energy scales are less than the Thouless energy of the diffusive metal -the general form such action may take is [11] where γ l are coefficients which depend on the details of the interface (in our case the ratio G T /G D ). The largeness of the Thouless energy is necessary to disregard the retardation effects and thus possible time non-local coupling between Q S and Λ N . Under such condition, the proximity action, Eq. (16), is completely equivalent to those given by Eqs. (13) and (14) upon the proper choice of the set, γ l [11].

III. THE INSTANTON APPROACH
We are interested in the limit of strong coupling between the dot and the leads, meaning G D , G T > 1 (the weak coupling limit may be treated in the spirit of Refs. [22][23][24]). For G D ≫ 1 the fluctuations of the Q-field around its optimal value are suppressed. One may employ therefore the stationary phase treatment of the NLσM for the dot-lead interface [15]. Taking the variation of the action Eqs. (13), (14) under the condition Q 2 = 1 one obtains the Usadel equation This equation is to be solved for a fixed Q S = Q S [Φ] given by Eq. . The later may then be investigated using the instanton approach applied to the Φ(τ )-field.

A. Zero winding number
As a warm-up exercise, consider the zero winding number sector of the theory, W a = 0 for a = 1, . . . 2p. Obviously it does not produce an oscillatory dependence of F (q), cf. Eq. (12), and therefore serves only an axillary purpose. The lowest energy configuration in zero winding number sector is Φ a (τ ) = 0 and therefore Q S on the dot is simply given by the BCS Λ S , Eq. (5). Solution of the Usadel equation (17) may be written as where u(r) is the normalized "voltage drop" inside the normal region: u(r) = (L − r)/L in 1D and u(r) = ln(L/r)/ ln(L/d) in 2D. The absolute value of the voltage drop, θ, is coordinate independent diagonal (in replica and Matsubara space) matrix. It has a physical meaning of the SC rotation angle of the normal metal in the direct proximity to the dot. One may substitute the solution back into the action to obtain (for a single replica and Matsubara component) where the subscript "0" stresses that we work with zero winding numbers and θ n is defined by Eq. (6). For small energy, ǫ n ≪ ∆, one has θ n ≈ π/2 and therefore the corresponding action takes the form where t ≡ G T /G D . This action is minimized when θ = θ(t) satisfies the equation The lowest energy solution of this equation smoothly interpolates between θ = 0 for t ≪ 1 and θ = π/2 for t ≫ 1.
Finally, this solution has to be substituted into Eq. (20) to find the action cost, S 0 = S 0 (t), for the zero winding number configuration.

B. Non-zero winding numbers
To calculate the oscillatory component of the free energy, F (q), one has to consider the Φ-field configurations with non-zero winding numbers, cf. Eq. (12). Consider, thus, the simplest even configuration of winding numbers with W 1 = 2 and all others W a = 0 in the remaining 2p − 1 replica. With the exponential accuracy it is sufficient to consider the "straight" windings: Φ 1 (τ ) = 2πT τ . The saddle point of the SC Q-field on the dot is given by Eq. (4) and takes the form where the 2 × 2 structure is the Nambu space. In all other replica except of a = 1 the Q S -matrix has the form Eq. (5). One may check that the Usadel equation is solved by exactly the same O(2) rotation as in W = 0 case, Eq. (18), performed in each of the following 2 × 2 Nambu blocks the Nambu indexes are explicitly stated on the l.h.s. In this 2 × 2 block Λ N = σ 3 and therefore the matrix from the r.h.s of Eq. (23) may indeed be rotated into Λ N as in Eq. (18). The only exception are the two lowest Matsubara components n = 0 and n = 1, which form 4 × 4 block of the form The corresponding 4 × 4 block of the Λ N -matrix on the normal lead is the unit matrix. Obviously the unit matrix cannot be unitary rotated into the block Eq. (24) and therefore solution of the Usadel equation in this block is impossible. The difficulty originates from the fact that, due to the random phase of the CB oscillations, the average free energy, F (q) , is not an oscillatory function. We need therefore to consider winding number configuration of the form W 1 = 2 and W p+1 = −2, while all other W a = 0. The contribution to the correlation function from such a configuration is proportional to exp{2πi(q − q ′ )}, that is the lowest normal harmonic of the correlation function F (q)F (q ′ ) . The "dangerous" 4 × 4 block in the p + 1-st replica is given by while the corresponding 4 × 4 block of the Λ N -matrix is minus one times the unit matrix. Here as well the unitary rotation between the points r = L and r = 0 is impossible. However, if one combines a = 1 and a = p + 1 replica and allows rotation between them -then the unitary rotation may be found [15]. Indeed, combining both "dangerous" blocks into the single 8 × 8 block (i.e. combine Eq. (24) with Eq. (25) on the dot and unit matrix with the minus unit matrix on the normal lead), one readily see that they may be unitary connected (since they possess the same set of eigenvalues).
The calculations are simplified in the low-temperature case, T ≪ ∆, where θ n ≈ π/2. In this case the 8 × 8 block takes the form on the dot and on the normal lead (σ 0 is the unit matrix in the Nambu space). Here the outer 2 × 2 unit matrix represents the space of n = 0 and n = 1 Matsubara components, while the inner one -the replica space of a = 1 and a = p + 1 (Pauli matrices act in the Nambu space). We seek thus for the solution of the Usadel equation in the 8 × 8 sub-space, having in mind that the solution for all other Matsubara components and replica is exactly the same as in W = 0 sector, Eqs. (18)- (21). The general solution for r = 0 may be written as where u(r) function was defined after Eq. (18) and B is a coordinate-independent Nambu matrix. Employing the singular value decomposition, it may be written as where U = e iu3σ3 e iu2σ2 e iu1σ1 and V = e iv3σ3 e iv2σ2 e iv1σ1 are SU (2) matrices and Θ 1 ≤ Θ 2 are real singular values. Substituting the solution back into the action, one obtains (for each of the two involved Matsubara frequencies) where the subscript ±2 stays for the corresponding winding numbers and the coefficient two on the l.h.s. reminds that two replica were involved. The next step is to minimize the action over Θ 1,2 , u 2 and v 2 angles. Three of the four equations for the minima have only trivial (parameter independent) solutions: Θ 1 = 0; u 2 = −π/4 and v 2 = π/4. The action in terms of the single non-trivial angle Θ 2 finally takes the form where, as above, t ≡ G T /G D . The corresponding saddle point equation for Θ 2 = Θ 2 (t) is For t ≤ 1 the only solution of this equation is Θ 2 = 0, while for t ≥ 1 the angle Θ 2 (t) interpolates between zero (for t = 1) and π (for t ≫ 1).

C. CB suppression
We are now on the position to discuss the suppression of the CB. Consider first the component of the correlation function Z p (q)Z p (q ′ ) which is proportional to cos 2π(q −q ′ ). As was explained above the relevant field configurations are those having a single replica with W a = ±2, where a ∈ [1, p] and a single replica with W a = ∓2, where a ∈ [p+1, 2p]. The corresponding action is given by 4S ≡ 2S ±2 + (2p − 2)2S 0 ; (here S 0 is multiplied by the number of replica with zero winding number -2p−2, and by factor of two, because two Matsubara components are different between W = ±2 and W = 0). Taking the replica limit p → 0 and employing Eqs. (20), (31), one finds where θ(t) and Θ 2 (t) are the solutions of Eqs. (21) and (32) correspondingly. As a result, the contribution to the correlation function with the unit period has the form: exp{−4S} cos 2π(q − q ′ ), where the factor of 4 stays for the fact that two replica are involved in the correlation function and the winding number is two. We shall discuss this result in more details in section VI, after deriving it using other methods. For the later reference we need to calculate t 2 ∂ t (2S/t); employing the saddle point equations (21) and (32), we obtain As was explained in section II, the parity effects sets in at small temperature T < T * ≪ ∆. Technically it manifests itself in the appearance of the odd winding numbers. The calculations for the field configuration with W 1 = 1, and W p+1 = −1 are exactly parallel to the one for the even winding numbers. The only difference is that there is the single "dangerous" Matsubara frequency, n = 0, consequently there is no outer 2 × 2 Matsubara structure as in Eqs. (26) and (27) (the replica and Nambu structures are exactly the same). As a result, the action is exactly twice smaller than for the oscillations with the unit period. We thus obtain for the component of the correlation function with doubled gate voltage period: exp{−2(S + A(T ))} cos π(q − q ′ ), where the factor of two in the exponent originates from the two replica involved.
Finally, let us discuss the technicalities of the replica limit, p → 0 [15]. There are altogether p 2 distinct configurations of the winding numbers which contribute to the correlation function Z p (q)Z p (q ′ ) with the same oscillatory factor, say exp{2πi(q − q ′ )}. They correspond to the p possible choices for W a = 2 and p independent possible choices of W a = −2. There is thus p 2 combinatorial factor in the Z p (q)Z p (q ′ ) that is cancelled against the same factor in the denominator for F (q)F (q ′ ) = lim p→0 p −2 ( Z p (q)Z p (q ′ ) − 1) .

IV. PROXIMITY ACTION
In this Section we discuss the derivation and simple applications of the proximity action approach to treat the interface between the SC dot and normal lead. This method was recently suggested [11,25] for the analysis of superconductive-normal metal transition in a 2D proximity-coupled array [26,27]. Notice that our definition of the running parameter ζ and charges {γ l } differ from those used in Refs. [11,25,27], see Appendix A for details.

A. Derivation of the proximity action
Our immediate goal is to derive the proximity action, Eq. (16), starting from Eqs. (13) and (14). To this end, let us imagine adding the diffusive metal step by step with infinitesimal shells having conductance δG ≫ G D > 1. At some intermediate step of this procedure one has the tunneling barrier with the part of the normal metal incorporated into the set of coefficients γ l (ζ). Here ζ ≡ G D /G D (ζ) ∈ [0, 1] is the running parameter such that γ l (0) = δ l,1 G T /G D , and S(0) = S T is the bare tunneling barrier action; while γ l (1) = γ l and S(1) = S T D is the full proximity action of the barrier "dressed" with the diffusive metal. At each step one adds another thin shell of the diffusive metal. The action of the entire system takes the form where Q is the field right on the boundary between the new shell and the already integrated region. The action of the newly added shell has the simple form, (δG/16) Tr(Q − Λ N ) 2 , since δG ≫ G D > 1 and therefore the matrix Q has to be rather close to Λ N (in the leading order in (δG) −1 one may disregard all higher order terms (Q − Λ N ) l with l > 2). The next step is to integrate out the Q field and obtain the new set γ l (ζ + δζ). We perform this procedure in the saddle point approximation. Taking variation of the action, Eq. (35), under the condition Q 2 = 1 one finds the saddle point (Usadel) equation In the leading order in G D /δG ≪ 1 one obtains for the saddle point solution Putting it back into the action, one finally finds where we employed the fact that δG = G D /δζ. One ends up with the following evolution equations for γ l (ζ): To solve this set of equations it is convenient to define the function u(ζ, Θ) ≡ supplemented with the initial condition u(0, Θ) = t sin Θ, where, as before, t ≡ G T /G D . Employing the method of characteristics, one may write the solution of the latter problem in the implicit form: As a result, the function u 0 (Θ) ≡ u(1, Θ), that has lγ l as it's Fourier coefficients, may be found as the solution of the following differential or algebraic equations The last expression solves the problem of writing the action of the barrier and the diffusive region in the form of Eq. (16). Note, that to derive this result we have not used a specific form of Q S and Λ N (other than the fact that Q 2 S = Λ 2 N = 1). Therefore the action and the coefficients γ l are equally applicable to any other setup. The only two conditions that were employed are G D ≫ 1 to disregard localization effects and perform the saddle point calculations and E Th ≫ E C to disregard the frequency term in the action and justify the structure of Eq. (16). The condition E Th ≫ E C also allows to neglect Cooper-channel interaction in the normal conductor, while deriving proximity action (cf. [11] for the discussion of general situation).

B. Physical observables
In this subsection we repeat briefly derivation presented originally in [11]. Normal state conductance. Consider for a moment the dot in the normal state with a small applied bias V (τ ) = V cos ω m τ . The dot's Green function takes the form Substituting such Q S in the action and taking the second variation with respect to V (τ ) at V = 0, one finds for the linear normal state conductance For small Θ one has u 0 (Θ) ≈ Θ∂ Θ u 0 (0); employing Eq. (42), one finds ∂ Θ u 0 (0) = t/(1 + t). As expected, one obtains the resistance addition rule Andreev conductance. We shall concentrate now on the temperature region ∆ ≫ T ≫ δ, where the SC Q-matrix may be written in the form Eq. (4) with θ n = π/2, that corresponds to infinite ∆ (the leading effect of the finite ∆ is to renormalize E −1 C → E −1 C + G/2π∆ [28,24]). In this case the Q S matrix is off-diagonal in the Nambu space and time-local. Since the Λ N -matrix is diagonal in the Nambu space, only even powers of (Q S Λ N ) may contribute to the action. As a result, one finds for the proximity action is the Green function of the normal lead (the Fourier transform of Λ N (ǫ n ) to the time domain is Λ N (τ, τ ′ ) = G(τ − τ ′ )σ 3 ). If the external voltage V (τ ) is applied to the SC dot, the phase is to be understood as Φ = dτ V (τ ). Calculating the second variation of the above action with respect to V (τ ), one finds for the linear (Andreev) conductance Employing Eq. (42) for Θ ≈ π/2 one finds where θ = t cos θ = u 0 (π/2). One therefore obtains

C. Instanton treatment of the proximity action
To calculate the proximity action on the simplest instanton trajectory Φ 1 (τ ) = πT τ W and Φ a (τ ) = 0 for a ∈ [2, 2p], one may employ e.g. the Matsubara basis (see Appendix A for calculations in the time domain). Employing Eq.
where the last equality is a direct consequence of the definition u 0 (Θ) = ∞ l=1 lγ l sin(lΘ). Employing the differential equation (42), one finds Comparing with Eq. (34), we find an exact coincidence of the real-space result, Eq. (33), and the proximity action, Eq. (49), S(t) = S 1 (t). Notice that, since the (random) diffusive region was integrated out upon derivation of the proximity action, the phase of the CB oscillations is not random. In a sense the proximity action represents a typical diffusive region (mesoscopic fluctuations are omitted in the derivation given above). One may therefore calculate directly typical Z(q), without resorting to the correlation function. From the two smallest winding numbers W = 1 and W = 2 one finds Z(q) ≈ exp{−(S 1 + A(T ))} cos πq + exp{−S 2 } cos 2πq + . . ., in agreement with the result of the previous section.
It is worth mentioning here that the general expression for the instanton action in the case of the normal island, derived in Ref. [14], may be presented in the form similar to Eqs. (33) and (49): and only even W are allowed for the normal case. It is easy to see, using Eq. (32), that in the entire interval G D > G T the value of the action S N 2 is the same as in the tunnelling limit: S N 2 = G T /2. As t = G T /G D approaches unity, a kind of phase transition occurs [29], which reveals itself as a non-analytic behavior of S N 2 (t) at t = 1.

D. Landauer approach
Employing Nazarov's techniques [29], one may show that the function u 0 (Θ) is directly related to the generating function of the transmission coefficients where T α are eigenvalues of the N × N matrix t † t describing transmission between the dot and the lead. Factor 1/2 on the l.h.s. takes into account that G D is defined for the two spin components. Performing the Fourier transform one finds where T α = cosh −2 µ α . As a result, the proximity action takes the compact form: In the last expression we have omitted a constant which goes to zero in the replica limit. Exactly the same interface action is known in the supersymmetric formalism [30,6] (notice that here we deal with spin 1/2 electrons). One may calculate now the action on the instanton trajectory to find the renormalized charging energy. Employing Eq. (53) and recalling that S 1 = −(G D /2) k γ 4k−2 , one finds This expression is to be compared with the one for the Andreev conductance obtained from Eqs. (46) and (53) [17,18]: (c.f. with the normal conductance G N = G D l l 2 γ l = 2 α T α , where 2 stays for spin). Factor 4 on the r.h.s. of Eq. (56) stems from the fact that Andreev particles carry charge 2e and no spin. One may thus identify the combination T α = T 2 α /(2 − T α ) 2 as the Andreev conductance in channel α. The renormalisation of the charging energy is therefore given by This may be compared with the corresponding factor for the normal spinless particles [14]: α (1 − T α ) 1/2 . Indeed, in our notations the normal instanton action is S N 2 = (G D /2) l γ 2l−1 = − α ln(1 − T α ); for spinless particles one has exp{−S N 2 /2} for the CB suppression. There is thus a perfect analogy between Andreev and normal spinless charge fluctuations mechanisms, provided that the normal transmission coefficients, T α , are substituted by the Andreev ones, T α . From Eq. (57) one may extract some general results. In the tunnelling limit: T α ≪ 1 in all channels, one finds: S 1 = G A /8. In the diffusive metal interface the transmission eigenvalues are distributed according to the Dorokhov distribution [16,29]: P (T ) = G D /(4T √ 1 − T ). The corresponding distribution for the Andreev transmissions is given by P (T ) = G D /(8T √ 1 − T ), as a result G A = G D . For the typical action of the diffusive interface one finds therefore S 1 = π 2 G A /32.

E. Fluctuation determinant and summation over many-instanton configurations
We now expand the action, Eq. (45), to the second order in deviations from the instanton trajectory Φ a (τ ) = πT τ W a + δΦ a (τ ), where δΦ a (0) = δΦ a (β). After diagonalization of the resulting quadratic form one finds the following spectrum of the small fluctuations [12,13,31]: where n ∈ [−∞, ∞]. For T ≪ E C there are 2|W a | + 1 (almost) zero modes labelled by n ≤ |W a |. One of them n = 0 is the trivial shift of Φ a (τ ) by a constant, whereas the remaining 2|W a | zero modes are associated with the deformation of the instantons. The general solution of the saddle point equations (for T ≪ E C ) known as the Korshunov instanton [32] may be written as where z k is a set of |W a | complex numbers (instanton coordinates) which parametrize the 2|W a | dimensional zeromode manifold. The instanton coordinates, z k , are complex numbers from inside of the unit circle. One may thus characterize any deviation δΦ a (τ ) by 2|W a | zero-mode coordinates, z k , and remaining transversal modes having nonzero masses, Eq. (58). The Jacobian of such transformation may be found in the standard way [33] and is given by where l, k ∈ [1, |W a |]. The fact that the Jacobian vanishes if two of the coordinates coincide signals the repulsive interaction between the instantons. The Gaussian integration over the massive degrees of freedom normalized by those in the zero winding number sector results in up to the terms of the order T /(G A E C ) ≪ 1. Combining together all the factors, one finds for the one instanton, W a = 1, trajectory The cutoff of the logarithmically divergent z-integral originates from the following consideration. In presence of the charging energy term the Korshunov instantons, Eq. (59), are not true zero-modes of the action. The z-dependence of the action brings the factor exp{−(T /2E C )|z| 2 /(1 − |z| 2 )} to the integral. As a result the effective integration range shrinks to |z| 2 < 1 − T /E C < 1. The renormalized charging energyẼ C is defined as It is also instructive to look at the W a = 2 contribution to the partition function The leading power of the large logarithm comes from the diagonal term of the Jacobian, Eq. (60). This term corresponds to the non-interacting instantons approximation, which is justified therefore by the large parameter ln(E C /T ) ≫ 1. Notice, that the (repulsive) interaction correction (the off-diagonal part of the Jacobian, Eq. (60)) comes with ln 0 and not ln 1 , as may be naively expected. Due to the large logarithmic factor, originating from the integration over sizes of each individual instanton (i.e. neglecting their interactions), one may treat positive and negative instantons (anti-instantons) on the equal footing. Strictly speaking, a combination of instantons and anti-instantons is not a saddle point solution. However, relative weakness of instanton interactions (no log-factor in the second term in Eq. (64)) makes it possible to consider an ideal gas of instantons and anti-instantons. The contribution to the partition function with a given winding number W is given by all configurations having m + |W | instantons and m anti-instantons in an arbitrary order. Taking only the terms with the leading power of the logarithms (diagonal term of the Jacobian) one finds [13] This ugly looking series may be summed up into an unexpectedly simple expression: The simplest way to check it is to expand back Eq. (66) in a double series in powers ofẼ C and e iπq . The role of the anti-instantons, therefore, is to convert cos(W πq), which may be expected from the instantons only, into cos W (πq).
Under the conditionẼ C ≤ T ≪ E C the interaction between instantons can be treated perturbatively; for the lowest-order correction we use the last term in Eq. (64) to find for the gate voltage -dependent free energy We retain the second term in (67) since it will be important at temperatures comparable with the parity-effect temperature T * , see Sec. VI. It also demonstrates that the instanton-instanton interaction corrections are small only ifẼ C < T .

V. RG TREATMENT OF THE QUANTUM FLUCTUATIONS
We next approach the problem of the CB from the different perspective. Instead of calculating the action on the instanton field configurations, we shall perform the RG analysis of the coefficients γ l upon integrating out fast fluctuations of Φ(τ ). To this end we write Φ(τ ) = Φ s (τ ) + Φ f (τ ), where superscripts s and f stay for slow (ω m < Ω) and fast (ω m > Ω) component of the field correspondingly. The running cutoff Ω runs from Ω ∼ G A E C down to the lowest Matsubara frequency Ω ∼ T . We next substitute the Φ-field into the SC proximity action, Eq. (45), expand to the second order in Φ f and integrate out the fast field fluctuations. The bare propagator of the Φ f -fields may be read out from the action Eq. (45): where G A (ζ Ω ) = G D ∞ l=1 (−1) l (2l) 2 γ 2l (ζ Ω ) is the running value of the Andreev conductance, and ζ Ω ≡ ln(G A E C )/Ω. (Note that this definition of ζ differs from that used to derive the proximity action in Section IV.) The further calculation is essentially similar to that in the Keldysh formalism [27]. One obtains the following set of the evolution equations for γ l (ζ Ω ): where we have omitted the subscript Ω. Since only the even coefficients are involved it is convenient to define the functionũ(ζ, Θ) = [u(ζ, Θ) − u(ζ, π − Θ)]/2 which contains only even harmonics. In terms of this function the RG equations (69) take the formũ where the Andreev conductance is G A (ζ) = G DũΘ (ζ, π/2). The initial condition for Eq. (70) isũ(0, Θ) =ũ 0 (Θ), with the u 0 (Θ) function evaluated above. Let us mention for completeness, that in the normal case the corresponding RG equation takes the form u ζ = (2/G N (ζ))[u(Θ) cot(Θ/2)] Θ , where the running normal conductance according to Eq. (43) is given by 0). The key parameter that determines the strength of phase fluctuations is the Andreev conductance G A (ζ). We consider the large bare value, G A (0) ≫ 1. Upon integration over fluctuations of Φ(τ ), the effective value of G A (ζ) decreases (the physical reason of this phenomena is a partial loss of the phase coherence between multiple Andreev reflections), and becomes comparable to unity at some time scale Ω −1 c = (G A E C ) −1 e ζc . At longer times coupling between the dot and the lead is weak, thus the phase fluctuations grow and autocorrelation function C(τ ) = e i(Φ(τ )−Φ(0)) decays fast. It is natural to associate Ω c with an effective Coulomb energy of the island connected to the wire,Ẽ C . We will show now that the estimate for Ω c that follows from the RG equation coincides (within exponential accuracy) with the results of the instanton analysis.
Equation (70) may be solved with the method of characteristics. We define the function v(Θ) =ũ(Θ) tan Θ and then introduce new "space" and "time" variables: ξ = − ln sin Θ and τ which is defined via the relation In terms of these new variables the RG equation (70) takes the simple form: v τ = v ξ , yielding v(ξ, τ ) = v 0 (ξ + τ ). As a result,ũ (ζ(τ ), Θ) = e −τ cos Θ 1 − e −2τ sin 2 Θũ 0 (arcsin(e −τ sin Θ)) , The dependence ζ(τ ) can be now found from Eq. (71) with the Andreev conductance given by G A (ζ) = G DũΘ (ζ, π/2). Our goal is to find a scale ζ c , where the Andreev conductance G A becomes small. To find ζ c we notice that in terms of the RG time τ the smallness of G A means the limit τ → ∞. Then, using definition of τ given by Eq. (71), we find Recalling the definition ofũ 0 (Θ) and comparing this result with Eq. (49) we conclude that ζ c = S 1 . As a result, the amplitude of the lowest oscillatory component (∼ cos πq) of the free energy, from the RG point of view, is given by the precisely the same exponential factor as in the framework of the instanton calculation. Yet another instructive way to approach the problem is to follow the renormalization of the instanton action.
where in the last equality we have used that G A (ζ) = G DũΘ (ζ, π/2). (In the normal case one finds dS N W /dζ = −|W |/2; only even W are allowed.) From the dimensional analysis one expects F (q) ∼ Ω exp{−S 1 } cos πq + . . ., where Ω ∼ E C G A is the cutoff energy. We may now renormalize down the cutoff energy Ω → Ω(ζ) ≡ Ωe −ζ , simultaneously following renormalization of the action, S 1 → S 1 (ζ). According to Eq. (74), S 1 (ζ) = S 1 − ζ. As a result, Ω(ζ) e −S1(ζ) = const . (75) This means that the lowest oscillatory component of the free energy is expected to be temperature independent, in agreement with Eq. (67) (up to log-factor, originating from the zero modes). The instanton calculation is done at ζ = 0, alternatively the RG calculation presented above looks for ζ c , such that S W (ζ c ) ≈ 0 (the scale where the instantons do not cost anything). In view of Eq. (75) it is not surprising that they yield the same result. Moreover, one can renormalize down to any intermediate scale ζ (such that G A (ζ) > 1) and then calculate the instanton action with the current parameters {γ l (ζ)} -the result is still guaranteed to be the same. Notice also that the renormalized action S W (ζ) = S W − |W | ln(G A E C /T ) coincides with the result that comes from the instanton action together with the massive Gaussian fluctuations, c.f. Eqs. (61)-(64). The effect of the zero modes (ln E C /T in the prefactor) is missed in the one-loop RG calculation. Finally we comment upon the comparison of the approach used in this Section with Ref. [27] where the same kind of RG approach was employed in the case of 2-dimensional normal conductor with very low Thouless energy, E Th ≪ E C . In the last case the RG procedures involves simultaneous integration over phase fluctuations Φ(τ ) and diffuson/Cooperon modes in the normal conductor. As a result, the full functional RG equation contains three terms in the r.h.s.: one coming from Eq. (40), the second is similar to Eq. (70), and the third contribution accounts for the Cooper-channel repulsion in the normal conductor (cf. [11,25]). The resultant RG equation cannot be reduced to the set of even harmonicsũ only; we are not aware of any method to solve it apart from "brut-force" numerical integration. Considerable simplification of the problem treated here stems from the fact that, due to the condition E Th ≫ E C all diffuson/Cooperon modes may be integrated out before (means at larger energies than) the phase fluctuations Φ(τ ) become relevant. This procedure leads to the action functional S T D that depends upon Φ(τ ) trajectory only, as defined in (45). Thus the results presented in this Section can be considered as "minimal generalization" of those obtained in Ref. [26].

VI. DISCUSSION OF THE RESULTS
The main result of our study, given by Eq. (67), applies in the temperature rangẽ where T * is the parity-effect temperature given by Eq. (10). The condition T ≪ T * ensures that both even and odd winding numbers contribute to the partition function, since A(T ) ≪ 1, c.f. Eq. (9). The two lowest oscillatory components of the free energy are given then by Eq. (67). At higher temperature T ∼ T * an addition of odd number of electrons to the dot becomes possible, changing the relative amplitude of the harmonics. On the level of the partition function the components with odd winding numbers acquire the temperature dependent factor exp{−A(T )}. As a result, one finds for the gate voltage dependent free energy: Finally at T > T * the parity effect disappears and only the normal (unit) oscillation period remains in the free energy 4T cos(2πq) + . . .
Notice that if ln E C /T * > π/ √ 6 there is the sign change of the cos(2πq) component at T ≈ T * . Moreover, unlike the low temperature case, where the amplitude of the dominant harmonics was weakly (logarithmically) temperature dependent, Eq. (67), at larger temperature there is the stronger dependence (∼ T −1 ln 2 (E C /T )) of the amplitude. The characteristic crossover temperature is determined from the relation A(T † ) ≈ ln(T † /Ẽ C ) and is given by where we have used the fact that S 1 ∼ G A .
We emphasize an interesting feature of the weak Coulomb blockade compared to the usual one: although the role of effective Coulomb energy is taken by the effective energyẼ C ≪ E C , the oscillation amplitude depends relatively weakly upon temperature within the range given by Eq. (76). This is to be contrasted with the usual case of an island connected by highly resistive contacts: at temperatures T > E C oscillation amplitude vanishes exponentially fast with T /4E C .
Next we shall discuss the quantitative value of the renormalized charging energyẼ C , which determines both the magnitude of the CB and the region of applicability of our results. To this end we need to evaluate S 1 (t) (where t = G T /G D ), which is given by one of the two equivalent expressions, Eqs. (33) or (49). In the two limiting cases t ≪ 1 (the tunneling barrier limit) and t ≫ 1 (the diffusive metal limit) the answer may be found analytically. One obtains in the tunneling limit, π 2 /4 ; in the diffusive limit, where the Andreev conductance is given by Eq. (47). Recall that for a normal dot the corresponding exponent is given by G N /2 in the tunnelling limit and G N π 2 /8 in the diffusive limit. We observe that upon the same dissipative conductance the SC dot exhibits factor of four (in the two limiting cases) smaller action. Therefore one may observe a sizeable CB in the superconductive state, while in the analogous normal dot the CB is practically suppressed. The crossover behavior for the SC action S 1 (t), the normal action S N 2 (t) and the Andreev conductance G A (t) as function of t = G T /G D is depicted in Fig. 2. Dependence of the same quantities on the resistance G −1 D of the normal region at fixed G T is demonstrated in Fig. 3.
In the caseẼ C ≪ T † the crossover between T < T † and T > T † regimes leads to the sharp drop of the oscillation amplitude by the factorẼ C 4T † ln EC T † ≪ 1 The same drop in the residual Coulomb blockade may occur as function of magnetic field, due to the suppression of the ∆ and thus of T † by the magnetic field. This may result in a strong negative magnetoresistance at low temperatures of a granular media made out of small superconductive grains (cf. e.g. [34] for the discussion of relevant experiments). This effect can occur if the low-temperature oscillation amplitudẽ E C is sufficiently large to destroy Josephson coupling between grains, so transport of Cooper pairs between grains is blocked byẼ C . However, quantitative theory of such an effect is still to be developed.
Another open question concerns the behavior of F (q) in the zero temperature limit, T <Ẽ C . As temperature decreases belowẼ C , all approximate methods, we used, run out of their applicability range: the renormalized conduc- tance G A drops below unity, fluctuations become strong and there is no apriori reason to treat them within Gaussian approximation. Whereas the overall scale of the oscillation amplitude is probably given by its first harmonic, Eq. (67) at T ∼Ẽ C , and is of the orderẼ C ln EC EC , the precise shape of the oscillations is still to be determined. It is possible that in the T → 0 limit the function F (q) becomes nonanalytic at the degeneracy point q = 1/2. An extreme case of such a nonanalytic behavior -finite steps in dF/dq at half-integer q -was found in Ref. [24] where the Andreev conductance was completely neglected. To proceed we need to establish a correspondence between our notations and those adopted in Refs. [35,11,25]. The latter notations will be designated by a prime. Firstly, our conductance quantum G Q = e 2 /2πh is different from G ′ Q = e 2 /h used in those papers. Secondly, the running parameter ζ and charges {γ n } are related by Thus, the Keldysh multicharge proximity action may be written in our notations as [36] where the trace is taken over time indexes, Nambu and Keldysh spaces. HereQ S is a matrix in the superconducting dot acting in the Nambu space asQ where Φ ↔ = diag(Φ > , Φ < ) is a matrix in the Keldysh space, with Φ > (t) and Φ < (t) being the fields residing on the forward and backward branches of the Keldysh contour. In the normal lead,Λ N (E) =Λ 0 (E) ⊗ σ 3 , where σ 3 is the Pauli matrix in the Nambu space, and the matrixΛ N (E) acts in the Keldysh space: f (E) being the distribution function. Tracing over the Nambu space reduces the action (A2) to the form: where the trace is taken only over the Keldysh space. Transition to the imaginary time is achieved by the deformation of the Keldysh contour C which initially run over the time axis from −∞ to ∞ and then in the backward direction. The desired deformation introduces a vertical segment of length β at some time t = t 0 so that the contour originates at −∞ + iβ/2, runs through the points t 0 + iβ/2, t 0 to ∞ and then back through t 0 , t 0 − iβ/2 to −∞ − iβ/2, see Fig. 4. For the purpose of the evaluation of the instanton action the choice of the point t 0 is somewhat arbitrary. Note, however, that were we to consider, e.g., the correlation function of the form Ñ (q)Ñ (q ′ ) , it would be necessary to allow for adiabatic evolution of the gate voltage q = q(t), q ′ = q(t ′ ) and to introduce two vertical segments at times t and t ′ .
The analysis is simplified in the zero-temperature case. Choosing the position of the vertical segment at t 0 = 0, we deform the forward segment (−∞, 0) of the Keldysh contour to the upper imaginary half-axis and the backward segment (0, −∞) to the lower imaginary half-axis. The remaining horizontal appendix of the Keldysh contour running from 0 to ∞ and then back to t = 0 can be neglected as it does not contribute to thermodynamic quantities.