{'fraction_non_alphanumeric': 0.042957719965362746, 'fraction_numerical': 0.013120740935958736, 'mean_word_length': 4.304343484772841, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relaxation dynamics of the electron distribution in the Coulomb blockade problem

We study the relaxation dynamics of electron distribution function on the island of a single electron transistor. We focus on the 'interaction without coherence regime in which an electron coherence can be neglected but quantum fluctuations of charge are strong due to Coulomb interaction. The quantum kinetic equation governing the evolution of the electron distribution function due to escape of electrons to the reservoirs is derived. Analytical solutions for time-dependence of the electron distribution are obtained in the regimes of weak and strong Coulomb blockade. We find that usual exponential in time relaxation is strongly modified due to the presence of Coulomb interaction.


I. INTRODUCTION
Recently, the phenomenon of Coulomb blockade in single electron devices 1-5 has come into the focus of the field of thermoelectricity. 6 Among significant experimental achievements one can list development of the Coulomb blockade thermometer, 6 the thermal rectifier on the basis of a quantum dot, 7 new technique to measure temperature gradients across a quantum dot, 8 etc. The standard characteristic of thermoelectric performance is the figure of merit which involves the product of conductance, thermopower squared and inverse thermal conductance. Measurements of thermopower and thermal conductance in single electron transistors (SET) and quantum dots have been performed during the last decade at different temperature regimes. [9][10][11] The theory of the thermoelectric effects in electron devices has been put forward in Refs. [12,13]. During the last decade the thermopower and thermal conductance have been studied in single electron transistors and quantum dots, [14][15][16][17][18][19][20] and in granular metals [21][22][23][24][25] in various regimes. However, the thermopower and thermal conductance are linear response parameters and, therefore, describe the equilibrium properties of a system only.
In contrast, our work is focused on properties of single electron devices in the out-of-equilibrium regime which has attracted a lot of theoretical interest recently. In particular, the conductance of a quantum dot under ac pumping in the stationary non-equilibrium state was obtained in Ref. [26], the current noise of an ac-biased quantum dot was studied in Ref. [27], the non-equilibrium dephasing rate and zero-bias anomaly in single electron transistor was computed in Ref. [28], the statistics of temperature and current fluctuations in the fully outof equilibrium single electron transistor was investigated in Refs. [29,30], and the extension of the P (E)-theory 31 to the out-of-equilibrium regime has been developed in Refs. [32,33]. However, these works dealt with regimes when the distribution function of electrons in a quantum dot or an island of single electron transistor were fixed by external sources, e.g., ac or dc bias voltage.
In this paper we address a different question: how does an electron distribution function once being prepared relaxes toward the equilibrium state in the Coulomb blockade problem. Apart from general physical interest in understanding of a non-equilibrium regime, the answer to this question is important for the field of electron thermometry. 6 We consider the simplest system: single electron transistor. The set-up is shown schematically in Fig. 1. Metallic island is coupled to an equilibrium electron reservoirs via tunneling junctions. Depending on the task, the reservoirs and the island may be kept at different temperatures (T l , T r respectively), different chemical potentials (constant or varying in time: µ l (t), µ r (t)) or quasi-stationary gate voltage (U g (t) = U 0 + U ω cos ωt) may be applied to the system. The physics of the system is governed by several energy scales: the Thouless energy of an island E Th , the charging energy E c , and the mean single-particle level spacing δ. Throughout the paper the Thouless energy is considered to be the largest scale in the problem. This allows us to treat the metallic island as a zero dimensional object with vanishing internal resistance. The dimensionless total conductance (in units e 2 /h) of tunneling junctions g is an essential control parameter. The junctions are assumed to have a large number of channels but the conductance of each one is assumed to be small g (l,r) ch ≪ 1. The temperature is assumed to be low enough: T ≪ max{1, g}E c in order to keep electrons strongly correlated due to Coulomb interaction. At low temperatures the interplay of Coulomb interaction and electron coherence dominates the physics of single electron devices. Account of both effects is a formidable undertaking. Therefore, we are going to restrict ourselves to the regime of not very low temperatures in which an electron coherence can be neglected but quantum fluctuations of charge are strong due to Coulomb interaction. Then, the physics of the system is Island Lead Lead g l g r V l , T l V r , T r T d adequately described in the framework of Ambegaokar-Eckern-Schön (AES) effective action. 34 The effective action describes the system in terms of bosonic field ϕ, which is usually termed as the plasmon field. Its time derivative is interpreted as fluctuating electric potential of electrons inside the island. AES-approach has wellknown limitations. Deriving AES action one assumes that the products of electron Green's functions averaged over disorder are substituted with products of disorderaveraged Green's functions in every calculation. That is why the processes of phase-coherent multiple impurity scattering inside the island are left out. The limitations in the regime g ≫ 1 and g ≪ 1 were discussed in detail in Refs. [35] and [36], respectively. It was shown that at temperatures T ≫ max{1, g}δ, AES-action approach is justified. Following Ref. [37], we shall term the temperature range, max{1, g}E c ≫ T ≫ max{1, g}δ as an interaction without coherence regime. This 'interaction without coherence' regime is an attainable experimental reality, e.g. in experiments reported in Refs. [9] and [10] the necesary conditions were satisfied.
In the case of strong Coulomb blockade (g ≪ 1), the theoretical study of relaxation of an electron distribution in the interaction without coherence regime has been done before for a single quantum dot [38] and for an 1D array of quantum dots [39]. However, the considerations of Ref. [38] have been restricted by assumptions that i) the electron distribution is the Fermi function with some temperature different from the equilibrium one; ii) transport is dominated by co-tunneling processes (Coulomb valley regime); iii) temperatures of the island and the reservoirs are close to each other.
In the present paper, we undertake the analysis of relaxation of an electron distribution function which is free of above-mentioned restrictions.
Since we are going to capture non-equlibrium physics, we employ the formalism of AES-action in its out-ofequilibrium form throughout the paper. We supplement it with quantum kinetic equation to explore relaxation dynamics of electron distribution. For a SET with large number of tunneling channels we derive the quantum kinetic equation with the collision integral due to escape of electrons to the reservoirs. It is valid in the entire span of values of g and generalizes the one obtained in Ref. [26] for sequential tunneling (first order in g) and cotunneling (second order in g) approximations in the framework of the orthodox theory of the Coulomb blockade. In fact, our collision integral is always an infinite series in powers of g. Indeed, each tunneling event is accompanied by the radiation of a plasmon. That is why the collision integral becomes of the infinite order in the distribution function of electrons inside the island. This situation is entirely different from the one in Fermi liquid and leads to non-trivial relaxation.
As a test of the quantum kinetic equation, in the regime of linear response we derive analytical expressions for transport coefficients: conductance, thermal conductance and the response of electric current to temperature difference. In the regime of weak Coulomb blockade (g ≫ 1) we establish the following new results for the transport coefficients: i) the conductance and thermal conductance violate Wiedemann-Franz law, and deviation of the Lorentz ratio L from value π 2 /3e 2 demonstrates weak periodic dependence on the gate voltage; ii) the thermopower weakly oscillates with the gate voltage around zero value. Weak oscillations of the Lorentz ratio and thermopower with the gate voltage found in the regime g ≫ 1 are manifestation of the known gatevoltage dependence of these quantities [12][13][14][15][16][17][18][19][20] in the strong Coulomb blockade regime, g ≪ 1.
In weak and strong Coulomb blockade regimes we have employed the quantum kinetic equation to solve the relaxation of the electron distribution in two cases: i) the distribution of electrons inside the island is the Fermifunction with some temperature; ii) the distribution function of electrons inside the island is arbitrary. In the former case we have managed to extract the relaxation dynamics of the electron temperature; in the latter case we have obtained evolution of a distribution function itself. In both cases we assumed that electron escape to reservoirs is the primary relaxation mechanism. In general, the collision integral in the quantum kinetic equation is non-local in energy due to inelastic nature of tunneling processes: the radiation of plasmon always accompanies the tunneling event. In a number of wide parametric regimes: weak Coulomb blockade and Coulomb peak in the strong Coulomb blockade, the kernel of the quantum kinetic equation acquires a quasi-elastic form. However, the collision integral remains non-local in energy due to renormalization effects in these cases. The co-tunneling regime is qualitatively different: the kernel of the collision integral is entirely inelastic.
Our new result is that despite quasi-elastic form of the collision integral, strong Coulomb interaction dramatically changes the relaxation laws comparing to simple exponential ones expected from golden-rule type arguments. They suggest that electron relaxation rate is to be proportional to the width of electrons' levels inside the island, gδ, prompting simple exponential relaxation. The renormalization effects due to Coulomb interaction make the width of electrons' levels dependent on the electron distribution and lead to the non-exponential relaxation laws. For example, in the regime of the sequential tunneling, we have discovered that there is a time regime in which relaxation of the electron temperature in a SET island is independent of the tunneling conductance g.
The paper is organized as follows. In Sec. II we introduce the hamiltonian and essential parameters of the problem. Sec. III is devoted to the out-of-equilibrium AES-model and to derivation of the quantum kinetic equation. Sec. IV is devoted to derivation of general expressions for the linear response coefficients. The relaxation dynamics of electrons in the island is explored in the weak (g ≫ 1) and strong (g ≪ 1) Coulomb blockade regimes in Sec. V-VI. Discussion of the results, comparision with other relaxation mechanisms, different from electron escape to reservoirs and conclusions are presented in Sec. VII.

