text:
# Coulomb blockade in metallic grains at large conductance

## I Introduction

The study of electron-electron interactions in mesoscopic systems has been at the focus of experimental and theoretical interest over the past two decades. One of the most striking manifestations of electron interactions at low temperatures is the phenomenon of Coulomb blockade.  It can be observed by measuring, say the charge of a metallic grain which is connected by a tunneling contact to a metallic lead and is capacitively coupled to a metallic gate that is maintained at voltage $`V_g`$, as in Fig. 1.

For a poorly conducting contact the charge in the grain is quantized at low temperatures and exhibits a characteristic step-like dependence on the gate voltage $`V_g`$. This leads to oscillatory gate voltage dependence of all physical quantities, referred to as the Coulomb blockade oscillations. The amplitude of these oscillations decreases with increasing transparency of the contact and depends not only on the total conductance of the contact but on the individual transparencies of the tunneling channels.

In the case of a few transmission channels in the contact and in the limit of vanishing mean level spacing, $`\delta 0`$, the Coulomb blockade oscillations of the thermodynamic, transport and thermoelectric quantities are completely smeared at perfect transmission.  For a finite mean level spacing in the grain, the Coulomb blockade oscillations remain finite even at perfect transmission and exhibit strong mesoscopic fluctuations.

In the case of a contact with many weakly transmitting tunneling channels, the Coulomb blockade oscillations at large conductance have been studied in the limit of vanishing mean level spacing in the grain, $`\delta `$, and shown to be exponentially suppressed  $`\mathrm{exp}(g_T/2)`$. Here $`g_T`$ is the conductance of the tunneling contact measured in units of $`e^2/(2\pi \mathrm{})`$. In the $`\delta 0`$ limit for a non-random contact the mesoscopic fluctuations in the Coulomb blockade oscillations vanish. At finite mean level spacing the mesoscopic fluctuations in the thermodynamic quantities for a weakly transmitting tunneling contact were studied in Ref.  and for a multi-channel contact near perfect transmission in Ref. . In the case of a random diffusive contact and vanishing mean level spacing in the grain the Coulomb blockade oscillations were studied in Refs. .

Here we study the Coulomb blockade oscillations in the thermodynamic quantities of a disordered metallic grain with finite mean level spacing $`\delta `$ which is connected to a metallic lead by a non-random multi-channel tunneling contact and is capacitively coupled to a gate, as in Fig. 1. The conductance of the tunneling contact is assumed to be large, $`g_T1`$. We assume that electron-electron interactions in the metallic lead may be neglected due to screening and treat the Coulomb interaction of electrons in the grain in the framework of the constant interaction model:

$$\widehat{H}_C=E_C\left(\widehat{N}q\right)^2.$$
(1)
Here $`E_C`$ is the charging energy, $`\widehat{N}`$ is the operator of the number of electrons in the grain, and $`q`$ is the dimensionless parameter which is proportional to the gate voltage $`V_g`$ and has the meaning the number of electrons that minimizes the electrostatic energy of the grain. We consider a disordered metallic grain in which electrons move in the presence of a random impurity potential and study both ensemble averaged thermodynamics quantities of the grain and their mesoscopic fluctuations. We assume that the Thouless energy of the grain, $`E_T=D/L^2`$, where $`D`$ is the diffusion coefficient and $`L`$ is the grain size, is greater than the charging energy, $`E_T>>E_C`$. In addition, we assume that the time reversal symmetry is broken by a magnetic field $`H`$ and that the cooperon gap $`DeH/(\mathrm{}c)`$ is greater than the charging energy $`E_C`$. Under these conditions the ensemble of disorder potentials is equivalent to the unitary random matrix ensemble.

We use the replica formalism to treat disorder averaging in the presence of interactions.  The application of the Itzykson-Zuber  integral enables us to evaluate the resulting functional integral over the $`Q`$-matrix exactly. The method used here may have useful applications to the treatment of non-perturbative interaction effects in granulated disordered systems.

The main results of the paper are the analytic expressions for the average oscillatory part of the grand canonical potential, Eq. (66) and for the oscillatory part of the correlator of the grand canonical potentials at different values of the gate voltage, Eq. (70).

The paper is organized as follows: In section II we introduce the formalism by considering a closed disordered metallic grain with non-interacting electrons in the zero dimensional limit, and obtain the ensemble average and the fluctuations of the thermodynamic potential $`\mathrm{\Omega }(q,T)`$. In section III we apply this method to compute the average and the mesoscopic fluctuations of the amplitude of the Coulomb blockade oscillations in the thermodynamic potential $`\mathrm{\Omega }(q,T)`$. We discuss our results in section IV.

## II Thermodynamics of non-interacting electrons

First, to illustrate the formalism we consider the thermodynamic quantities of an isolated disordered metallic grain. In this section we neglect the electron-electron interactions in the grain. For simplicity, in this section we consider spinless electrons.

We are working within the grand canonical ensemble with the chemical potential $`\mu `$. All thermodynamic quantities can be obtained from the grand canonical potential $`\mathrm{\Omega }_0(\mu ,T)=T\mathrm{ln}Z_0(\mu ,T)`$. Here the subscript $`0`$ indicates that the quantities pertain to non-interacting electrons and $`Z_0(\mu ,T)`$ denotes the partition function which can be written as a functional integral over the fermionic variables $`\psi `$, $`\overline{\psi }`$ as

$$Z_0(\mu ,T)=\underset{x=1}{\overset{N}{}}\underset{n}{}d\overline{\psi }_{x,n}d\psi _{x,n}\mathrm{exp}\left[\underset{n,m}{}\underset{x,y=1}{\overset{N}{}}\overline{\psi }_{x,n}\left(i\widehat{\epsilon }\delta _{nm}(\widehat{H}_{0,xy}\mu \delta _{xy})\right)\psi _{y,m}\right].$$
(2)
In this equation the fermion variables are labeled by the Hilbert space indices $`x`$ and $`y`$ and by the Matsubara indices $`n`$ and $`m`$, the diagonal matrix $`\widehat{\epsilon }=\delta _{xy}\delta _{nm}\epsilon _n`$ is expressed through the fermionic Matsubara frequencies $`\epsilon _n=(2n+1)\pi T`$, $`\widehat{H}_0`$ is the single particle Hamiltonian, and $`\mu `$ is the chemical potential.

We consider the case of broken time-reversal symmetry. Therefore, in the zero dimensional limit  the single particle Hamiltonian $`\widehat{H}_0`$ can be modeled by a random $`N\times N`$ Hermitian matrix belonging to the Gaussian Unitary Ensemble (GUE) with the probability distribution

$$P(H_0)\mathrm{exp}\left(\frac{N}{2}\mathrm{Tr}\widehat{\mathrm{H}}_0^2\right).$$
(3)
As is well known , the average single particle density of states (DOS) $`\nu (E)`$ for such a distribution is given by the semi-circle law

$$\nu (E)=\frac{N}{\pi }\sqrt{1\frac{E^2}{4}}.$$
(4)
Note we normalized the distribution (3) in such a way that the width of the semi-circle is equal to unity. Thus, throughout this paper all energies are measured in units of the width of the semi-circle. For simplicity, we assume that the chemical potential $`\mu `$ lies near the middle of the semi-circle, where the mean level spacing in the dot is given by

$$\delta =\frac{\pi }{N}.$$
(5)
Below we will consider the ensemble averaged thermodynamic properties of the dot. For this purpose we resort to the replica trick  and find the averaged replicated partition function

$$Z_0^\alpha (\mu )=\mathrm{exp}\left(\alpha \frac{\mathrm{\Omega }_0(\mu )}{T}\right),$$
(6)
where $`\mathrm{}`$ denotes the averaging over the ensemble of random matrices $`\widehat{H}_0`$ with the probability distribution function defined in Eq. (3), and $`\alpha `$ is the number of replicas. The function $`Z_0^\alpha (\mu )`$ is initially determined for the positive integer number of replicas $`\alpha `$ and then analytically continued to $`\alpha 0`$. The $`n`$-th cumulant of the thermodynamic potential $`\mathrm{\Omega }_0(\mu )^n`$ with respect to the distribution Eq. (3) is then found from $`Z_0^\alpha (\mu )`$, Eq. (6) by using the formula

$$\mathrm{\Omega }_0(\mu )^n=T^n\frac{d^n\mathrm{ln}Z_0^\alpha (\mu )}{(1)^nd\alpha ^n}|_{\alpha =0}.$$
(7)
We express the replicated partition function $`Z_0^\alpha (\mu )`$ through the functional integral over the replicated fermionic fields $`\psi ^j`$ as in Eq. (2) and average $`Z_0^\alpha (\mu )`$ with respect to the probability distribution function Eq. (3). As a result we obtain