A SET is described by the Hamiltonian
where describes free electrons in the leads and the island, H c describes Coulomb interaction of carriers in the island, and H t describes the tunneling. Here operators a create a carrier in the i-th lead (island). Then, the tunneling hamiltonian is The charging Hamiltonian of electrons in the box is taken in the capacitive form: Here E c = e 2 /(2C) denotes the charging energy, andn d is an operator of a particle number in the island: To characterize the tunneling it is convenient to introduce the following hermitean matrices: The first of them acting in the Hilbert space of the states of the lead, the second -in the space of the islands states.
The energies ε α are accounted from the Fermi level, and the delta-functions should be smoothed on the scale δE, such that max{δ, δ (l,r) } ≪ δE ≪ T . Here, δ and δ (l,r) stand for mean level spacing of single-particle states on the island and reservoirs, respectively. The classical dimensionless conductance (in units e 2 /h) of the junction between a reservoir and the island can be expressed as follows 5 Therefore, each non-zero eigenvalue ofĝ (i) orǧ (i) corresponds to the transmittance of some 'transport' channel between a reservoir and the island. 40 The effective number of these 'transport' channels (N (i) ch ) is given by The effective dimensionless conductance g ch of a 'transport' channel can be written as follows The dimensionless conductance g then becomes In what follows we will always assume Notice that under these circumstances the conductances g l,r can still be large provided the effective number of channels N (l,r) ch is sufficiently large.

III. ACTION AND KINETIC EQUAITONS
A. AES-action To tackle the system which is out of equilibrium we have to employ essentially non-equilibrium formalism. Keldysh technique is thus the only way through. We employ Keldysh form of AES-action (we sketch the known details of derivation in Appendix A): 1,41 where Here ϕ c,q = (ϕ + ± ϕ − )/2 with ϕ ± denoting bosonic field on both branches of Keldysh contour. Physically, the bosonic field is associated with the fluctuating electric potential on the island. In terms of classic and quantum boson exponents the dissipative part of AES-action reads: Here Π R,A,K are corresponding components of electron polarization operator in the Keldysh space. They are given by a standard formulae presented for reference in Appendix A. In a case of constant density of states (DOS) in the island and leads the kernel of the AES action can be simplified: Here we define a slow time τ = (t + t ′ )/2. Function F ε (τ ) is given in terms of the Wigner transform f ε (τ ) of the electron distribution function f (t, t ′ ): F ε (τ ) = 1−2f ε (τ ). Π R,A,K ω (τ ) are Wigner transforms of corresponding functions in time domain.
As seen from the structure of the r.h.s. of (19), it is suitable to introduce function F r ε = ( α g α F α ε )/g which may be called the effective distribution function of two reservoirs. It is that combination of reservoir distribution functions that enters all the following equations of the paper.
Although Eq. (18) is exact we neglect all derivatives with respect to slow time τ in Eq. (19). It is also convenient to introduce function B ω (τ ) in accordance with relating Keldysh and retarded (advanced) components of polarization operator. The function B ω (τ ) plays a role of a distributiion function for electron-hole excitations. In the equilibrium it is given by coth(ω/2T ).
In what follows we assume that electrons in the leads are locally thermalized such that f l,r ε (τ ) are the Fermifunctions. Depending on the parameters of the model, both quasi-equilibrium and non-equilibrium regimes can exist in the island. Therefore, we assume F d ε (τ ) to be an arbitrary function slow varying with time τ , and derive the kinetic equation for F d ε (τ ) by which the AES action should be supplemented.
is an averaged single particle density of states in the island and Σ K,R,A are the components of self-energy in Keldysh space. To the second-order in tunneling Hamiltonian H T (lowest order in 1/N ch ) the Wigner transform of the self-energy shown in Fig. 2a Here we perform Wigner transform of the exact selenergies and introduce the correlation functions of boson exponents: It is convenient to parametrize D K (t, t ′ ) via the boson distribution function B(t, t ′ ): It is worthwhile to mention that the next ( fourth order in H T ) contribution to the self-energy which is shown in Fig. 2b is of the order g 2 /N ch . This correction to the selfenergy is of the same order as terms omitted in the course of derivation of the AES action (16). In the considered limit N ch ≫ 1 it can be safely neglected.
Performing Wigner transform of Eq. (21) and neglecting all slow time derivatives in its r.h.s. we obtain the quantum kinetic equation for the distribution function of electrons on the island of a SET: This quantum kinetic equation constitutes one of the main results of the present paper. It describes evolution of the distribution function F d ε (τ ) of electrons in the island due to interaction with boson field ϕ and tunneling to the leads and back. The quantum kinetic equation (25) is derived for any values of g r and g l ; the r.h.s. of Eq. (25) can be written as the series in powers of g due to the presence of Im D R ω (τ ) and B ω . The boson distribution B ω is determined by electron distribution function F d ε (τ ) and should be found from the solution of the AES-action (13).
At g ≫ 1 the kernel ( Im D R ω ) of the collision integral of the quantum kinetic equation (25) resembles the kernel of the collision integral in the quantum kinetic equation for disordered electron liquid [43][44][45] for energy transfers ω ≫ gδ as it is expected. 35 At g ≪ 1 the quantum kinetic equation (25) which takes into account the renormalization effects via Im D R ω generalizes the kinetic equation derived in Ref. [26] in the framework of the orthodox theory 46 for sequential tunneling and inelastic cotunneling approximations.