$`Z_0^\alpha (\mu )`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{\alpha }{}}\underset{x=1}{\overset{N}{}}\underset{n}{}d\overline{\psi }_{x,n}^jd\psi _{x,n}^j\mathrm{exp}\left[\underset{n}{}\underset{x=1}{\overset{N}{}}\underset{j=1}{\overset{\alpha }{}}\overline{\psi }_{x,n}^j\left(i\epsilon _n+\mu \right)\psi _{x,n}^j+\frac{1}{2N}\underset{x,y=1}{\overset{N}{}}\left|\underset{n}{}\underset{j=1}{\overset{\alpha }{}}\overline{\psi }_{x,n}^j\psi _{y,n}^j\right|^2\right]}.`$ (8)
Next, we introduce a Hermitian $`\alpha 2M\times \alpha 2M`$ matrix $`\widehat{Q}`$ to decouple the quartic term in Eq. (8) via the Hubbard-Stratonovich transformation:

$`\mathrm{exp}\left[{\displaystyle \frac{1}{2N}}{\displaystyle \underset{x,y=1}{\overset{N}{}}}\left|{\displaystyle \underset{n}{}}{\displaystyle \underset{j=1}{\overset{\alpha }{}}}\overline{\psi }_{x,n}^j\psi _{y,n}^j\right|^2\right]=c_{\alpha \times 2M}{\displaystyle d[\widehat{Q}]\mathrm{exp}\left\{\frac{N}{2}\mathrm{Tr}\widehat{\mathrm{Q}}^2+\mathrm{i}\underset{\mathrm{x}=1}{\overset{\mathrm{N}}{}}\underset{\mathrm{n},\mathrm{m}}{}\underset{\mathrm{i},\mathrm{j}=1}{\overset{\alpha }{}}\overline{\psi }_{\mathrm{x},\mathrm{n}}^\mathrm{i}\mathrm{Q}_{\mathrm{nm}}^{\mathrm{ij}}\psi _{\mathrm{x},\mathrm{n}}^\mathrm{j}\right\}},`$ (9)
where

$$c_{\alpha \times 2M}=\left(\frac{N}{2\pi }\right)^{(\alpha 2M)^2/2}2^{\alpha M(\alpha 2M1)},$$
(10)
and the trace in Eq. (9) is taken over both the Matsubara and the replica indices. Here $`2M`$ is the number of Matsubara frequencies in each replica. We keep it finite at this point to make the sums over the Matsubara frequencies well defined. Eventually the limit $`2M\mathrm{}`$ will be taken to obtain the final expressions. The elements of the $`\widehat{Q}`$-matrix are labeled by four indices: $`Q_{nm}^{ij}`$, two of which refer to the replica space $`i,j`$, and two others, $`n,m`$ refer to the Matsubara space. After the integration over the fermionic variables in Eq. (9) we obtain

$$Z_0^\alpha (\mu )=c_{\alpha \times 2M}d[\widehat{Q}]\mathrm{exp}\left\{\frac{N}{2}\mathrm{Tr}\widehat{\mathrm{Q}}^2+\mathrm{NTr}\mathrm{ln}\left(\mathrm{i}\widehat{\mathrm{Q}}+\mathrm{i}\widehat{\epsilon }+\mu 𝟙\right)\right\},$$
(11)
where we introduced the notation $`𝟙=\delta _{𝕚𝕛}\delta 𝕞𝕟`$. Following Guhr  we make the following change of variables in Eq. (11): $`\widehat{Q}=\widehat{Q}^{}\widehat{\epsilon }+i\mu 𝟙`$. As a result we obtain the following expression

$$Z_0^\alpha (\mu )=c_{\alpha \times 2M}d[\widehat{Q}^{}]\mathrm{exp}\left\{\frac{N}{2}\mathrm{Tr}(\widehat{\mathrm{Q}}^{}\widehat{\epsilon }+\mathrm{i}\mu 𝟙)^\mathrm{𝟚}+\mathrm{Tr}\mathrm{ln}\left(\mathrm{i}\widehat{\mathrm{Q}}^{}\right)\right\}.$$
(12)
The matrix $`\widehat{Q}^{}`$ can be diagonalized using a unitary matrix $`\widehat{U}`$

$`\widehat{Q}^{}`$ $`=`$ $`U^1\widehat{\mathrm{\Lambda }}U,`$ (13)
$`\widehat{\mathrm{\Lambda }}`$ $`=`$ $`\delta _{ij}\delta _{nm}\lambda _n^j.`$ (14)
In terms of the new variables $`\widehat{U}`$ and $`\widehat{\mathrm{\Lambda }}`$ the integration measure becomes

$$d[\widehat{Q}^{}]=\mathrm{\Delta }^2(\widehat{\mathrm{\Lambda }})d[\widehat{\mathrm{\Lambda }}]d[\widehat{U}],$$
(15)
where $`d[\widehat{\mathrm{\Lambda }}]=\underset{i,n}{}d\lambda _n^i`$ and $`d[\widehat{U}]`$ denotes the invariant integration measure in the unitary group. The Jacobian of the transformation (13) is expressed through the Vandermonde determinant $`\mathrm{\Delta }(\widehat{\mathrm{\Lambda }})^{}(\lambda _m^i\lambda _n^j)`$, where the prime indicates that the product is taken over non-coinciding pairs of indices $`\{i,m\}\{j,n\}`$. In the variables $`\widehat{\mathrm{\Lambda }}`$ and $`d[\widehat{U}]`$ Eq. (12) acquires the form

$$Z_0^\alpha (\mu )=c_{\alpha \times 2M}\mathrm{\Delta }^2(\widehat{\mathrm{\Lambda }})d[\widehat{\mathrm{\Lambda }}]d[\widehat{U}]\mathrm{exp}\left\{\frac{N}{2}\mathrm{Tr}(\widehat{\mathrm{U}}^1\widehat{\mathrm{\Lambda }}\widehat{\mathrm{U}}\widehat{\epsilon }+\mathrm{i}\mu 𝟙)^\mathrm{𝟚}+\mathrm{Tr}\mathrm{ln}\left(\mathrm{i}\widehat{\mathrm{\Lambda }}\right)\right\}$$
(16)
The integration over the unitary group $`d[\widehat{U}]`$ in Eq. (16) can be carried out exactly using the integral first computed by Harish-Chandra  and known in the physics community as the Itzykson-Zuber integral.  As a result we obtain

$$Z_0^\alpha (\mu )=c_{\alpha \times 2M}\left(\frac{\pi }{N}\right)^{\alpha M(\alpha 2M1)}d[\widehat{\mathrm{\Lambda }}]\left(\frac{\mathrm{\Delta }(\widehat{\mathrm{\Lambda }})}{\mathrm{\Delta }(\widehat{\epsilon })}\right)\mathrm{exp}\left\{N\mathrm{Tr}\left(\mathrm{ln}\left(\mathrm{i}\widehat{\mathrm{\Lambda }}\right)+\frac{1}{2}(\widehat{\mathrm{\Lambda }}\widehat{\epsilon }+\mathrm{i}\mu 𝟙)^\mathrm{𝟚}\right)\right\}.$$
(17)
The Vandermonde determinant $`\mathrm{\Delta }(\widehat{\epsilon })=^{}(\epsilon _n\epsilon _m)`$ vanishes if the number of replicas $`\alpha `$ is greater than unity because of the presence of identical Matsubara frequencies in different replicas. Therefore the expression in Eq. (17) appears to diverge. This is not so however, since the integral over $`d[\widehat{\mathrm{\Lambda }}]`$ also vanishes when eigenvalues of the matrix $`\widehat{\epsilon }`$ become degenerate. To regularize the integral in Eq. (17) we introduce a diagonal regulator matrix $`\widehat{V}=\delta _{ij}\delta _{mn}V_i`$ with $`V_iV_j`$, and calculate the expression

$$Z_0^\alpha (\widehat{V})=c_{\alpha \times 2M}\left(\frac{\pi }{N}\right)^{\alpha M(\alpha 2M1)}d[\widehat{\mathrm{\Lambda }}]\left(\frac{\mathrm{\Delta }(\widehat{\mathrm{\Lambda }})}{\mathrm{\Delta }(\widehat{\epsilon }+\widehat{V})}\right)\mathrm{exp}\left\{N\underset{\mathrm{j},\mathrm{m}}{}\mathrm{F}(\lambda _\mathrm{m}^\mathrm{j},\mathrm{x}_{\mathrm{m},\mathrm{j}})\right\},$$
(18)
where the function $`F(\lambda _m^j,x_{m,j})`$ is defined by the equation

$`F(\lambda _m^j,x_{m,j})`$ $`=`$ $`\mathrm{ln}(i\lambda _m^j)+{\displaystyle \frac{1}{2}}(\lambda _m^jx_{m,j})^2,`$ (20)
$`x_{m,j}`$ $``$ $`\epsilon _mi\mu +V_j.`$ (21)
The ratio of the two Vandermonde determinants in Eq. (18) is rendered finite due to the presence of the regulators $`V_i`$ in Eq. (21):

$$\frac{\mathrm{\Delta }(\widehat{\mathrm{\Lambda }})}{\mathrm{\Delta }(\widehat{\epsilon }+\widehat{V})}=\underset{\{i,m\}\{j,n\}}{}\frac{\lambda _m^i\lambda _n^j}{x_{m,i}x_{n,j}}.$$
(22)
We therefore carry out the integration over $`d[\widehat{\mathrm{\Lambda }}]`$ in Eq. (18) at finite $`\widehat{V}`$ and take the limit $`\widehat{V}0`$ at the end of the calculations.

For $`N1`$ the integration over $`d[\widehat{\mathrm{\Lambda }}]`$ in the r.h.s. of Eq. (18) can be performed by the saddle point method.  To find the saddle point value of $`\lambda _m^i`$ at $`N1`$ we may ignore the ratio of the Vandermonde determinants in the pre-exponential factor in the right hand side of Eq. (18). This can be done even after we let the number of Matsubara frequencies $`2M`$ go to infinity, as we shall do eventually. We show in Appendix A that the shift of the saddle point equation which arises from the pre-exponential factor is small as $`\delta /T`$ and may be neglected. We denote the saddle point value of $`\lambda _m^i`$ by $`\overline{\lambda }_m^i`$. Minimization of the right hand side of Eq. (20) leads to the following saddle point equation

$$\overline{\lambda }_m^i\left(\overline{\lambda }_m^ix_{m,i}\right)=1.$$
(23)
This equation has two solutions for each $`\overline{\lambda }_m^i`$.

$$\overline{\lambda }_m^i=\lambda _\pm (x_{m,i})=\frac{x_{m,i}}{2}\pm \sqrt{1+\frac{x_{m,i}^2}{4}}.$$
(24)
Therefore there are $`2^{2M\alpha }`$ saddle points of the integrand in Eq. (18).

It is clear from Eq. (18) that the fluctuations of $`\lambda _m^i`$ about its saddle point value, Eq. (24) are of order $`1/\sqrt{N}`$. Let us assume for a moment that both the temperature $`T`$ and $`V_i`$ satisfy the condition $`1/\sqrt{N}T,V_iV_j1`$. Then the ratio of the Vandermonde determinants in the pre-exponential factor in Eq. (18) can be replaced by its value at the saddle point

$$\frac{\mathrm{\Delta }(\widehat{\mathrm{\Lambda }})}{\mathrm{\Delta }(\widehat{\epsilon }+\widehat{V})}|_{\mathrm{sp}}=\underset{\{i,m\}\{j,n\}}{}\frac{\overline{\lambda }(x_{m,i})\overline{\lambda }(x_{n,j})}{x_{m,i}x_{n,j}}.$$
(25)
Next, we formally extend the product in this equation to include the diagonal terms $`\{i,m\}=\{j,n\}`$ by replacing the ratio in the r.h.s. of Eq. (25) by $`d\overline{\lambda }(x_{m,i})/dx_{m,i}`$ to obtain

$$\frac{\mathrm{\Delta }(\widehat{\mathrm{\Lambda })}}{\mathrm{\Delta }(\widehat{\epsilon }+\widehat{V})}|_{\mathrm{sp}}=\mathrm{exp}\left\{\frac{1}{2}\underset{im,nj}{}\mathrm{ln}\left(\frac{\overline{\lambda }(x_{m,i})\overline{\lambda }(x_{n,j})}{x_{m,i}x_{n,j}}\right)\right\}\underset{i,m}{}\left(\frac{d\overline{\lambda }(x_{m,i})}{d\epsilon _m}\right)^{1/2},$$
(26)
where the sum over the replica and the Matsubara indices is unrestricted and includes the terms with $`\{i,m\}=\{j,n\}`$. The factor $`1/2`$ in the exponent in Eq. (26) appears due to double counting of terms in the summation.

We now substitute Eq. (26) into Eq. (18) and expand the function $`F(\lambda _m^i,\epsilon _m,V_i)`$ in Eq. (20) to second order in the fluctuations around the saddle point Eq. (24). It is easy to show that

$$\frac{^2F(\lambda _m^i,x_{m,i})}{\lambda _{m}^{i}{}_{}{}^{2}}|_{\overline{\lambda }}=\left(\frac{\overline{\lambda }(x_{m,i})}{\epsilon _m}\right)^1.$$
(27)
Therefore the Gaussian integration over the fluctuations around the saddle point, Eq. (24) results in the product $`\underset{i.m}{}[(2\pi /N)d\overline{\lambda }(x_{m,i})/d\epsilon _m]^{1/2}`$ which partially cancels the last product in Eq. (26), and we obtain

$`Z_0^\alpha (\widehat{V})`$ $`=`$ $`{\displaystyle \underset{\{\overline{\lambda }\}}{}}\mathrm{exp}\left\{N{\displaystyle \underset{m,i}{}}F(\overline{\lambda }(x_{m,i}),x_{m,i})+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j;n,m}{}}\mathrm{ln}\left({\displaystyle \frac{\overline{\lambda }(x_{m,i})\overline{\lambda }(x_{n,j})}{x_{m,i}x_{n,j}}}\right)\right\},`$ (28)
where $`\underset{\{\overline{\lambda }\}}{}`$ denotes the sum over all the saddle points given by Eq. (24).

Although, to justify the saddle point procedure and to arrive at Eq. (28) we assumed $`V_iV_j1/\sqrt{N}`$, the saddle point result (28) is exact for the unitary ensemble and holds the way to $`V_i0`$. It was shown by Zirnbauer  that this special feature of the unitary ensemble is a consequence of the Duistermaat-Heckman theorem.

Equation (28) is free of divergences in the $`\widehat{V}0`$ limit. For the replica-symmetric saddle points it is obvious because in the “dangerous” terms with $`m=n`$ the argument of the logarithm in Eq. (28) remain finite and equal to $`\overline{\lambda }(\epsilon _mi\mu )/\epsilon _m`$ in the $`V_i0`$ limit. For the saddle points with broken replica symmetry it is necessary to sum the contributions of all replica permutations of a given saddle point in order to obtain a finite result.

Next, we take the number $`2M`$ of the Matsubara frequencies to infinity. The resulting infinite sums over the Matsubara frequencies in Eq. (28) diverge and need to be regularized. The method of regularization can be inferred from the observation that the sums over the Matsubara frequencies should be regarded as traces of operators in the Matsubara space. Let us consider the one particle Matsubara Green function as an example

$$G_{x,y}(\tau _1\tau _2)=\psi _x(\tau _1)\overline{\psi }_y(\tau _2),$$
(29)
where the average is taken over the functional integral as in Eq. (2). It is important to remember that the Fermion operators $`\widehat{\psi }`$ and $`\widehat{\overline{\psi }}`$ appearing in the Hamiltonian which enter the operator expression for the partition function $`Z_0(\mu ,T)`$ are normal ordered. Therefore, in the time-discretized version of the functional integral for the partition function, Eq. (2) the Hamiltonian $`\widehat{H}_0`$ is coupled to fermionic fields at different moments of time: $`\overline{\psi }(t_{i+1})\widehat{H}_0\psi (t_i)`$. This forces us to understand the traces of the Greens function as

$$\mathrm{Tr}\widehat{\mathrm{G}}=\underset{0}{\overset{\beta }{}}d\tau \underset{\eta +0}{lim}\underset{\mathrm{x}=1}{\overset{\mathrm{N}}{}}\mathrm{G}_{\mathrm{x},\mathrm{x}}(\eta ),$$
(30)
which coincides with the total number of particles in the system. In the frequency representation the regularized trace of the Green function is given by

$$\mathrm{Tr}\widehat{\mathrm{G}}=\mathrm{T}\underset{\eta +0}{lim}\underset{\epsilon _\mathrm{n}=\mathrm{}}{\overset{+\mathrm{}}{}}\underset{\mathrm{x}=1}{\overset{\mathrm{N}}{}}\mathrm{G}_{\mathrm{x},\mathrm{x}}(\epsilon _\mathrm{n})\mathrm{e}^{\mathrm{i}\epsilon _\mathrm{n}\eta }.$$
(31)
The exponential factor $`e^{i\epsilon _n\eta }`$ in Eq. (31) is essential for the convergence of the sum over the Matsubara frequencies. Since the Hamiltonian $`\widehat{H}_0`$ in Eq. (2), over which the averaging is performed, is normal ordered the traces of all the ensuing operators should be regularized by the introduction of the factors $`e^{i\epsilon _n\eta }`$ as in Eq. (31).

Thus, Eq. (28), where the summations over Matsubara frequencies should be understood according to the rule $`_nf(n)_{n=\mathrm{}}^{\mathrm{}}f(n)\mathrm{exp}(i\eta \epsilon _n)`$, represents an exact expression for the averaged replicated partition function for $`N\mathrm{}`$.

Below we will evaluate Eq. (28) at $`T\delta `$. In this case the sum in Eq. (28) is dominated by the saddle point which corresponds to the lowest value of the exponent in Eq. (28). The contributions of all the other saddle points are exponentially small in $`T/\delta `$. The leading term in the sum in Eq. (28) corresponds to the saddle point

$$\overline{\lambda }_m^i=\lambda ^0(x_{m,i})=\frac{x_{m,i}}{2}+\sqrt{1+\frac{x_{m,i}^2}{4}}.$$
(32)
Here the function $`\lambda ^0(z)=z/2+\sqrt{1+z^2/4}`$ is understood as an analytic function of $`z`$ in the complex plane with a branch cut from $`2i`$ to $`2i`$ as in Fig. 2, in particular the square root in Eq. (32) changes sign when the real part of $`z`$, $`\mathrm{}z`$, crosses zero.

To obtain the expression for $`Z_0^\alpha (\mu ,T)`$ we take the limit $`V_i0`$ in Eq. (28) and perform the summation over the replica indices. As a result we obtain

$`Z_0^\alpha (\mu )`$ $`=`$ $`\mathrm{exp}\{\underset{\eta +0}{lim}[N\alpha {\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}F(\lambda ^0(\epsilon _mi\mu ),\epsilon _mi\mu )e^{i\epsilon _m\eta }`$ (33)
$`+`$ $`{\displaystyle \frac{\alpha ^2}{2}}{\displaystyle \underset{m,n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{i(\epsilon _m+\epsilon _n)\eta }\mathrm{ln}\left({\displaystyle \frac{\lambda ^0(\epsilon _mi\mu )\lambda ^0(\epsilon _ni\mu )}{\epsilon _m\epsilon _n}}\right)]\}.`$ (34)
The first term in the exponent in this equation is linear in the number of replicas $`\alpha `$ and corresponds to the average grand canonical potential $`\mathrm{\Omega }_0(T,\mu )`$ which is evaluated in section II A. The second term in the exponent of Eq. (34) is quadratic in $`\alpha `$ and describes the mesoscopic fluctuations of $`\mathrm{\Omega }_0(T,\mu )`$ which are studied in section II B.

### A Average grand canonical potential

In this section we evaluate the average grand canonical potential $`\mathrm{\Omega }_0(\mu )`$. It is easy to show using Eqs. (7) and (34) that it is given by

$$\mathrm{\Omega }_0(\mu )=NT\underset{\eta 0}{lim}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}F(\lambda ^0(\epsilon _mi\mu ),\epsilon _mi\mu )e^{i\epsilon _m\eta }.$$
(35)
The function $`F(\lambda ^0(\epsilon _mi\mu ),\epsilon _mi\mu )`$ is defined through Eqs. (20) and (32). It has a branch cut from $`2i+i\mu `$ to $`2i+i\mu `$ in the complex plane of the variable $`\epsilon _m`$ and is analytic everywhere else. It is convenient to rewrite this sum in terms of the contour integral

$$\mathrm{\Omega }_0(\mu )=\frac{N}{4\pi i}\underset{\eta 0}{lim}\underset{C}{}𝑑\epsilon F(\lambda ^0(\epsilon _mi\mu ),\epsilon _mi\mu )e^{i\epsilon \eta }\left(\mathrm{tan}\left(\frac{\epsilon }{2T}\right)+i\right),$$
(36)
where the contour $`C`$ is shown in Fig. 2. Taking into account the analytic properties of $`F(\lambda ^0(\epsilon _mi\mu ),\epsilon _mi\mu )`$ we can deform the integration contour to $`C^{}`$ and integrate once by parts to obtain

$$\mathrm{\Omega }_0(\mu )=\frac{NT}{2\pi i}\underset{C^{}}{}𝑑\epsilon \left[\frac{\epsilon i\mu }{2}\sqrt{1+\frac{(\epsilon i\mu )^2}{4}}\right]\mathrm{ln}\left(\frac{1+e^{i\epsilon /T}}{2}\right).$$
(37)
Note that the two terms in the square brackets coincide with the trace (in the Hilbert space) of the averaged single particle Green function. The latter is equal to the retarded Green function on the right side of the branch cut and to the advanced Green function on its left side. Therefore, upon the change of variables $`\epsilon =i(E\mu )`$ we can express the integral in Eq. (37) through the average single particle density of states defined in Eq. (4)

$$\mathrm{\Omega }_0(\mu )=T\underset{2}{\overset{2}{}}𝑑E\nu (E)\mathrm{ln}\left(\frac{1+e^{(E\mu )/T}}{2}\right).$$
(38)
This result is precisely what one expects for the average grand canonical potential of non-interacting electrons.

### B Mesoscopic fluctuations of the grand canonical potential

In this section we consider the mesoscopic fluctuations of $`\mathrm{\Omega }_0(\mu )`$ at $`T\delta `$. From Eqs. (7) and (34) we immediately find that the second cumulant $`\mathrm{\Omega }_0(\mu )`$ is determined by the second term in the exponent of Eq. (34)

$$\mathrm{\Omega }_0^2(\mu )=T^2\underset{m,n=\mathrm{}}{\overset{\mathrm{}}{}}e^{i(\epsilon _m+\epsilon _n)\eta }\mathrm{ln}\left(\frac{\lambda ^0(\epsilon _mi\mu )\lambda ^0(\epsilon _ni\mu )}{\epsilon _m\epsilon _n}\right)$$
(39)
As in the previous section it is convenient to rewrite the summation over the Matsubara frequencies in the last formula in terms of the contour integral

$$\mathrm{\Omega }_0^2(\mu )=\frac{1}{(4\pi i)^2}\underset{C_1}{}\underset{C_2}{}𝑑\epsilon _1𝑑\epsilon _2\mathrm{ln}\left(\frac{\lambda ^0(\epsilon _mi\mu )\lambda ^0(\epsilon _ni\mu )}{\epsilon _m\epsilon _n}\right)\left(\mathrm{tan}\frac{\epsilon _1}{2T}+i\right)\left(\mathrm{tan}\frac{\epsilon _2}{2T}+i\right)e^{i(\epsilon _1+\epsilon _2)\eta }$$
(40)
where the contours of integration $`C_1`$ and $`C_2`$ are shown in Fig. 3.

In Eq. (40) the function $`\lambda ^0(zi\mu )`$ was defined in Eq. (32) and is understood as a function of $`z`$ which is analytic everywhere in the complex plane except for a branch cut from $`2i+i\mu `$ to $`2i+i\mu `$ as in Fig. 3. We therefore can deform the integration contours $`C_1`$ and $`C_2`$ to $`C_1^{^{}}`$, $`C_2^{^{}}`$. After that we integrate by parts once with respect to both $`\epsilon _1`$ and $`\epsilon _2`$ and obtain

$$\mathrm{\Omega }_0^2(\mu )=\frac{T^2}{(2\pi )^2}\underset{C_1^{^{}}}{}\underset{C_2^{^{}}}{}𝑑\epsilon _1𝑑\epsilon _2G_c(i\epsilon _1+\mu ,i\epsilon _2+\mu )\mathrm{ln}\left(\frac{1+e^{i\epsilon _1/T}}{2}\right)\mathrm{ln}\left(\frac{1+e^{i\epsilon _2/T}}{2}\right),$$
(41)
where we introduced the function

$$G_c(i\epsilon _1+\mu ,i\epsilon _2+\mu )=\frac{^2}{\epsilon _1\epsilon _2}\mathrm{ln}\left(\frac{\frac{\epsilon _1i\mu }{2}+\sqrt{1+\frac{(\epsilon _1i\mu )^2}{4}}\frac{\epsilon _2i\mu }{2}\sqrt{1+\frac{(\epsilon _2i\mu )^2}{4}}}{\epsilon _1\epsilon _2}\right).$$
(42)
The function $`G_c(\epsilon _1,\epsilon _2)`$ is equal to the smooth part (i.e. not containing the oscillations on the scale of the mean level spacing $`\delta `$) of the connected two point correlator of Green functions at energies $`i\epsilon _1+\mu `$ and $`i\epsilon _2+\mu `$ which was obtained by Brezin and Zee.

Using the fact that the square root in $`\lambda _0(\epsilon _ii\mu )`$ has a different sign on the different sides of the contour $`C_i^{^{}}`$ we can express $`\mathrm{\Omega }_0(\mu )^2`$ through the integral over, say the right side of the contours $`C_1^{^{}}`$ and $`C_2^{^{}}`$ only. Introducing the variables $`E_1=i\epsilon _1+\mu `$ and $`E_2=i\epsilon _2+\mu `$ we obtain

$$\mathrm{\Omega }_0^2(\mu )=T^2\underset{2}{\overset{2}{}}\underset{2}{\overset{2}{}}𝑑E_1𝑑E_2\stackrel{~}{\rho }_c(E_1,E_2)\mathrm{ln}\left(\frac{1+e^{(E_1\mu )/T}}{2}\right)\mathrm{ln}\left(\frac{1+e^{(E_2\mu )/T}}{2}\right),$$
(43)
where $`\stackrel{~}{\rho }_c(E_1,E_2)`$ denotes the smooth part of the connected density of states correlator

$$\stackrel{~}{\rho }_c(E_1,E_2)=\frac{E_1E_24}{2\pi ^2(E_1E_2)^2\sqrt{(4E_1^2)(4E_2^2)}}.$$
(44)
This expression has a singularity when $`E_1E_20`$ which leads to the divergence of the integral in Eq. (43). If we replace $`\stackrel{~}{\rho }_c(E_1,E_2)`$ in Eq. (43) by the exact connected density correlation function $`\rho _c(E_1,E_2)`$ which does not diverge at $`E_1E_20`$ the integral in Eq. (43) becomes finite and gives the exact expression for the second cumulant of the grand canonical potential. Since we have retained only one term corresponding to the lowest energy saddle point $`\{\lambda _0\}`$ in the sum over the saddle points in Eq. (28), the expression in Eq. (43) is missing the oscillatory part of the two point level density correlator. These oscillatory terms arise from the other saddle points in Eq. (28) and cut off the divergence of the integral at the scale of the mean level spacing $`E_1E_2\delta 1/N`$. Therefore the integral in Eq. (43) should be cut off at $`E_1E_21/N`$ which leads to the variance $`\mathrm{\Omega }_0^2(\mu )N^2\mathrm{ln}N`$. This result is to be expected since the main contribution to the thermodynamic potential arises from levels deep below the Fermi surface. The fluctuation of the number of such levels is proportional to $`\mathrm{ln}N`$ and their typical energy is of order $`N`$.

## III Coulomb interaction and tunneling

In this section we consider a disordered metallic grain coupled by a tunneling contact to a clean metallic lead and capacitively coupled to a metallic gate as in Fig. 1. The Thouless energy of the grain is assumed to be greater than the charging energy, $`E_T>>E_C`$. In addition, we assume that the dot is placed in a magnetic field, such that the cooperon gap exceeds the charging energy $`DeH/(\mathrm{}c)E_C`$. It was shown by Efetov  that under these conditions the ensemble of disorder potentials in the grain can be described by the unitary random matrix ensemble. At the same time we assume that the Zeeman splitting of electron energy levels, $`\mathrm{}eH/mc`$ is smaller than the temperature $`T`$ and may be neglected. We therefore assume that each orbital level is doubly degenerate.

The single particle levels in the grain are broadened due to tunneling into the lead. We assume that the tunneling contact is broad, so that each single particle state in the dot is coupled to many lead states. In this case the probability distribution $`P(\mathrm{\Gamma }_x)`$ of level half-widths $`\mathrm{\Gamma }_x`$ is sharply peaked about the mean value $`\mathrm{\Gamma }`$. We can therefore neglect the fluctuations of $`\mathrm{\Gamma }_x`$ and consider them to be equal to the mean value $`\mathrm{\Gamma }`$ which can be expressed through the dimensionless conductance of the contact $`g_T`$ as

$$\mathrm{\Gamma }=\frac{g_T\delta }{8\pi }.$$
(45)
We assume that the temperature satisfies the conditions $`\delta TE_C`$ and focus on the Coulomb blockade oscillations in the thermodynamic quantities of the grain at $`g_T1`$ and at a finite mean level spacing $`\delta `$. The quantity of particular interest is the differential capacitance of the dot $`C(q)`$ which can be expressed through the grand canonical potential $`\mathrm{\Omega }(q)`$ of the dot at gate voltage $`q`$ as

$$C(q)=\kappa \frac{^2\mathrm{\Omega }(q)}{q^2},$$
(46)
where $`\kappa `$ is a factor depending on the geometry of the dot. In addition to the average differential capacitance and the grand canonical potential, Eq. (7), we study correlations of thermodynamic quantities at different values of the gate voltage, $`q`$ and $`q^{}`$. For this purpose we calculate the product of replicated partition functions at different values of the gate voltage $`Z^\alpha (q)Z^\alpha ^{}(q^{})`$ averaged over the ensemble of random Hamiltonians with the distribution (3). The correlators $`\mathrm{\Omega }^n(q)\mathrm{\Omega }^m(q^{})`$ can then be obtained from the formula

$$\mathrm{\Omega }^n(q)\mathrm{\Omega }^m(q^{})=(T)^{n+m}\left(\frac{}{\alpha }\right)^n\left(\frac{}{\alpha ^{}}\right)^mZ^\alpha (q)Z^\alpha ^{}(q^{})|_{\alpha ,\alpha ^{}0}.$$
(47)
### A Average thermodynamic potential

In the absence of the Coulomb interaction but at a finite tunneling rate $`\mathrm{\Gamma }`$ we can write the partition function $`Z_\mathrm{\Gamma }^0(\mu ,T)`$, where $`\mu `$ is the chemical potential of electrons in the lead, as a functional integral over the fermionic fields $`\psi `$

$$Z_\mathrm{\Gamma }^0(\mu ,T)=\underset{x=1}{\overset{N}{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}d\overline{\psi }_{x,n}^\sigma d\psi _{x,n}^\sigma \mathrm{exp}\left[\underset{n,m=\mathrm{}}{\overset{\mathrm{}}{}}\underset{x,y=1}{\overset{N}{}}\underset{\sigma }{}\overline{\psi }_{x,n}^\sigma \left(i\widehat{\epsilon }+i\mathrm{\Gamma }\mathrm{sgn}(\widehat{\epsilon })\delta _{\mathrm{nm}}(\widehat{\mathrm{H}}_{0,\mathrm{xy}}\mu \delta _{\mathrm{xy}})\right)\psi _{y,m}^\sigma \right].$$
(48)
Here the index $`\sigma `$ denotes the spin of the particle.

In the presence of the Coulomb interaction, Eq. (1), the replicated partition function $`Z^\alpha (q)`$ may be written as a functional integral over the fermionic fields similar to that in Eq. (48). However, in this case the exponent acquires a quartic term in the fermion fields $`\psi `$. A convenient way to proceed is to decouple the interaction via the Hubbard-Stratonovich transformation by introducing an auxiliary field $`V(\tau )`$. Following Finkelstein  we introduce such auxiliary fields $`V_j`$ in each replica,

$$\mathrm{exp}\left\{\underset{i=1}{\overset{\alpha }{}}\underset{0}{\overset{\beta }{}}E_C(N_j(\tau )q)^2𝑑\tau \right\}=d[V]\mathrm{exp}\left\{\underset{j=1}{\overset{\alpha }{}}\left(\underset{0}{\overset{\beta }{}}𝑑\tau \frac{V_j^2(\tau )}{4E_C}i\underset{0}{\overset{\beta }{}}V_j(\tau )(N_j(\tau )q)𝑑\tau \right)\right\},$$
(49)
where $`N_j(\tau )_{x,\sigma }\overline{\psi }_{x,j}^\sigma (\tau )\psi _{x,j}^\sigma (\tau )`$. Next we denote the static part of $`V_j(\tau )`$ by $`V_j`$ and write the integral over $`V_j`$ as a sum of integrals labeled by a set of winding numbers $`\{W\}`$

$`{\displaystyle \underset{j}{}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑V_j={\displaystyle \underset{\{W\}}{}}{\displaystyle \underset{j}{}}{\displaystyle \underset{\pi T(2W_j1)}{\overset{\pi T(2W_j+1)}{}}}𝑑V_j.`$For each set of winding numbers we express $`V_j(\tau )`$ as

$$V_j(\tau )=V_j2\pi TW_j+\dot{\varphi }_j(\tau ),$$
(50)
where the phases $`\varphi _j(\tau )`$ satisfy the periodicity condition $`\varphi _j(\beta )=\varphi _j(0)+2\pi W_j`$. Then, making the gauge transformation $`\psi (\tau )\psi (\tau )\mathrm{exp}(i\varphi (\tau ))`$, using Eqs. (49), (48) and averaging over the random matrix $`\widehat{H}_0`$ with the distribution function defined in Eq. (3) we write the ensemble averaged replicated partition function $`Z^\alpha (q)`$ as a sum over the winding numbers

$$Z^\alpha (q)=\underset{\{W\}}{}\underset{j=1}{\overset{\alpha }{}}\underset{\pi T(2W_j1)}{\overset{\pi T(2W_j+1)}{}}\frac{dV_j}{\sqrt{4\pi E_cT}}e^{\frac{V_j^2}{4TE_c}iq\frac{V_j}{T}}D\widehat{Q}d[\varphi ]e^{\frac{N}{2}\mathrm{Tr}\widehat{Q}^2+N\mathrm{Tr}\mathrm{ln}\left(i\widehat{Q}+i\widehat{J}\right)\underset{0}{\overset{\beta }{}}\frac{(\dot{\varphi }_j2\pi TW_j)^2d\tau }{4E_C}}.$$
(51)
Here the trace is taken over the spin, the replica and the Matsubara space. The $`\widehat{Q}`$-matrix was introduced to decouple the quartic interaction of fermions arising from the averaging over $`\widehat{H}_0`$, $`D\widehat{Q}`$ denotes the normalized integration measure, and we introduced the operator

$`\widehat{J}=\widehat{\epsilon }i\mu +\left[V_j2\pi TW_j+\mathrm{\Gamma }e^{i\varphi _j(\tau )}\mathrm{\Lambda }_{\tau ,\tau ^{}}e^{i\varphi _j(\tau ^{})}\right]\delta _{ij},`$where $`\widehat{\epsilon }=i_\tau `$, and $`\mathrm{\Lambda }_{\tau ,\tau ^{}}=i/\mathrm{sin}[\pi T(\tau \tau ^{})]`$.

To perform the integration over $`\widehat{Q}`$ in Eq. (51) we proceed as in section II: we make a change of variables $`\widehat{Q}\widehat{Q}^{}\widehat{J}`$ and parameterize the matrix $`\widehat{Q}^{}`$ in terms of the unitary matrix $`\widehat{U}`$ and the diagonal matrix $`\widehat{\mathrm{\Lambda }}`$ as $`\widehat{Q}^{}=\widehat{U^1}\widehat{\mathrm{\Lambda }}\widehat{U}`$. The ensuing integral over the unitary matrices $`\widehat{U}`$ is of the Itzykson-Zuber type and can be evaluated exactly. The result is expressed through the eigenvalues $`J_m^i`$ of the operator $`\widehat{J}`$. The operator $`\widehat{J}`$ is not diagonal in the Matsubara space but may still be diagonalized using a unitary transformation. The index $`m`$ can be identified with the Matsubara index only if $`\dot{\varphi }=2\pi TW`$. We then perform the integration over $`\widehat{\mathrm{\Lambda }}`$ by the saddle point method which for the unitary ensemble gives an exact result for $`N1`$. The saddle point equation for the eigenvalues $`\lambda _{m,\sigma }^i`$ is given by $`\overline{\lambda }_{m,\sigma }^i(\overline{\lambda }_{m,\sigma }^iJ_m^i)=1`$ and has two solutions $`\lambda _\pm (J_m^i)`$ defined in Eq. (24). For $`T\delta `$ the integration over the static components $`V_j`$ can be carried out in the saddle point approximation.  The saddle point equation reads

$$iNT\underset{n,\sigma }{}(\lambda _{n,\sigma }^jJ_n^j)e^{i\eta \epsilon _n}=q\frac{iV_j}{2E_C}.$$
(52)
The left hand side of Eq. (52) may be interpreted as the number of particles in the dot with the chemical potential $`\mu +i(V_j2\pi TW_j)`$. Note that in order to have a reasonably small charging energy the imaginary part of the expression in the r.h.s. of this equation should be at most of order unity. Since the prefactor in the l.h.s. $`NT=\pi T/\delta `$ is much greater than unity according to our assumptions, the saddle point must have an equal number of $`\lambda _+`$ and $`\lambda _{}`$ in order to satisfy (52). Such a saddle point with the lowest free energy is given by

$$\lambda _0(J_m^i)=\frac{J_m^i}{2}+\mathrm{sgn}\epsilon _m\sqrt{1+\frac{(J_m^i)^2}{4}}.$$
(53)
All other saddle points satisfying Eq. (52) are given permutations between $`\lambda _+`$ and $`\lambda _{}`$ in (53) such that the numbers of $`\lambda _+`$ and of $`\lambda _{}`$ are equal. These saddle points provide only small contributions to the partition function $`\mathrm{exp}(T/\delta )`$ and will be neglected.

We first solve Eq. (52) for $`V_j`$ at $`\varphi _j(\tau )=2\pi W_j\tau `$. In this case the operator $`\widehat{J}`$ is diagonal in the Matsubara basis and its eigenvalues are given by

$`J_m^j=\epsilon _mi\mu +V_j2\pi TW_j+\mathrm{\Gamma }\mathrm{sgn}\left(\epsilon _m2\pi TW_j\right).`$Substituting this expression and (53) into Eq. (52) we find in the leading order in $`\delta /E_C`$ and $`\delta /T`$

$$V_j=2\pi TW_ji\frac{\delta }{2}\left[qN_0(\mu )\right],$$
(54)
where $`N_0(\mu )`$ is the number of particles in the dot at the chemical potential $`\mu `$ given by the l.h.s. of (52) at $`V_j=W_j=0`$. To arrive at Eq. (54) we have used that $`N_0(\mu )/\mu =2/\delta `$. It is easy to convince oneself that Eq. (54) holds for an arbitrary phase configuration $`\varphi _j`$ in the topological class $`W_j`$. We therefore can write the operator $`\widehat{J}`$ at the saddle point (53) as

$$\widehat{J}=\widehat{\epsilon }i\frac{\delta }{2}\left[qN_0(\mu )\right]+\mathrm{\Gamma }e^{i\varphi _j(\tau )}\mathrm{\Lambda }_{\tau ,\tau ^{}}e^{i\varphi _j(\tau ^{})}\delta _{ij}.$$
(55)
Retaining only the leading saddle point (53) for each set of winding numbers $`\{W\}`$ and integrating over the gaussian fluctuations of $`V_j`$ about the saddle point solutions (54) and using the fact that $`\underset{0}{\overset{\beta }{}}(\dot{\varphi }_j2\pi TW_j)^2𝑑\tau =\underset{0}{\overset{\beta }{}}\dot{\varphi }_j^2𝑑\tau +4\pi ^2TW_j^2`$ we can write the replicated partition function as a sum over the winding numbers

$`Z^\alpha (q)`$ $`=`$ $`\left({\displaystyle \frac{\delta }{4E_C}}\right)^{\alpha /2}{\displaystyle \underset{\{W\}}{}}{\displaystyle d[\varphi ]Z_{\{W\}}^\alpha e^{\underset{0}{\overset{\beta }{}}\frac{\dot{\varphi }_j^2d\tau }{4E_C}}},`$ (56)
where $`Z_{\{W\}}^\alpha `$ and $`\mathrm{\Delta }_{\{W\}}[\varphi ]`$ are given by the expressions

$`Z_{\{W\}}^\alpha `$ $`=`$ $`\mathrm{\Delta }_{\{W\}}[\varphi ]e^{_j\left[2\pi iqW_j+_{\{W\}}(\varphi _j)/T\right]},`$ (57)
$`\mathrm{\Delta }_{\{W\}}[\varphi ]`$ $`=`$ $`{\displaystyle \underset{i,j;n,m}{}}\left({\displaystyle \frac{\lambda _0(J_m^i)\lambda _0(J_n^j)}{J_m^iJ_n^j}}\right)^2.`$ (58)
Here we introduced the notation $`_{\{W\}}(\varphi _j)`$ for the free energy of the dot at the saddle point (53)

$$_{\{W\}}(\varphi _j)=NT\underset{n,\sigma }{}F(\lambda _0(J_n^j),J_n^j),$$
(59)
where $`F(\lambda _0(J_n^j),J_n^j)`$ was defined in Eq. (20). The Vandermonde determinant in Eq. (58) is generally finite. This is so because the eigenvalues $`J_m^i`$ are non-degenerate due to the presence of the auxiliary fields $`\varphi _i`$.

For the winding number $`W_j`$ the minimum of the free energy $`_{\{W\}}^j(\varphi )`$ is achieved on the instanton phase configuration $`\varphi _z(\tau )`$ which can be expressed using complex variable $`u=\mathrm{exp}(i2\pi T\tau )`$ as ,

$$\mathrm{exp}\left[i\varphi _z(\tau )\right]=\underset{\rho =1}{\overset{W_j}{}}\frac{1u^1z_\rho }{1uz_\rho ^{}}.$$
(60)
Here the complex instanton parameters $`z_\rho `$ satisfy the inequality $`|z_\rho |<1`$ for $`W_j>0`$, and $`|z_\rho |>1`$ for $`W_j<0`$. The free energy $`_{\{W\}}(\varphi _z)`$ on the instanton configuration (60) is independent of $`z_\rho `$ and is given by

$$\frac{_{\{W\}}(\varphi _z)}{T}=\frac{g_T}{2}|W_j|.$$
(61)
It is clear from Eqs. (61), (57) that the largest term in the sum over $`\{W\}`$ in Eq. (56) has a monotonic dependence on the gate voltage, $`q`$ and corresponds to $`\{W\}=\{0\}`$ with all $`W_j=0`$. The leading oscillatory contributions to this sum arise from the terms with $`\{W\}=\{l\}`$, $`l=\pm 1`$, having only one non-zero $`W_j=l=\pm 1`$ and can be chosen in $`\alpha `$ ways by permutations between the replicas. Retaining only these terms we obtain

$$\frac{Z^\alpha (q)}{Z^\alpha (q)_{\{0\}}}=1+\alpha \underset{l=\pm 1}{}e^{2\pi ilq}\frac{d[\varphi ]\mathrm{\Delta }_{\{l\}}[\varphi ]\mathrm{exp}\left\{_j\left(\frac{_{\{l\}}(\varphi _j)}{T}+_0^\beta \frac{\dot{\varphi }_j^2d\tau }{4E_C}\right)\right\}}{d[\varphi ]\mathrm{\Delta }_{\{0\}}[\varphi ]\mathrm{exp}\left\{_j\left(\frac{_{\{0\}}(\varphi _j)}{T}+_0^\beta \frac{\dot{\varphi }_j^2d\tau }{4E_C}\right)\right\}},$$
(62)
In Eq. (62) we need to integrate over the fluctuations of $`\varphi `$ around the instanton configuration Eq. (60). To this end we write $`\varphi =\varphi _z+\stackrel{~}{\varphi }`$, where $`\stackrel{~}{\varphi }`$ represents the massive modes of the dissipative action in Eq. (62) The integration over the zero modes of the dissipative action (62) should be performed with the measure $`d^2z/(1|z|^2)`$ for $`|z|<1`$, and $`d^2z/[|z|^2(|z|^21)]`$ for $`|z|>1`$. All the other fluctuations $`\stackrel{~}{\varphi }`$ have a large mass of order $`g_T`$. Therefore the integration over them may be performed in the gaussian approximation. Moreover, in carrying out this integration one may replace the determinants $`\mathrm{\Delta }[\varphi ]`$ in Eq. (62) by their values $`\mathrm{\Delta }[\varphi _z]`$ on the instanton configurations (60) due to their weak dependence on $`\stackrel{~}{\varphi }`$. In this approximation the gaussian integrals over the massive modes $`\stackrel{~}{\varphi }`$ in different replicas in Eq. (62) factorize. For the replicas with vanishing winding numbers the integrals in the numerator cancel with those in the denominator. We are therefore left with the expression involving only one replica in which the instanton is present,

$$\frac{Z^\alpha (q)}{Z^\alpha (q)_{\{0\}}}=1+2\alpha \mathrm{}e^{2\pi iq}\underset{|z|<1}{}\frac{d^2zf(z)}{(1|z|^2)}\frac{\mathrm{\Delta }_{\{l\}}[\varphi _z]}{\mathrm{\Delta }_{\{0\}}[0]},$$
(63)
where $`f(z)`$ denotes the ratio of the integrals over the massive modes with and without the instanton

$$f(z)\mathrm{exp}\left\{\frac{2\pi ^2T|z|^2}{E_C(1|z|^2)}\right\}\frac{d[\stackrel{~}{\varphi }]\mathrm{exp}\left\{\frac{_{\{1\}}(\varphi _z+\stackrel{~}{\varphi })}{T}+_0^\beta \frac{(\dot{\stackrel{~}{\varphi }})^2d\tau }{4E_C}\right\}}{d[\varphi ]\mathrm{exp}\left\{\frac{_{\{0\}}(\varphi )}{T}+_0^\beta \frac{\dot{\varphi }^2d\tau }{4E_C}\right\}}.$$
(64)
To arrive at Eqs. (63), (64) we used that $`_0^\beta \frac{\dot{\varphi }_z^2d\tau }{4E_C}=\frac{2\pi ^2T|z|^2}{E_C(1|z|^2)}`$ and the fact that the terms with $`l=\pm 1`$ in Eq. (62) are complex conjugates of each other.

In order to evaluate the ratio of the Vandermonde determinants in Eq. (63) we need to find the eigenvalues of the operator (55) on the instanton configurations. This is done in Appendix B. We then show in Appendix C that in the limit $`\alpha 0`$ the ratio $`\mathrm{\Delta }_{\{1\}}[\varphi _z]/\mathrm{\Delta }_{\{0\}}[0]`$ is equal to unity.

The logarithmic divergence of the integral over $`z`$ for short instantons, $`|z|1`$, in Eq. (63) is cut off at $`1|z|^2T/E_C`$ by the first term in $`f(z)`$, Eq. (64). The ratio of the integrals over the massive modes in (64) is evaluated in Appendix D. For very long instantons, $`|z|0`$, the result is given by $`\frac{g_T^2E_C}{2\pi (\pi T+\mathrm{\Gamma })}`$, whereas for short instantons, $`|z|1`$, and $`\mathrm{\Gamma }T`$ the result is $`\frac{g_T^2E_C}{2\pi ^2T}(1\frac{\mathrm{\Gamma }}{\pi T}\left(1|z|^2\right)\mathrm{ln}(1+\mathrm{\Gamma }/\pi T))`$. To evaluate the integral (63) with logarithmic accuracy we may interpolate $`f(z)`$ for the intermediate instanton lengths between the two asymptotics as

$$f(z)=\frac{g_T^2E_C\mathrm{exp}\left\{\frac{2\pi ^2T|z|^2}{E_C(1|z|^2)}\right\}}{2\pi ^2T(1+\frac{\mathrm{\Gamma }}{\pi T}\left[1|z|^2\right])}.$$
(65)
Using this and integrating over $`z`$ in Eq. (63) we obtain for the oscillatory part of the ensemble averaged thermodynamic potential

$$\mathrm{\Omega }(q)_{\mathrm{osc}}=\frac{g_T^2E_C}{\pi }\mathrm{ln}\left[\frac{E_C}{T+\mathrm{\Gamma }}\right]\mathrm{exp}\left(\frac{g_T}{2}\right)\mathrm{cos}(2\pi q).$$
(66)
Eq. (66) is the main result of this section. It is applicable in the temperature range $`\delta TE_C`$. For $`\mathrm{\Gamma }T`$ it coincides with the one instanton approximation of Ref. . The result of Eq. (66) was obtained in the single instanton approximation which holds for $`Tg_T^2E_C\mathrm{exp}(\frac{g_T}{2})`$.

### B Correlation function $`\mathrm{\Omega }(q)\mathrm{\Omega }(q^{})`$

In this section we consider the correlation function $`\mathrm{\Omega }(q)\mathrm{\Omega }(q^{})`$ of the thermodynamic potentials at different values $`q`$ and $`q^{}`$ of the gate voltage. For this purpose we calculate the average product of the replicated partition functions at two values of the gate voltage, $`Z^\alpha (q)Z^\alpha ^{}(q^{})`$. Repeating the steps of section III A we can write $`Z^\alpha (q)Z^\alpha ^{}(q^{})`$ as a sum over the sets of winding numbers $`\{W\}`$ and $`\{W^{}\}`$

$$Z^\alpha (q)Z^\alpha ^{}(q^{})=\underset{\{W\}}{}\underset{\{W^{}\}}{}d[\varphi ]e^{\underset{j=1}{\overset{\alpha +\alpha ^{}}{}}\underset{0}{\overset{\beta }{}}\frac{\dot{\varphi }_j^2d\tau }{4E_C}}Z_{\{W\}}^\alpha Z_{\{W^{}\}}^\alpha ^{}\underset{n,m}{}\underset{i=1}{\overset{\alpha }{}}\underset{j=1}{\overset{\alpha ^{}}{}}\left(\frac{\lambda _0(J_m^i)\lambda _0(J_n^j)}{J_m^iJ_n^j}\right)^4,$$
(67)
where $`Z_{\{W\}}^\alpha `$ and $`Z_{\{W\}}^\alpha ^{}`$ are given by Eq. (57). The primed quantities in Eq. (67) refer to the set of replicas pertaining to the partition function at the gate voltage $`q^{}`$. If one replaces the product (Vandermonde determinant) in Eq. (67) by unity this equation reproduces the product of the averaged replicated partition functions. The deviations of this product from unity describe correlations of the partition functions in the different sets of replicas.

For a given set of winding numbers $`\{W\}`$ the functional integral over the phases is dominated by configurations of $`\varphi `$ close to the instanton trajectories, Eq. (60). As in the previous section, when integrating over the massive fluctuations $`\stackrel{~}{\varphi }`$ about the instanton configurations we may neglect the dependence of the Vandermonde determinant in Eq. (67) on $`\stackrel{~}{\varphi }`$ and evaluate it on the instanton configuration. The dominant contribution to the sum over $`\{W\}`$ and $`\{W^{}\}`$ in Eq. (67) has a monotonic dependence on the gate voltages $`q`$ and $`q^{}`$ and arises from the term with all $`W_j=0`$. The largest oscillatory contributions to this sum arise from the terms with $`\{W\}=\{l\}`$ and $`\{W^{}\}=\{l^{}\}`$, $`l,l^{}=\pm 1`$, having only one non-zero $`W_j=\pm 1`$ in each set of replicas. Such terms can be chosen in $`\alpha \alpha ^{}`$ ways. It can be shown, see Appendix C, that the Vandermonde determinant in Eq. (67) on the instanton configurations with $`l=l^{}=\pm 1`$ goes to unity in the limit $`\alpha ,\alpha ^{}0`$. Such terms therefore do not contribute to the irreducible correlator of the thermodynamic potentials. The contributions of the terms with $`l=l^{}=\pm 1`$ are complex conjugates of each other. Using this property, denoting the Vandermonde determinant in the $`\alpha ,\alpha ^{}0`$ limit by $`\mathrm{\Delta }_{l,l^{}}`$, and using Eq. (47) we can write the leading oscillatory term in the irreducible correlator of thermodynamic potentials as

$$\mathrm{\Omega }(q)\mathrm{\Omega }(q^{})_{\mathrm{osc}}=2\mathrm{}e^{2\pi i(qq^{})}\underset{|z^{}|,|z|<1}{}\frac{d^2zd^2z^{}f(z)f(z^{})}{(1|z|^2)(1|z^{}|^2)}\left(\frac{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}{\mathrm{\Delta }_{\{0\}}[0]}1\right).$$
(68)
Here $`f(z)`$ is given by Eq. (65) and we have made an inversion transformation for the anti-instanton variable $`z^{}1/z^{}`$. Therefore the anti-instanton coordinate also obeys $`|z|<1`$.

The ratio of the Vandermonde determinants $`\frac{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}{\mathrm{\Delta }_{\{0\}}[0]}`$ is determined by the product over the pairs $`n,m`$ such that $`\epsilon _n\times \epsilon _n<0`$, see Appendix C. We are unable to evaluate the correlator (68) for arbitrary values of the instanton parameters and the gate voltage difference $`qq^{}`$ However, for large differences of the gate voltage, $`qq^{}\mathrm{\Gamma }/\delta `$, the correlator (68) may be evaluated with logarithmic accuracy. Indeed, at such gate voltage differences even the maximal possible instanton correction to the eigenvalue difference in the denominator in Eq. (67) is relatively small, $`\delta (J_m^iJ_n^j)\mathrm{\Gamma }\delta (qq^{})`$ (see Eq. (B8)), and may be evaluated perturbatively. In this case $`\frac{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}{\mathrm{\Delta }_{\{0\}}[0]}`$ only weakly depends on the instanton parameters $`z`$ and $`z^{}`$. Therefore the integral is again dominated by short instantons $`1|z|^2T/\mathrm{\Gamma }`$ for which the eigenvalues $`J_n^j`$ are accurately described by Eq. (B8). We show in appendix C that for $`\delta (qq^{})T,\mathrm{\Gamma }`$

$$1\frac{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}{\mathrm{\Delta }_{\{0\}}[0]}=\frac{\left(\frac{4\mathrm{\Gamma }}{\pi T}\right)^2}{\left(b+\frac{1}{1|z|^2}\right)\left(b+\frac{1}{1|z^{}|^2}\right)},$$
(69)
where $`b=i\delta (qq^{})/4\pi T`$.

Substituting Eqs. (69) and (65) into (67) and integrating over $`z`$ and $`z^{}`$ we obtain with logarithmic accuracy the long range asymptotics of the oscillatory part of the irreducible correlator of the thermodynamic potentials

$$\mathrm{\Omega }(q)\mathrm{\Omega }(q^{})_{\mathrm{osc}}=\frac{g_T^4E_C^2}{2\pi ^2}\mathrm{exp}(g_T)\left(\frac{16\mathrm{\Gamma }}{\delta (qq^{})}\right)^2\mathrm{ln}^2\left(\frac{\delta (qq^{})}{\mathrm{\Gamma }+\pi T}\right)\mathrm{cos}[2\pi (qq^{})].$$
(70)
This expression represents the main result of this section.

## IV Summary and discussion

We studied the effects of Coulomb interaction on the statistics of the thermodynamic quantities in an ensemble of weakly disordered metallic grains with broken time reversal symmetry. We assumed that the Thouless energy, $`E_T`$, in the grain exceeds the charging energy, $`E_C`$, and considered the case when the grain is connected to a metallic lead by a tunneling contact with a large conductance $`g_T1`$. We found expressions for the oscillatory parts of the average thermodynamic potential $`\mathrm{\Omega }(q)_{\mathrm{osc}}`$, Eq. (66) and of the correlation function $`\mathrm{\Omega }(q)\mathrm{\Omega }(q^{})_{\mathrm{osc}}`$ at $`\delta (qq^{})T,\mathrm{\Gamma }`$, Eq. (70). The results (66) and (70), apply in the temperature range $`\delta \mathrm{\Gamma }T`$. For $`\mathrm{\Gamma }T`$ equation (66) coincides with the single-instanton contribution of Ref. . In deriving these results we have neglected the multi-instanton contributions. This approximation is valid for $`Tg_T^2E_C\mathrm{exp}(g_T/2)`$. The correlator of mesoscopic fluctuations, Eq. (70), decays as $`(qq^{})^2\mathrm{ln}^2(qq)`$ at $`\delta |qq^{}|T,\mathrm{\Gamma }`$. Eqs. (66) and (70) represent the main results of the paper.

The expressions for the oscillatory parts of the average differential capacitance $`\delta C(q)`$ and of its irreducible correlator $`\delta C(q)\delta C(q^{})_{\mathrm{osc}}`$ may be easily obtained from Eqs. (66) and (70) with the aid of Eq. (46).

To derive our main results we employed a method based on the extension of the Itzykson-Zuber integral  to the infinite dimensional unitary group. In this method the $`Q`$-matrix degrees of freedom (which describe the mesoscopic fluctuations) are integrated out exactly. As in the case $`\mathrm{\Gamma }0`$, at $`T\delta `$ the replicated partition function and its correlator can be written over the topological classes of the phase $`\varphi `$ which decouples the Coulomb interaction term, see Eqs. (56) and (67). Each topological class is characterized by the winding numbers in different replicas . At large conductance, $`g1`$ the functional integral over the phase $`\varphi `$ within each topological class can be taken by the saddle point method. The saddle points configurations of the phase are characterized by the complex instanton parameters $`z`$. At finite escape rate from the dot $`\mathrm{\Gamma }`$ the functional integral over the massive fluctuations $`\stackrel{~}{\varphi }`$ around the instanton configurations depends on the instanton length $`(1|z|^2)/T`$, see Eq. (65). This is in contrast to the $`\mathrm{\Gamma }=0`$ case, where it is independent of the instanton length. This leads to the fact that the logarithmic divergence of the integral over the zero modes of the instantons at low temperatures is cut off at the escape rate $`\mathrm{\Gamma }`$, see Eq. (66). The study of mesoscopic fluctuations of Coulomb blockade oscillations in this approach amounts to evaluating the Vandermonde determinants (58) on the instanton configurations of the phase (60). In the case of the average partition function they are equal to unity when the number of replicas is taken to zero. In the case of the correlator of the thermodynamic potentials the Vandermonde determinant differs from unity and depends on the length of the instantons. The dependence of the Vandermonde determinant on the instanton parameters $`z`$ and $`z^{}`$ results in the presence of the logarithmic factor $`\mathrm{ln}^2\frac{\delta (qq^{})}{\mathrm{\Gamma }+\pi T}`$ in Eq. (70). In the $`\sigma `$-model formulation the mesoscopic fluctuations are described by diffusons. The dependence of the diffusons on the instanton parameters is determined by the $`z`$-dependent eigenvalues of the operator $`\widehat{J}`$ given by the solutions of Eq. (B7) and is easily generalizable to time-reversal invariant non-zero-dimensional systems. The method used here may have useful applications to the $`\sigma `$-model treatment of interaction effects in granulated disordered metals.

###### Acknowledgements.

We are thankful to K. Efetov for stimulating discussions at the early stage of this work, and to I. Aleiner, A. Altland, L. Glazman, A. Kamenev, A. Larkin, K. Matveev, E. Mishchenko and Yu. Nazarov for valuable discussions. We gratefully acknowledge the support of the Graduiertenkolleg 384 and the Sonderforschungsbereich 237 and the warm hospitality of the Max Planck Institute in Dresden where part of this work was performed. This research was sponsored by the Grants DMR-9984002 and BSF-9800338, and by the A.P. Sloan and the Packard Foundations.

## A Saddle point equation

In this appendix we demonstrate that the pre-exponential factor in Eq. (18) gives rise only to a small correction to the saddle point equation, Eq. (23) and may be neglected. To show this it is convenient to rewrite the integrand in Eq. (18) in the following form

$$\mathrm{exp}\left\{N\underset{m,i}{}F(\lambda (x_{m,i}),x_{m,i})+\frac{1}{2}\underset{i,j;n,m}{}\mathrm{ln}\left(\frac{\lambda (x_{m,i})\lambda (x_{n,j})}{x_{m,i}x_{n,j}}\right)\right\}.$$
(A1)
Taking a derivative of the function in the exponent in Eq. (A1) with respect to $`\lambda (x_{m,i})`$ we obtain the saddle point equation for $`\lambda `$

$`N\left({\displaystyle \frac{1}{\lambda (x_{n,i})}}+\lambda (x_{n,i})x_{n,i}\right)`$ (A2)
$`+{\displaystyle \frac{\alpha }{2}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\lambda (x_{n,i})\lambda (x_{m,i})}}=0.`$ (A3)
The saddle point value for $`\lambda `$, Eq. (32) has been obtained under the assumption that the second term in Eq. (A3) can be neglected. To check the self-consistency of this assumption we now evaluate the second term in Eq. (A3) at the saddle point $`\lambda =\lambda ^0(x_{m,i})`$, where $`\lambda ^0`$ is defined by Eq. (32). To this end we rewrite the sum in Eq. (A3) in terms of contour integral

$`{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\lambda ^0(x_{n,i})\lambda ^0(x_{m,j})}}`$ (A4)
$`={\displaystyle \frac{1}{4\pi iT}}\underset{\eta 0}{lim}{\displaystyle \underset{C}{}}𝑑\epsilon e^{i\epsilon \eta }\left(\mathrm{tan}\left({\displaystyle \frac{\epsilon }{2T}}\right)+i\right){\displaystyle \frac{1}{\lambda ^0(x_{n,i})\lambda ^0(\epsilon )}},`$ (A5)
where the contour $`C`$ is shown in Fig. 2, and $`\lambda ^0(\epsilon )=\epsilon /2+\sqrt{1+\epsilon ^2/4}`$. Taking into account the analytic properties of $`\lambda (\epsilon )`$ we can deform the integration contour to $`C^{}`$. After a straightforward calculation we obtain

$$\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{\lambda ^0(x_{n,i})\lambda ^0(x_{m,j})}=\frac{1}{2\pi iT\lambda ^0(x_{n,i})}\underset{2}{\overset{2}{}}𝑑x\left(\mathrm{tanh}\frac{x}{2T}+1\right)\frac{\sqrt{1+x^2/4}}{x+i\frac{\lambda ^0(x_{n,i})^21}{\lambda ^0(x_{n,i})}}\frac{1}{T}\mathrm{ln}T.$$
(A6)
To obtain the last expression in Eq. (A6) we used the fact that $`\lambda ^0(x_{n,i})11`$.

A direct comparison shows that the second term in Eq. (A3), given by Eq. (A6), is smaller than the first one in the ratio $`\frac{\alpha \delta }{T}\mathrm{ln}T`$ (to obtain this estimate we used the fact that $`N=\pi /\delta `$). Therefore we may neglect the second term in Eq. (A3) and obtain the saddle point value for $`\lambda `$ given by Eq. (32).

## B Eigenvalues of $`\widehat{J}`$ on the instanton configuration

To find the eigenvalues of the operator $`\widehat{J}`$ it is convenient to express the operator $`\widehat{J}`$, Eq. (55), in terms of the complex variable $`u=\mathrm{exp}(i2\pi T\tau )`$. Introducing the notation $`\mu _q=\frac{\delta }{2}[qN_0(\mu )]`$ we write

$$\frac{\widehat{J}+i\mu _q}{2\pi T}=u_u\frac{i\mathrm{\Gamma }}{2\pi ^2T}\frac{du^{}}{u^{}}\frac{1uz^{}}{uz}\frac{\sqrt{uu^{}}}{uu^{}}\frac{u^{}z}{1u^{}z^{}}.$$
(B1)
In Eq. (B1) we used that the matrices $`\mathrm{\Lambda }_{\tau ,\tau ^{}}`$ may be written using the complex variable $`u`$ as follows: $`\mathrm{\Lambda }=\frac{2\sqrt{uu^{}}}{uu^{}}`$. The integral in Eq. (B1) is taken over the unit circle, $`|u^{}|=1`$, and the operator should be understood as acting in the space of functions spanned by the fermionic Matsubara frequencies $`u^{n1/2}`$.

For $`|z|<1`$ the matrix elements of the $`\widehat{J}`$ in the Matsubara basis are given by

$`J_{nm}`$ $`=`$ $`\left(\epsilon _ni\mu _q+\mathrm{\Gamma }\mathrm{sgn}\epsilon _n\right)\delta _{nm}`$ (B3)
$`2\mathrm{\Gamma }(1|z|^2)z^nz^m\theta (\epsilon _n)\theta (\epsilon _m),`$
with $`\theta (x)`$ being the step function. It is clear that the functions $`u^{n1/2}`$ with negative Matsubara frequencies, $`n<0`$, are eigenfunctions of the operator (B1) with eigenvalues $`J_n=\epsilon _n\mathrm{\Gamma }i\mu _q`$.

For $`n0`$ the eigenfunctions of $`\widehat{J}`$ may be written as

$$\frac{g^{(n)}(u)}{\sqrt{u}}=\underset{m0}{}g_m^{(n)}\frac{u^m}{\sqrt{u}},$$
(B4)
where $`g^{(n)}(u)`$ is a non-singular function of $`u`$ outside the unit circle, $`|u|>1`$. Using this property and acting with the operator $`\widehat{J}`$, Eq. (B1) on (B4) we can write the eigenvalue equation as

$$\left[\frac{1}{2}u_u+\frac{\mathrm{\Gamma }i\mu _qJ_n(z)}{2\pi T}\right]g^{(n)}(u)=\frac{\mathrm{\Gamma }(1|z|^2)}{\pi T(1u^1z)}g^{(n)}\left(\frac{1}{z^{}}\right).$$
(B5)
Next we introduce the Green’s function of the operator on the left hand side of this equation,

$$G^{(n)}(u,u^{})=\underset{m0}{}\frac{u^mu^m}{J_n(z)+i\mu _q\epsilon _m\mathrm{\Gamma }},$$
(B6)
and rewrite Eq. (B5) as $`g(u)=\frac{1}{2\pi i}\frac{du^{}}{u^{}}G(u,u^{})f(u^{})`$, where $`f(u^{})`$ is the right hand side of Eq. (B5). Then setting $`u=1/z^{}`$ we obtain the eigenvalue equation in the form

$$\frac{1}{1|z|^2}=\underset{k=0}{\overset{\mathrm{}}{}}\frac{2\mathrm{\Gamma }|z|^{2k}}{J_n(z)+i\mu _q\epsilon _k\mathrm{\Gamma }}.$$
(B7)
In the $`z0`$ limit we obtain that the $`n=0`$ eigenvalue is given by $`J_0(z0)=\epsilon _0i\mu _q\mathrm{\Gamma }`$ whereas all the other eigenvalues are unchanged: $`J_n(z0)=\epsilon _ni\mu _q+\mathrm{\Gamma }`$. In the opposite limit of very short instantons, $`1|z|^2\mathrm{\Gamma }/T`$ (as well as for arbitrary $`z`$ at $`\mathrm{\Gamma }T`$ ) the sum in Eq. (B7) is dominated by a single term in which the denominator is small. We then obtain for the eigenvalues with positive $`\epsilon _n`$,

$$J_n(z)=\epsilon _n+\mathrm{\Gamma }i\mu _q2\mathrm{\Gamma }|z|^{2n}(1|z|^2).$$
(B8)
## C Evaluation of the Vandermonde determinant on the instanton configurations

In this appendix we calculate the Vandermonde determinant $`\mathrm{\Delta }_{\{W\}}[\varphi ]`$, Eq. (58) on the instanton configuration of the fields $`\varphi `$, Eq. (60).

We start with the case of the average thermodynamic potential, when the instanton with a unit winding number is present only in one replica, $`j=1`$. The correction $`\delta J_{n,1}`$ to the eigenvalues of $`\widehat{J}`$ in the $`j=1`$ replica due to the presence of the instanton, Eq. (B8) is small in comparison with the band width (which is of order unity in our notations). Therefore we may write the correction to the Vandermonde determinant Eq. (58), due to the presence of the instanton as

$$\delta \mathrm{ln}\mathrm{\Delta }=4\underset{j,n,m}{}\left(\frac{\lambda _0(J_{n,1}^0)/J_{n,1}^0}{\lambda _0(J_{n,1}^0)\lambda _0(J_{m,j}^0)}\frac{1}{J_{n,1}^0J_{m,j}^0}\right)\delta J_{n,1},$$
(C1)
where $`\lambda _0(J_{m,j}^0)`$ is given by Eq. (53). The terms with $`\epsilon _m,\epsilon _n>0`$ and with $`\epsilon _m,\epsilon _n<0`$ in the right hand side of this equation vanish. Therefore, to find the dependence of the Vandermonde determinant on the instanton parameters $`z`$ we may restrict the product over $`n`$ and $`m`$ in Eq. (58) only to the terms with $`\epsilon _m\epsilon _n<0`$.

To evaluate the logarithm of the Vandermonde determinant in Eq. (58) in the $`\alpha 0`$ limit we can use the following considerations. For a given pair of Matsubara frequencies $`n`$ and $`m`$ in the sum over the replica indices $`i`$ and $`j`$ there are: i) one term where both eigenvalues $`J_n^i`$ and $`J_m^j`$ of the operator $`\widehat{J}`$ depend on the instanton parameter $`z`$; ii) $`2(\alpha 1)`$ terms where one eigenvalue depends on $`z`$ and the other one does not; iii) $`(\alpha 1)^2`$ terms where both eigenvalues are independent of $`z`$ and equal to $`J_n^0=\epsilon _n+\mathrm{\Gamma }\mathrm{sgn}\epsilon _n`$ and $`J_m^0`$. Thus in the $`\alpha 0`$ limit we obtain

$$\mathrm{ln}\mathrm{\Delta }=4\underset{n,m}{\overset{}{}}\mathrm{ln}\left\{\frac{[J_n(z)J_m(z)](J_n^0J_m^0)}{[J_n^0J_m(z)][J_n(z)J_m^0]}\right\},$$
(C2)
where the summation goes over $`n0`$, $`m<0`$, $`J_n^0=\epsilon _n+\mathrm{\Gamma }\mathrm{sgn}\epsilon _n`$ are the eigenvalues of $`\widehat{J}`$ in the absence of instantons, and $`J_n(z)`$ are given by Eq. (B8). By observing that for negative Matsubara frequencies $`J_m(z)=J_m^0`$ it is easy to see that the right hand side of this equation vanishes and $`\mathrm{\Delta }_1[\varphi _z]=1`$.

Next we turn to evaluating the Vandermonde determinant for the correlator of thermodynamic potentials, i.e. the last product in Eq. (67). In the limit $`\alpha ,\alpha ^{}0`$ it can be done using the considerations employed above. In the case when all the winding numbers $`W_j`$ have the same sign it is not difficult to see that the Vandermonde determinant is equal to unity. In the case when the instantons with both positive and negative winding numbers are present this is no longer the case. For the case when one instanton with a positive winding number, $`|z|<1`$ and one instanton with a negative winding number, $`|\zeta |>1`$ are present in the limit $`\alpha ,\alpha ^{}0`$ we obtain

$$\frac{\mathrm{\Delta }_{\{0\}}[0]}{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}=\underset{n,m}{}\left\{\frac{[J_n(z)J_m(\zeta )](J_n^0J_m^0)}{[J_n^0J_m(\zeta )][J_n(z)J_m^0]}\right\}^4.$$
(C3)
For $`\mathrm{\Gamma }T`$ the eigenvalues $`J_n(z)`$ for $`|z|<1`$ and $`n0`$ are given by Eq. (B8). The expression for $`J_m(\zeta )`$ for $`|\zeta |>1`$ and $`m<0`$ is given by $`J_m(\zeta )=\epsilon _m\mathrm{\Gamma }i\mu _q^{}+2\mathrm{\Gamma }|\zeta |^{2m}(|\zeta |^21)`$. Substituting these eigenvalues into Eq. (C3) and introducing the new variable $`z^{}=1/\zeta `$ we obtain the following expression for the exponential factor in Eq. (67)

$$\frac{\mathrm{\Delta }_{\{0\}}[0]}{\mathrm{\Delta }_{\{1,1\}}[\varphi _z]}=e^{\underset{n,m0}{}\frac{4|z|^{2n}|z^{}|^{2m}(1|z|^2)(1|z^{}|^2)}{\left[\frac{\pi T}{\mathrm{\Gamma }}(n+m+1)+i(qq^{})\frac{\delta }{4\mathrm{\Gamma }}+1\right]^2}}.$$
(C4)
The sum in the exponent can be evaluated in the limits when the lengths of the instantons $`(1|z|^2)/T`$ and $`(1|z^{}|^2)/T`$ are either shorter or longer than the inverse chemical potential difference $`1/\delta (qq^{})`$. To evaluate the correlator (70) with logarithmic accuracy we can interpolate between the two limits and obtain Eq. (69).

## D Integration over the massive modes

In this appendix we evaluate the ratio of the integrals over the massive fluctuations of $`\varphi `$ around the instanton configuration in Eq. (62). The massive modes may be integrated out in the Gaussian approximation.

We start with the limit of very long instantons, $`|z|0`$, or $`|z|\mathrm{}`$. In this case it is convenient to work in the Fourier basis. We introduce the Fourier components $`\varphi _k`$ of $`\varphi (\tau )`$ as

$$\dot{\varphi }(\tau )=i\underset{k0}{}\omega _ke^{i\omega _k\tau }\varphi _k,\omega _k=2\pi Tk,$$
(D1)
and expand the eigenvalues $`J_m^j`$ of the operator $`\widehat{J}`$, Eq. (55), to second order of perturbation theory in $`\varphi _k`$,

$$\delta J_{n,j}^{(2)}=\underset{k0}{}\frac{|\dot{\varphi }_k^j|^2}{\epsilon _n\epsilon _{nk}+\mathrm{\Gamma }(\mathrm{sgn}\epsilon _\mathrm{n}\mathrm{sgn}\epsilon _{\mathrm{n}\mathrm{k}})}.$$
(D2)
Expanding the free energy $`_{\{W\}}^j(q,\varphi )`$ given by Eqs. (59) and (18) in $`\delta J_{n,j}^{(2)}`$,

$$\delta _{\{W\}}^j=2NT\underset{n}{}\left(\frac{F(\lambda _0(J_{n,j}),J_{n,j})}{\lambda _0(J_{n,j})}\frac{\lambda _0(J_{n,j})}{J_{n,j}}+\frac{F(\lambda _0(J_{n,j}),J_{n,j})}{J_{n,j}}\right)\delta J_{n,j}^{(2)}=2TN\underset{n}{}\left(\lambda _0(J_{n,j}^0)J_{n,j}^0\right)\delta J_{n,j}^{(2)}.$$
(D3)
To arrive at Eq. (D3) we used the fact that $`F(\lambda ,J_{n,j})/\lambda =0`$ at the saddle point. Substituting Eq. (D2) into Eq. (D3) we obtain for the second order correction to the free energy

$$\delta _{\{W\}}^j=2NT\underset{n,m}{}\frac{|\dot{\varphi }_{nm}^j|^2\left(\lambda _0(J_{n,j}^0)J_{n,j}^0\right)}{\omega _{nm}+\mathrm{\Gamma }(\mathrm{sgn}\epsilon _\mathrm{n}\mathrm{sgn}\epsilon _\mathrm{m})},$$
(D4)
In Eq. (D4) the summation over $`n`$ and $`m`$ goes over $`nm`$. Next we collect the terms with $`|\dot{\varphi }_k^j|^2`$, $`k0`$. To this end we retain only the terms with $`m=n\pm k`$ in the sum over $`m`$ in the right hand side of Eq. (D4). As a result we obtain

$$\delta _{\{W\}}^j=2NT\underset{n}{}|\dot{\varphi }_k^j|^2\left(\lambda _0(J_{n,j}^0)J_{n,j}^0\right)X_{n,k},$$
(D5)
where we introduced the notation

$$X_{n,k}=\frac{1}{2\pi Tk+\mathrm{\Gamma }(\mathrm{sgn}\epsilon _\mathrm{n}\mathrm{sgn}\epsilon _{\mathrm{n}\mathrm{k}})}+\frac{1}{2\pi Tk+\mathrm{\Gamma }(\mathrm{sgn}\epsilon _\mathrm{n}\mathrm{sgn}\epsilon _{\mathrm{n}+\mathrm{k}})}.$$
(D6)
The function $`X_{n,k}`$ vanishes for $`|\epsilon _n|>2\pi Tk`$. Therefore the sum over $`n`$ in Eq. (D5) is restricted to $`|\epsilon _n|<2\pi Tk`$ and can be straightforwardly evaluated leading to the following second order correction to the free energy

$$\delta _{\{W\}}^j=2\underset{k1}{}\frac{(2\pi T)^2\mathrm{\Gamma }|k|\left(|k||W_j|\right)}{\delta (\pi T|k|+\mathrm{\Gamma })}|\varphi _k^j|^2\theta (|k||W_j|).$$
(D7)
Using Eq. (D7) we obtain for the ratio of the integrals over the Fourier components $`\varphi _k`$ with $`k>1`$ in Eq. (62),

$$\left(\frac{2\pi ^2T}{E_C}+\frac{\pi Tg_T}{\pi T+\mathrm{\Gamma }}\right)\underset{k=2}{\overset{\mathrm{}}{}}\left(\frac{\frac{2\pi ^2Tk^2}{E_C}+\frac{\pi Tg_Tk^2}{\pi T|k|+\mathrm{\Gamma }}}{\frac{2\pi ^2Tk^2}{E_C}+\frac{\pi Tg_T|k|\left(|k|1\right)}{\pi T|k|+\mathrm{\Gamma }}}\right)=\frac{g_T^2E_C}{2\pi (\pi T+\mathrm{\Gamma })}.$$
(D8)
In the limit $`\mathrm{\Gamma }T`$, Eq. (D8) reproduces the result of Ref. .

Next we evaluate the ratio of the integrals over the massive fluctuations of $`\varphi `$ in Eq. (62) in the limit of short instantons, $`1|z|^2T/\mathrm{\Gamma }`$. To this end we use Eq. (D3) to rewrite the correction to the free energy in the form

$$\delta F=NT\underset{n,\sigma }{}\mathrm{sgn}\epsilon _\mathrm{n}\delta \mathrm{J}_\mathrm{n},$$
(D9)
where $`\delta J_n`$ up to second order in $`\varphi `$ is given by

$$\delta J_n=\underset{mn}{}\frac{||\phi _n|\dot{\varphi }|\phi _m|^2}{J_n(z)J_m(z)},$$
(D10)
here $`|\phi _n`$ are the eigenfunctions of the operator $`\widehat{J}`$ and $`J_n(z)`$ are its eigenvalues. For the positive frequencies $`\epsilon _n`$ the latter are given by Eq. (B8).

For a finite length of the instanton the operator $`\widehat{J}`$ is no longer diagonal in the Fourier basis. However, in the limit of short instantons, $`1|z|^2T/\mathrm{\Gamma }`$, the non-diagonal matrix elements of $`\widehat{J}`$, Eq. (B3), are small and may be treated perturbatively. Treating the instanton correction to $`\widehat{J}`$ given by the last term in Eq. (B3), $`\gamma _{m,n}=2\mathrm{\Gamma }(z^{})^nz^m(1|z|^2)\theta (\epsilon _n)\theta (\epsilon _m)`$ as a perturbation we can write the eigenfunctions in Eq. (D10) as

$$|\phi _n=|n+\underset{mn}{}\frac{|m\gamma _{mn}}{J_n^0J_m^0},$$
(D11)
here $`J_n^0=\epsilon _n+\mathrm{\Gamma }\mathrm{sgn}\epsilon _\mathrm{n}`$ and $`|n`$ are the eigenvalues and the eigenfunctions of the operator $`\widehat{J}`$ in the absence of the instanton. Substituting Eq. (D11) into Eq. (D10) and expanding the denominator Eq. (D10) to linear order in $`\gamma _{n,n}`$ after a lengthy but straightforward calculation we obtain the correction to the free energy in Eq. (D9) in the form $`\delta F=\frac{Tg_T}{2}_{l,l^{}}\varphi _l^{}\varphi _l^{}\delta F_{l,l^{}}`$, where

$$\delta F_{l,l^{}}=\frac{\delta _{l,l^{}}l^2}{|l|+\mathrm{\Gamma }/\pi T}+\frac{|l||l^{}|\left(z^{|l|+|l^{}|}z^{\left||l||l^{}|\right|}\right)}{\left(|l|+\mathrm{\Gamma }/\pi T\right)\left(|l^{}|+\mathrm{\Gamma }/\pi T\right)}\left[\theta (l)\theta (l^{})+\theta (l)\theta (l^{})\right].$$
(D12)
Since the functional integral over the fluctuations depends only on the absolute value of $`z`$ we assumed $`z`$ to be real without loss of generality. Equation (D12) is written in the Matsubara basis. It is more convenient to use the basis of Ref.  which diagonalizes the matrix $`\delta \widehat{F}`$, Eq. (D12), at $`\mathrm{\Gamma }=0`$. Denoting the functions in this basis by $`|\stackrel{~}{\phi }_m`$, where

$$|\stackrel{~}{\phi }_m=\{\begin{array}{cc}\frac{\sqrt{1|z|^2}}{uz},\hfill & \hfill m=1\\ u^m\frac{1uz^{}}{1u^1z},\hfill & \hfill m2\\ u^m\frac{1u^1z}{1uz^{}},\hfill & \hfill m2\\ \frac{\sqrt{1|z|^2}}{u^1z^{}},\hfill & \hfill m=1,\end{array}$$
(D13)
we can express the matrix elements of $`\delta \widehat{F}`$ in this basis as

$$\delta F_{nm}=\underset{l,l^{}}{}\stackrel{~}{\phi _n}|l\delta F_{l,l^{}}l^{}|\stackrel{~}{\phi }_m,$$
(D14)
where $`l|\stackrel{~}{\phi }_m`$ are the elements of the unitary transformation matrix between the two bases,

$$l|\stackrel{~}{\phi }_m=\{\begin{array}{cc}z^{lm}\left(1|z|^2\right),\hfill & \hfill lm\\ z^{},\hfill & \hfill l=m1\\ z^{l1}\sqrt{1|z|^2},\hfill & \hfill m=1.\end{array}$$
(D15)
Substituting Eqs. (D15), (D12) into Eq. (D14) after straightforward calculations in the limit of short instantons, $`1|z|^2T/\mathrm{\Gamma }`$, for the matrix elements $`\delta F_{nm}`$ in the new basis we obtain

$$\frac{\delta F_{nm}}{T}=\frac{g_T}{2}\left[\frac{(n1)^2\delta _{n,m}}{n1+a}+\frac{a^2\left(z^{n+m}z^{nm}\right)}{(n1+a)(m1+a)}\right].$$
(D16)
Here $`n,m2`$ and $`a=\mathrm{\Gamma }/\pi T`$. For $`\mathrm{\Gamma }=0`$ the matrix $`\delta \widehat{F}`$ in Eq. (D16) becomes diagonal and independent of the instanton size, in agreement with Ref. . At finite $`\mathrm{\Gamma }`$ and in the limit of short instantons, $`1|z|^2T/\mathrm{\Gamma }`$, the off-diagonal matrix elements in (D16) are small and may be neglected for the purpose of evaluating the ratio of the two integrals over the fields $`\stackrel{~}{\varphi }`$ in Eq. (64). Proceeding this way we obtain the following expression for the ratio of the two integrals over the massive modes in Eq. (64) for $`1|z|^2T/\mathrm{\Gamma }`$

$$\frac{g_T^2E_C}{2\pi ^2T}\left[1\frac{\mathrm{\Gamma }}{\pi T}(1z^2)\mathrm{ln}\left(1+\frac{\mathrm{\Gamma }}{\pi T}\right)\right].$$
(D17)
Note that in the limit $`z1`$ this expression does not contain $`\mathrm{\Gamma }`$ and coincides with the result of Ref. .