IV. TRANSPORT COEFFICIENTS
Using quantum kinetic equation (25) we are able to derive general formulae for all linear response coefficients of the SET for any value of g. Voltage (∆V = V r − V l ) and temperature (∆T = T r − T l ) differences across the SET cause charge (I (e) ) and heat (I (q) ) currents. Electric and thermoelectric transport coefficients are defined as Here coefficients M and G T (the response of a heat current to voltage difference and the response of electric current to temperature difference, respectively) are related via Onsager relation M = G T T . 47 The thermal conductance κ is usually defined as κ = K − G V T S 2 where S = G T /G V stands for the thermopower. The electric and heat currents in the α-th reservoir can be found as The current conservation corresponds to the condition I e l + I e r = 0. It fixes the boson distribution function B ω to be equal to the electron-hole distribution function B ω introduced in Eq. (20): . (28) The heat current conservation I q l + I q r = 0 determines the equilibrium temperature of the island: A straightforward computation of charge and heat currents gives Introducing the quantities g ′ , g ′ T , and k ′ as we obtain  We stress that Eq. (32) is valid for any value of tunneling conductance g. It generalizes expressions for transport coefficients obtained in Refs. [20,48] for g ≪ 1 to arbitrary values of g.
It is worthwhile to express Eq. (32) in terms of the tunneling density of states of electrons in the island (see Appendix C): Substituting expression (33) for the tunneling density of states and performing standard integrals with Fermi and Bose distribution functions one can check that the results (32) are allowed to be exactly rewritten in the form  Eq. (34) for the transport coefficients resembles the corresponding expression in the Fermi-liquid. 47,49 However, contrary to Fermi-liquid, the tunneling density of states ν d (ε) has strong dependence on energy for ε → 0. In general, ν d (ε) = ν even is even/odd function of ε. It can be shown that ν even/odd d (ε) is even/odd function of the external charge q. Therefore, g ′ and k ′ are even functions of q whereas g ′ T is an odd function of the external charge.
For macroscopic samples of ordinary metals, the Wiedemann-Franz law provides a universal relation between the conductance and thermal conductance. It states that the Lorenz ratio L = κ/(G V T ) , is a constant given by the Lorenz number π 2 /3e 2 . As follows from Eq. (34), one can expect the violation of Wiedemann-Franz law in the presence of strong dependence of ν d (ε) on electron energy.
In the case g ≫ 1 one is able to perform perturbative expansion in 1/g and take into account non-perturbative corrections. The function Im D R ω acquires the following form in the equilibrium 50 Here, function ωδ(ω) can be understood as Im a/[π(ω + a+ i0)] where the limit a → 0 should be performed at the very end of all calculations, (this calculation can be, e.g. integration over ω). The non-perturbative in 1/g corrections (exponential terms exp(−g/2)) come from Korshunov instantons 51 of the AES-action. Then by using Eq. (35) we find from Eq. (32) The result for g ′ has been obtained in Ref. [37]. Equations (37)- (38) are new and valid for temperatures T ≫ g 2 E c exp(−g/2). We emphasize that g ′ T has only nonperturbative instanton) contribution. The same holds for the thermopower: At g ≫ 1 the violation of the Wiedemann-Franz law is weak and the Lorentz number is given as Due to the presence of the non-perturbative contribution, the Lorentz number is temperature dependent and oscillates as a function of the external charge q. Eqs. (39) and (40) constitute one of the main results of the present paper.
In the strong coupling regime g ≪ 1, Eq. (32) supplemented by the proper expression for Im D R ω (cf. Eqs. (75) and (114)) results in exactly the same expressions for the transport coefficients as obtained in Refs. [17,19,20,48]. We refer a reader to these works for details.
Next, we want to illustrate the ability of quantum kinetic equation (25) combined with fine field-theoretical scaling of essential physical quantities. We consider the problem of relaxation of electrons in the island towards the equilibrium due to the tunneling to the reservoirs and back. There are two possible scenarios. The first one can be refered to as a quasi-equilibrium regime. The electron distribution inside the island is given by the Fermifunction but with non-equilibrium temperature T d which slowly relaxes to its equilibrium value. The second scenario is fully non-equilibrium regime when electron distribution is arbitrary. Which scenario persists depends on the ratio τ E /τ ee where τ E stands for the energy relaxation time due to tunneling mechanism and τ ee for the energy relaxation time due to electron-electron interaction in the island. The non-equilibrium regime persists provided τ E ≪ τ ee and the quasi-equilibrium regime is possible if τ E ≫ τ ee . We will argue below (see Sec. VII) that both scenarios are possible.
There is one more relaxation time involved: τ RC which determines relaxation of the electric charge on the island. In the weak Coulomb blockade regime, τ RC is given by the following clasical estimate: τ RC ≃ 2π/gE c . As we shall see below, τ E ≫ τ RC . Therefore, it is allowed to assume that at first there is quick relaxation of the electric charge on the island and, then, slow relaxation of the electron distribution function or temperature towards the equilibrium. Technically, it means that initial elec- As was discussed in the Introduction, the renormalization of physical observables drastically changes the relaxation dynamics of the system. Therefore, before solving kinetic equation we need to establish the scaling of a theory's coupling constants under non-equilibrium conditions.
A. Renormalization of AES-action at g ≫ 1 The AES-action is renormalized due to its nonlinear form. In the equilibrium case renormalization of the action is well-known (see e.g., Ref. [1]). In our case non-equilibrium makes the problem non-trivial. As in equilibrium, we expect the necessary scaling of the coupling constant g. The additional question that inevitably arises is whether the structure of the kernel of AES action (the components of polarization operator Π ω in Keldysh space in Eq. (16)) is changed due to renormalization? The details of the calculation are presented in Appendix D. We prove that the structure of the bare action is fully restored, the kernel of the AES action being intact during renormalization group (RG) procedure. The coupling constant renormalizes according to Here high energy scale Λ is naturally set by the first term in Eq. (14): Λ ∼ gE c . To demonstrate that the integral in Eq. (41) is indeed logarithmic we explore the behavior of the integrand at ω → ∞. It is straightforward to get the following asymptotic for function B ω at ω → ∞: We expect that any physical distribution function obeys the condition F d ε → sgn ε at ε → ∞. Then Therefore, the high-energy asymptotic of function B ω is given by sgn ε as in the equilibrium. This way, the logarithmic behavior of integral in Eq. (41) is ascertained.
To get the renormalized action one has to integrate out all the frequencies down to the lowest scale ω 0 , at which the RG stops.This energy scale can be determined as Let ε d be a characteristic energy scale of the island distribution function (the scale at which electron distribution function F d ε becomes almost equal to sgn ε). Then one can easily check that the following estimate holds (see Appendix D for elaborate details) where T r , T l are temperatures of the reservoirs. Energy scale ω 0 serves as a natural lower cut-off, Λ = ω 0 , in the RG procedure. Finally, we find In the equilibrium, ε d = T r = T l = T and one finds ω 0 = T . Eqs. (45)-(46) describe renormalization of the AES-action under non-equilibrium conditions.

B. Non-equilibrium regime
The relaxation problem is formulated as follows. At t = 0 the island is heated and some electron distribution function F d ε (0) is created. The characteristic energy ε d of electrons in the island is larger than temperatures of the reservoirs, which are kept fixed and equal to each other ε d > T r = T l . The system is released and the island is cooling down due to the tunneling of electrons to the reservoirs and back.
Performing expansion to the second order in boson fields ϕ c,q , one straightforwardly finds (see Appendix D) We mention that this result generalizes the perturbative (independent of q) part of Eq. (35) to the nonequilibrium case. With the help of (47) one can compute the collision integral in the r.h.s. of the quantum kinetic equation (25) and obtain Here we neglect last term in Eq. (47) for the following reasons. It gives contribution of the order of unity whereas the first term in Eq. (47) involves dω B ω /ω ∼ ln gE c /ε d ≫ 1. Although, Eq. (48) has a quasi-elastic form, in fact, it is highly non-linear equation: G(τ ) involves information about the electron distribution at all energies. As was shown above, the quantity G(τ ) has meaning of the renormalized coupling constant of the theory. Simple algebra leads us to the differential equation for the function G(τ ): 52 The solution reads Now by using this result we integrate Eq. (48) and obtain the evolution of the electron distribution function F d ε : Eq. (51) demonstrates energetically uniform relaxation of the electron distribution function. This fact is a direct consequence of the quasi-elastic form of the kinetic equation (48). However, due to renormalization effects the form of the relaxation law is different from the exponential one.
Let us define the characteristic energy ε d as ε 2 Then, in the case ε d (0) ≫ T r and at not too long times τ ≪ 2π/δG r , one finds from Eq. (51) that the characteristic energy decreases according to the power-law: C. Quasi-equilibrium regime In the quasi-equilibrium regime, we need to take into account the collision integral I (ee) ε due to electronelectron interaction in the island. 43,44 As this term is added to the r.h.s. of Eq. (25), it makes the electron distribution F d ε to be the Fermi-function. Multiplying both parts of Eq. (25) by ε and integrating them over energy, we obtain the following equation (the well-known identity dε εI Here we use the leading (classical) part of Eq. (47) ( Im D R ω = −2πδ(ω)/B ω ). Equation (53) yields standard exponential relaxation towards the equilibrium. In the limit T d ≫ T r (G(0) ≫ G r ) it is also possible to compute collision integral using the entire one-loop expression (47) of its kernel. Naturally, one-loop correction reveals itself in the logarithmic renormalization of g in the r.h.s. of Eq. (53). By using Eq. (47), we find where c is a numerical constant of the order of unity which does not influence final results. The solution of (54) reads: The condition in the second line of Eq. (55) implies that solution holds for not too long times at which T d (τ ) ≫ T r (G(τ ) ≫ G r ). The logarithmic renormalization of the conductance changes the character of temperature relaxation. At long times 2π/δ ≪ τ ≪ (2π/δ) ln G(0)/G r , the cooling of the island slows down in comparison to the standard exponential decay which is developed at short times τ ≪ 2π/δ: It is instructive to compare the relaxation of temperature in the quasi-equilibrium regime and the characteristic energy ε d in the non-equilibrium regime given by Eqs. (56) and (52) at times τ ≪ 2π/δG r , respectively. While the former demonstrates exponential behavior, the latter decreases in accordance with the power-law.

VI. RELAXATION OF ELECTRONS IN THE ISLAND, STRONG COUPLING REGIME, g ≪ 1
In the strong coupling regime there are two possible scenarios for relaxation of electrons in the island of the SET. The first one persists if τ E ≪ τ ee . In this nonequilibrium case the carriers inside the island do not have time to thermalize and to form the Fermi-distribution with some temperature. In this case the the time evolution of distribution function itself becomes the main objective. This task is solved in section VI C below. The second scenario develops in the opposite limit, τ E ≫ τ ee . Namely, the relaxation rate due to electron-electron interaction inside the island is much faster than the rate due to electron tunneling through the contacts. Thus, the temperature of carriers in the island becomes a well defined characteristic of a system. Consequently, the relaxation of the island's temperature will be the focus of our analysis in section VI D.
As in the previous section we shall assume that the electric charge on the island quickly relaxes and only then slow relaxation of the electron distribution or temperature starts. In the strong Coulomb blockade regime this picture is well justified since 1/τ RC ≃ g max{T, ∆} ≫ 1/τ E .
We concentrate on the most interesting case: the vicinity of a degeneracy point: q = k + 1/2 where k is an integer. Following Ref. [53], the hamiltonian (1)-(4) can be simplified by truncating the Hilbert space of electrons on the island to two charging states: with Q = k and Q = k + 1 (see Fig. 3). The projected hamiltonian then takes a form of 2 × 2 matrix acting in the space of these two charging states. Denoting the deviation of the external charge from the degeneracy point by ∆: q = k + 1/2 − ∆/(2E c ) we write the projected hamiltonian as: 53 where H 0 is given by Eq.

A. Non-equilibrium pseudo-fermions
To deal with spin operators it is standard to use Abrikosov's pseudo-fermion technique. 54 We introduce two-component pseudo-fermion operators ψ † α , ψ α such that The out-of-equilibrium pseudo-fermions were tackled before. 55,56 As usual, one introduces Keldysh contour, doubling the number of fermions. The system is out of equilibrium and one has to be very cautious. The distribution function F ε of pseudo-fermions is not known a priori. Rather, it is to be defined self-consistently from corresponding kinetic equation. Pseudo-fermions are also subject to constraint on their number: Thus the state of a system ought to be projected on the state with N = 1 at any instant of time. The operator of particle number is conserved by Hamiltonian (57)- (58). Consequently, the operator of projection on to physical subspace N = 1 commutes with hamiltonian too. It means that the projection on to physical subspace is needed at a single point of Keldysh contour only. We insert the factor exp(η αψ α ψ α ) into density matrix and take the limit η → −∞ at the end of any diagrammatic calculation. Then Provided the operator O has zero expectation value in the sector with zero pseudo-fermion number, N = 0, Eq. (61) can be simplified as The dissipative action is to be rewritten in the Keldysh representation. We plug representation (59) into the hamiltonian (57) and integrate out electrons in the lead and the island. This leads to the following effective action Here σ stand for Pauli matrices, σ ± = (σ x ± iσ y )/2, and are matrices in Keldysh space. The pseudo-fermion operators ψ are understood as vectors in the tensor product of isospin and Keldysh space. Π ij stands for the matrix of polarization operator (16)- (17). Next, we write the Wigner transform of the quantum kinetic equation for pseudo-fermion distribution function Here as before, we neglected all derivatives with respect to slow time τ . All functions entering Eq. (65) are understood as matrices acting in the isospin space. From the appearance of Eq. (65) we conclude that characteristic relaxation time of pseudo-fermion distribution function F ε (τ ) is τ pf ∼ 1/(g max{∆, T }) and is much shorter than τ E . It allows us to consider pseudo-fermions to be in the stationary state. Then the l.h.s. of the kinetic equation (65) can be omitted and we obtain the equation for pseudo-fermion distribution function: With the help of Eq. (63) we write down equations for the pseudo-fermion self-energies (Fig. 7): Here G tt ′ ,σ stands for the pseudo-fermion Green functions corresponding to the first line in Eq. (63) and Σ ± are matrices in the Keldysh space. We will need the explicit expressions for their Wigner transforms: Here σ stands for ± and Combining Eqs. (66) and (68) we find the following equation for the pseudo-fermion distribution function where F σ = F σ ∆σ/2−η . Plugging ε = ∆σ/2 − η we arrive at the closed equation for F σ : Now we need to investigate asymptotic properties of functions F σ when pseudo-fermion chemical potential η → −∞. It is natural to expect that the equilibrium result survives in the non-equilibrium. As one can check this assumption satisfies Eq. (71). In order to solve the quantum kinetic equation (25), we need to compute Im D R ω in the strong coupling limit g ≪ 1. With the help of Eqs. (58) and (59), one easily finds in the zeroth order in g: Now we express the physical correlation function through pseudo fermion one Im D R ω,pf . By using the following zeroth order in g result for the pseudo-fermion number Eq. (75) is the generalization of the equilibrium result for correlation function Im D R ω (see Refs. [48,50]) over the non-equilibrium case.
Next, by using Eqs. (58) and (59), we find in the zeroth order in g: Expressing the physical correlation function through pseudo-fermion one D K ω,pf , we obtain This result implies that the boson distribution function B ω is determined by B ω in the same as in the weak coupling regime, Before proceeding with the solution of the quantum kinetic equation we prefer to perform one-loop renormalization of the theory. This is done to sum up all large logarithmic corrections (which otherwise arise in perturbative analysis) and absorb them into renormalized physical constants of the theory.

B. One-loop structure of the pseudo-fermion theory
In this section, we establish the out-of-equilibrium generalization of the scaling of fundamental parameters in the pseudo-fermion theory (the gap ∆, the coupling constant g), the Green's function and the average pseudofermion density N pf . We expect that the action (63) can be renormalized with only one scaling parameter Z like in the equilibrium case [57][58][59]. This is indeed the case and the obtained renormalized structure of the theory is a natural generalization of the equilibrium one. The renormalized pseudo-fermion Green's function becomes See Appendix E for details of the computation. It is straightforward to check that coupling constant g and gap ∆ are renormalized according to: To complete the renormalization picture we need to establish the scaling dimension of the pseudo-fermion number N pf . In complete analogy with Ref. [57], N pf happens to have no renormalization For completeness we present the rigorous proof of Eq. (82) via Callan-Symanzik equation in Appendix E. The integral in Eq. (80) runs over frequencies E c ≫ |ω| ≫ ω 0 = max{T r , ε d ,∆}. The energy scale ω 0 determines the natural scale at which the RG procedure has to be stopped. The Green's function (79) acquires the width whereF σ ≡ F σ ξσ . Therefore, the renormalized physical correlation function becomes

C. Electron distribution relaxation in the island
In this section we consider the relaxation in the nonequilibrium case, τ E ≪ τ ee . We focus on the most interesting case of the Coulomb peak: ∆ = 0. Then the quantum kinetic equation (25) is greatly simplified (cf. Eq. (48)): Here, we stress that the kinetic equation (85)-(86) is of the infinite order in the electron distribution function on the island. Indeed, λ involves F d ε via electron-hole distribution function B ω .
The formal solution reads The function G(τ ) obeys the differential equation The solution of Eq. (88) is given as f (z) = z + ln(1 − z).

By using the relation
which follows from Eq. (88), we obtain Since Eq. (89) can not be solved analytically with respect to G(τ ) it is instructive to investigate limiting cases.
Let us assume that the effective energy of electrons in the island ε d ≫ T r such that G(0) ≫ G r . Then, expanding f (z) in the series in z, we find Eq. (92) is valid provided G(τ ) ≫ G r , i.e., for not too long times: τ ≪ π/δG r . It is worthwhile to mention that standard exponential relaxation occurring at short time τ ≪ πG r /(δG 2 (0)) transforms into regime of slower relaxation at intermediate time πG r /(δG 2 (0)) ≪ τ ≪ π/(δG r ): At longer time τ ≫ π/(δG r ), function G(τ ) becomes almost equal to G r and we find again the regime of standard exponential relaxation: The same exponential relaxation as given by Eq. (95) holds if the effective energy of electrons in the island ε d is slightly larger than T r such that G(0) − G r ≪ G(0), G r .

D. Temperature relaxation in the island
Now we investigate the relaxation in the quasiequilibrium case, τ E ≫ τ ee .

Coulomb peak, ∆ = 0
We start from the regime of the Coulomb peak: ∆ = 0. In the quasi-equilibrium regime, one needs to add to the r.h.s. of Eq. (85) the collision integral I (ee) ε due to electron-electron interaction in the island. It is this term that makes the electron distribution to be a Fermifunction. By using the well-known identity dεεI (ee) ε = 0, we obtain the following equation: where G(τ ) is given by Eq. (86). In the quasi-equilibrium case, we can not derive closed equation for G(τ ) as it was done in the non-equilibrium case due to the presence of additional term I Assuming that T d (0) ≫ T r we can estimate λ with logarithmic accuracy as λ = ln E c /T d . Then, we find from Eq. (96) and The solution (98) is valid provided the condition T d (τ ) ≫ T r holds. If T d (0) ≫ T r exp(π 2 /G(0)), then the exponential relaxation developing during initial period τ ≪ 2π 3 /[δG 2 (0)] transforms into regime of slower relaxation at intermediate time: We mention that in this regime the temperature relaxation is independent of the quantity G(0) which determines the SET conductance. In the opposite case, T d (0) ≪ T r exp(π 2 /G(0)) the temperature T d (τ ) evolves according to Eq. (99) for τ ≪ (4π/δG(0)) ln T d (0)/T r . At longer times τ ≫ (4π/δG r ) ln T d (0)/T r the temperature T d (τ ) becomes of the order of T r : T d (τ ) − T r ≪ T r and we find the standard exponential relaxation: Evolution of the temperature of electrons in the island is presented in Fig. 4. We mention universality of the relaxation at long time when the difference between the electron distribution in the island and in the reservoirs becomes small. In non-equilibrium τ E ≪ τ ee and quasiequilibrium τ E ≫ τ ee regimes the relaxation is exponential with a rate of the order of δG r . The same exponential relaxation as in Eq. (101) holds if the temperature of electrons in the island T d (0) is slightly larger than T r , It is worthwhile to mention that there is a parametric region of time domain G r /(δG(0) 2 ) ≪ τ ≪ 1/(δG(0) 2 ), when the relaxation of the distribution function in the non-equilibrium regime is much slower ln τ than the relaxation of the (Fermi) distribution function in the quasi-equilibrium regime, i.e., relaxation of temperature, ln( . 4: The dynamics of temperature relaxation, g ≪ 1,

Coulomb valley,∆ ≫ T d (0)
Now we consider the relaxation of the electron temperature on the island in the regime of Coulomb valley, ∆ ≫ T d (0). By using Eq. (84), we rewrite the quantum kinetic equation (25) as We remind that we consider the quasi-equilibrium regime. Then we need to add to the r.h.s. of Eq. (102) the collision integral I (ee) ε which describes scattering due to electron-electron interaction in the island. In what follows we assume that the condition T d (0) ≫ T r holds. With the help of the following results which are valid for T r ≪ T d (li 2 (z) = ∞ n=1 z n /n 2 denotes the polylogarithmic function), we obtain from Eq. (102) We can estimate parameter λ with logarithmic accuracy and find λ = ln E c /∆ since the temperature of electrons in the island T d ≪∆. Therefore, bothḡ and∆ are independent of τ . Integration of Eq. (102) yields Here Ei(z) = − ∞ −z dt exp(−t)/t stands for the integral exponential. By using the asymptotic h(z) = exp(z)/z 2 at z ≫ 1, we obtain and The results (109) and (110) are valid at not too long times As expected, due to the exponentially small SET conductance in the sequential tunneling regime, the temperature relaxation is very slow, namely, logarithmical. Therefore, it is instructive to consider contribution to the temperature relaxation due to the electron co-tunneling.

Inelastic cotunneling regime
As known very well, due to exponential suppression of the sequential tunneling mechanism deep in the Coulomb valley, T d ≪∆, the higher order process of inelastic cotunneling dominates the transport. 60 Contrary to the case of sequential tunneling, the cotunneling contribution to the collision integral in the r.h.s. of Eq. (25) comes from frequencies of order ω ∼ T d ≪∆.
In the pseudo-fermion technique the inelastic cotunneling is revealed as the broadening of delta-peaks in the imaginary part of the retarded and advanced pseudofermion Green functions. 48,50 After taking into account Eq. (83), the integrand in Eq. (73) becomes of a complex pole structure. There are two pairs of proximal poles There is an additional series of Matsubara-type poles resulting from distribution functions F + ǫ and F − ǫ+ω . They lead to logarithmically divergent sums. The latter are controlled by the renormalization scheme. In our case all leading logarithms are absent. They have already been absorbed into renormalized constantsḡ and∆ by the proper choice of reference energy scale. Thus we can omit all divergent sums over Matsubara frequencies. Expanding in the ω/∆ we obtain Next we use the same arguments that led us to leading order espression (75). The function Im D R ω then reads Using Eq. (114), we rewrite the quantum kinetic equation (25) as We remind that we consider the quasi-equilibrium regime. We mention that Eq. (115) coincides with the kinetic equation derived for the cotunneling regime in Ref. [26]. Then we need to add to the r.h.s. of Eq. (115) the collision integral I (ee) ε which describes scattering due to electron-electron interaction in the island.
In the case of T d − T r ≪ T r , we obtain where G r =ḡ 2 T 2 r /(6∆ 2 ) stands for the equilibirum SET conductance in the cotunneling approximation. In the opposite case T d ≫ T r , by using Eq. (115), we find the following equations: which govern the temperature relaxation. It is worthwhile to mention that if one substitutes G r by G(τ ) in Eq. (116) then it becomes similar to Eq.(4) of Ref. [38] for V = 0 and in the absence of phonons. However, due to different numerical coefficients in the right hand side of Eqs. (116) and (117) and Eqs. (119) and (120) are valid at times (121)

VII. DISCUSSIONS AND CONCLUSIONS
We have studied the relaxation dynamics of the SET under essentially non-equilibrium conditions. The language of kinetic equations happened to be the most adequate for this task. Analytical results are procured in the limiting cases of weak g ≫ 1 and strong g ≪ 1 Coulomb blockade. All relaxation equations (see Eqs. (48),(85),(96),(105)) obtained in the course reveal a pleasant generality. Namely, Here X d is a relaxing physical quantity (temperature, distribution function), and G(X d ) is conductance of a SET which depends on X d . Equation (122) has a transparent intuitive interpretation. Namely, the characteristic time scale determined by the r.h.s. of Eq. (122) is simply a dwell time of a particle inside the metallic island, 27 i.e. τ −1 E ∼ Gδ. The inverse dwell time can be also estimated as the ratio of the thermal conductance κ to the heat capacitance of the island. The latter is proportional to T d /δ. The generality of Eq. (122) is, however, deceptive as it leads to drastically different evolution of physical quantities over time in the case of small and large values of g.
In the course of all our analysis we generally discarded the influence of electron-phonon (e-ph) interaction. The reasoning behind this is as follows. The e-ph scattering rate was well studied for 2-dimensional electron gas with disorder. 61 The following estimate has been found The electron-electron (e-e) scattering rate in mesoscopic systems is widely studied as well (see, e.g. [62]). For small diffusive electron systems and for T ≪ E Th there are two parametrically different situations 62,63 where L D = (k F l) inverse screening length.The upper one comes from momenta k ∼ 1/L and of diffusive origin. Let us address the question which kind of dissipation dominates in various parametric regimes. The relaxation due to electron tunneling can be roughly estimated as By comparing Eqs. (124), (125) and (126), one can see that both quasi-equilibrium and non-equilibrium regimes can occur for g ≫ 1 and g ≪ 1. The non-equilibrium regime prevails for g ≫ 1 while the quasi-equilibrium one dominates for g ≪ 1 (see Figs. 5 and 6).
To estimate e-e scattering rate we use some experimental data taken from the experiment by Pasquer et al. 64 where they studied the Coulomb blockade effects in a small island of two-dimensional electron gas. The experimental data were as follows: the level spacing δ ≈ 85 mK, the Fermi energy E F ≈ 47 K, the elastic mean free path l ≈ 15 µm, the size of an island L ≈ 1 µm. This allows us to estimate the Thouless energy E Th ≈ 8 K and e-e relaxation rate The typical temperature of the contemporary mesoscopic experiment is T 100 mK. As we see, with lowering temperature the e-ph scattering rate decays faster then the corresponding electron-electron (e-e) rate. On the other hand for the same metallic island the typical relaxation rate due to electron escape to reservoirs is Estimates (123)-(128) show that for all relevant experimental temperatures the phonons are frozen and e-ph interaction can safely be discarded. Next, comparing estimates (127) and (128) we conclude that varying g two different parametric regimes explored in this paper can indeed be realized in the experiment. Namely, fully non-equilibrium regime is realized when g is large enough and electron distribution function is arbitrary inside the island. The quasi-equilibrium regime persists in the opposite limit, when g is small enough.
In addition to the relaxation of the electron distribution in the island due to escape of electrons to the reservoirs which we consider in details above there is another mechanism of energy relaxation which is due to interaction U ir of electrons in the island with electrons in the reservoirs. For a sake of simplicity we assume that the typical interaction parameter r s ∼ 1/(k F a B ) ∼ 1 with a B standing for Bohr radius. In the case L ≪ L D , the energy relaxation rate due to interaction U ir of electrons in the island with electrons in the reservoirs can be estimated as Here U ir (k) denotes the Fourier transform of the interaction U ir (r). Provided the condition k F d ≫ 1 holds U ir (k F ) is strongly suppressed, ν d U ir (k F ) ≪ 1, and 1 τ (ir) ee In the opposite case of large island, L ≫ L D , and for a B ≪ L the estimate for the energy relaxation rate 1/τ Here k L ∼ 1/L and U (k) stands for the Fourier transform of interaction between electrons in the island k = 0 component of which leads to the charging term H c in the hamiltonian (1). As one can see, both cases of τ ee ≫ τ For d ≪ L where d stands for the typical size of the tunneling junction between the island and reservoir Eq. (131) can be simplified as For the experiments by Pasquer et al. 64 we estimate the Bohr radius a B ≈ 10 nm and assume typical d to be of the order of 100 nm. Therefore, we expect that the regime in which the main mechanism of the energy relaxation of electrons in the island is due to its escape to the reservoirs can be realized in a laboratory.
To summarize, we have explored heat transport and relaxation processes in a SET with large number of tunneling channels over a wide range of parameters. In the regime of linear response we obtained analytical expressions for transport coefficients (conductance, thermal conductance and the response of electric current to temperature difference) in the entire span of values of g. It is possible to shape the general relations for linear response coefficients into Fermi-liquid type form. There is however an important difference, namely: the tunneling density of states undergoes dramatic renormalization due to Coulomb interaction. The latter leads to violation of Wiedemann-Franz law: in the g ≫ 1 limit the Lorentz ratio L acquires weak periodic dependence on gate voltage (the precursor of Coulomb blockade). The method of quantum kinetic equation supplemented with non-equilibrium AES action has allowed us to treat Coulomb interaction exactly. We have obtained the time evolution of electron temperature (in the quasiequilibrium regime) and the distribution function (in the non-equilibrium regime) of a SET island due to particle escape to reservoir. The corresponding collision integral is always non-local in energy due to inelastic nature of tunneling processes: the radiation of plasmon ϕ always accompanies the tunneling event. In general, this leads to highly complicated integro-differential kinetic equations. Surprisingly we have shown that kinetic equations can be reduced to simple differential ones in a number of wide parametric regimes, namely: g ≫ 1 (weakly blockaded SET) and g ≪ 1 (strongly blockaded SET in sequential tunneling approximation with renormalization taken into account). This simplification is achieved due to the presence of strong scale separation in the problem gδ ≪ T d or gδ ≪ ε d . Indeed, the characteristic frequency at which the distribution function in the kinetic equation changes is ω ∼ gδ, while the scale at which the renormalization due to the presence of Coulomb interaction occurs is ω T d or ω ε d . This separation is that allows us at first to treat Coulomb interaction and secondly to study evolution of the distribution function.
Still, quantum fluctuations of charge significantly change the relaxation laws comparing to simple exponential ones which are characteristic of semi-classical physics for g ≫ 1 and of orthodox theory for g ≪ 1. The regime g ≪ 1, ∆ ≫ T is dominated by cotunneling process.
In the latter case the kinetic equation retains its integrodifferential structure and is to be solved numerically elsewhere. Measurements of the predicted relaxation dynamics are an experimental challenge. Performing integrations over fast frequencies ω, ω 1 , ω 2 we see that the first integral is ∼ 1/((t − t ′ )Λ) and the second one is ∼ 1/((t − t 1 )(t ′ − t 2 )Λ 2 ). Thus they are irrelevant for RG-analysis. It means that only term S II contains logarithmic in Λ corrections. As usual we are interested in the first non-vanishing ϕ-dependent correction: Working out the trace in Eq. (D8) we obtain .

(D9)
Substituting it into (D2) we see, that the structure of the AES-action is restored. The only difference is the change of the coupling constant given by Eq. (41). Finally, we mention that in the case of non-zero ∆Q Eq. (45)should be changed to ω 0 ∼ max{ε d , T r , T l , |∆Q|δ}. The term with β g always contains extra g and can be dropped in the leading order.
To find γ we need to find N pf in the next to the tree-level order. The diagram representing the correc-tion to pseudo-fermion particle number is presented in Fig. 8. Calculating with logarithmic accuracy and using extensively the fact that B ω → sgn ω at large ω we obtain Plugging Eq. (E10)) into Eq. (E8)) we find that γ = 0 that proves Eq. (82).