Seven Etudes on dynamical Keldysh Model

We present a comprehensive pedagogical discussion of a family of models describing the propagation of a single particle in a multicomponent non-Markovian Gaussian random field. We report some exact results for single-particle Green's functions, self-energy, vertex part and T-matrix. These results are based on a closed form solution of the Dyson equation combined with the Ward identity. Analytical properties of the solution are discussed. Further we describe the combinatorics of the Feynman diagrams for the Green's function and the skeleton diagrams for the self-energy and vertex, using recurrence relations between the Taylor expansion coefficients of the self-energy. Asymptotically exact equations for the number of skeleton diagrams in the limit of large N are derived. Finally, we consider possible realizations of a multicomponent Gaussian random potential in quantum transport via complex quantum dot experiments.


P relude. Keldysh model
The original Keldysh Model (KM) was proposed by Keldysh in 1965 [1]. He showed that a single-particle problem of electron's propagation in a single-component random Gaussian field V (r) with forward scattering which is given by the correlator can be solved exactly by the Feynman diagrammatic technique. The idea behind Keldysh's solution is as follows. According to the diagrammatics rules, the diagram of order N contains N lines of the Gaussian random field (denoted by wavy lines), N + 2 Green's functions (denoted by solid lines) and 2N vertices [2]. The total number of diagrams is given by the number A N corresponding to the total number of possibilities to connect 2N vertices with N -lines: An example of the diagrammatic expansion corresponding to the single electron Green's function G is shown in Fig.1.
For the chosen Gaussian random field Eq.(1) all Feynman graphs in the order N give the same contribution to the Green's function (GF). As a result, the GF is represented by the following form: Here we use the shorthand notation E = − p 2 /(2m) and G 0 (E) = 1/E is the bare Green function. Using the representation for the Euler Gamma-function: and summing up the geometric progression results in a very simple expression for the full GF: This equation has a very simple physical meaning. Propagation of an electron in a Gaussian random field corresponds to an ensemble Gaussian averaging of the single particle Green's function over Gaussian realizations of the random potential. This is so-called zero-dimensional realization of the Gaussian disorder. The next very important contribution to the physics of the Keldysh model was made by Efros in [3]. He built a theory of semiconductors with very shallow charge impurities located at random positions. The potential created by the charge impurities is the screened Coulomb potential characterized by the screening radius r 0 . If the electron's density n 0 satisfies the condition n 0 r 3 0 1, then the electrons experience the potential only at the point at which they are located. The correlations of the impurity potential can be represented as a superposition of the short-scale and large-range parts. The large-scale fluctuations are exactly accounted in the Gaussian field approximation. The short-range correlations are accounted for by the perturbation theory. It was shown that the single particle Green's function satisfies the Dyson equation which for this model becomes an ordinary differential equation: The Dyson equation is obtained using the Ward identity [3], which has a very simple form due to δ(q)-shape of the Gaussian random potential correlator (1). One of the direct consequences of the Keldysh-Efros theory is the emergence of a tail in the single-particle density of states (DOS) at < 0. The formation of such tails in DOS is one of the general properties of disordered systems. For a detailed description of the original Keldysh model and its properties we refer to a remarkable book by M.V. Sadovskii [4].
Recently, it has been suggested that the Keldysh model can be realized not only in systems with long-range spatial variations of disordered potentials, but also in a system characterized by slow temporal fluctuations of the confining potential [5]. Such behaviour can be realized, for example, in functional nano-structures [6]. Generalized dynamical Keldysh models pave the way towards understanding of systems characterized by multi-component Gaussian non-Markovian random fields. A few-component dynamical Keldysh model have been proposed in [5][6][7][8] to describe effects of a slow electric gate noise in double quantum dots. In this paper, we develop the ideas of [5] in connection with implementation of the multi-component Keldysh model in the complex quantum dots.
These Notes are organized as follows: In Intermezzo we discuss the realizations of dynamical (time-dependent) multi-component Keldysh models in quantum circuits containing multiple quantum dots. We demonstrate that fluctuations of two potentials: the confinement potential and the inter-nodal barrier potential, which arise due to fluctuations of the corresponding electric gates, can create all the necessary ingredients for a multicomponent  The goal of this Section is to discuss possible realizations of time-dependent Keldysh models in functional nano-structure devices (see Fig. 2 for an illustration). The key elements of the electric circuit consist of quantum dots (see details below), gates controlling the shape of confining potential, gates controlling the shape of the inter-dot barriers and a global back gate operating uniformly on the dot's levels. The nano-devices are connected by the leads (contacts) to the rest of the quantum circuit.
blockade effects. The classical Gaussian random potential λ is defined by its mean value and the two-point correlation function: We use a color noise to describe the two-point correlator of the classical Gaussian noise.
where γ = 1/t noise defines the characteristic time for the noise correlation function and W is the amplitude of the noise. There are two important limits: Here D 1 (ω) is the Fourier transform of D 1 (t − t ). The equation corresponds to a fast noise (it's "fastest" realization is given by the white noise) and describes Markovian processes. The second equation corresponds to the a slow noise (it's "slowest" realization is known as the Keldysh model) representing the Gaussian disorder with a long (infinite) memory. InÉtudes we focus on a discussion of the physics described by the Keldysh model. The single particle Green's function in given realization of the quasi-static disorder is: Averaging over the realizations of the Gaussian disorder can be done equivalently by means of the distribution function P (λ) = 1/ √ 2πW 2 exp(−λ 2 /(2W 2 ).

Double QD
The single-particle Hamiltonian describing a two-level system (double well potential, double quantum dot Fig. 2b) consists of two terms: the diagonal H Z and the off-diagonal (tunneling) part H tun . The diagonal part is defined as follows: where α = l, r stand for the left and right quantum well respectively, S z = n l − n r . Thus, the diagonal part consists of a uniform term acting on the levels of each dot similar to that discussed in the previous subsection and an alternating term acting similarly to the longitudinal  random synthetic magnetic field (see Fig.4 upper panel) 1 . Since the effects of the uniform noise λ have been discussed above, in this subsection we focus on the effects of the synthetic magnetic field. The correlator of the longitudinal random Gaussian potential is given by: The off-diagonal (tunneling) Hamiltonian accounting for the tunneling through the barrier separating the two quantum wells (see Fig.4 lower panel) is given by: where ∆ = ∆ 1 + i∆ 2 . We introduced the pseudo-spin operators S + = c † r c l , S − = c † l c r such that the tunneling from one quantum well to another represents the pseudo-spin flip process. The two-component Gaussian random field describing slow non-Markovian fluctuations of the barrier is defined as: with Thus, the Hamiltonian describing the double well potential is expressed in terms of the basis of the Pauli matrices -the fundamental representation of the group SU (2) 2 The single particle Green's function in a given realization of the off-diagonal disorder (only) is: The lines |∆| 2 = ∆ 2 1 + ∆ 2 2 = const define circles S 1 in the d = 2 parametric space. 1 There is no real magnetic field in the model since all fields are associated with noisy electric gates. 2 The sum of the Hamiltonians (13) and (15) can be written as where τ i are the Pauli matrices (i = x, y, z) and τ 0 is the unit 2 × 2 matrix.
The single particle Green's functions in the presence of both the diagonal and the offdiagonal disorders are given by: The surfaces |∆| 2 + ∆ 2 3 = ∆ 2 1 + ∆ 2 2 + ∆ 2 3 = const define spheres S 2 in the d = 3 parametric space. When both uniform and staggered diagonal noises are present in addition to the offdiagonal noise, the GF in the given realization of quasi-static (infinite memory) Gaussian random potential casts the form:

Triple QD
Triple quantum dots (TQD) consist of three electron gas islands and either two (serial TQD, Fig. 2c) or three (triangular TQD, Fig. 2d) barriers separating the QDs. For the purposes of ourÉtudes we consider a serial TQD only to avoid a discussion of magnetic flux effects. The triangular TQD will be considered elsewhere. A cartoon illustrating the triple well potentials is shown in Fig. 5. A natural basis for representing single-electron processes in the serial triple quantum dots is the basis of eight Gell-Mann matrices -the fundamental representation of the SU (3) group (see Appendix A). The classification of all possible contributions to the Hamiltonian is shown on Fig. 5. Note, that we use a "rotated" representation of the Gell-Mann matrices [11][12][13] different from the conventional SU (3) basis. This basis embeds three SU (2), S = 1 generators as the first three Gell-Mann matrices (see Appendix A). The last five Gell-Mann matrices represent quadrupole moments expressed as bilinear combinations of S x , S y and S z [11][12][13][14].
First, there are two types of the non-uniform diagonal processes described by two diagonal Gell-Mann matrices µ 3 and µ 8 (see Appendix A) 3 forming the Cartan basis: We introduce two scalar random Gaussian fields η 1 and η 2 to account for the diagonal noise (see first two rows in Fig. 5).
The off-diagonal random processes can be classified as follows: i) Processes which can be considered as independent tunneling processes in the double well potential. These terms describe slow symmetric (in-phase) fluctuations of the two barriers (see the third row in Fig. 5): Here , and we label the states in the dot c α with α = l, c, r for the left, central and right quantum well correspondingly. ii) Processes which can be considered as anti-symmetric (out-of-phase) fluctuations of the barriers (see the fourth row in Fig. 5).
. We defined the two-component random Gaussian field as ∆ 34 iii) Processes directly connecting the left and right quantum wells (see fifth and the last row in Fig. 5).
HereP ± = (P x ± iP y )/2, and P + = c † r c l , (P + ) † = P − (see Appendix A). We define the two-component random Gaussian field as ∆ 56 = ∆ 5 + i∆ 6 . For a linear TQD such processes require co-tunneling and therefore are proportional to square of the tunneling matrix element (and therefore, smaller compared to the off-diagonal processes described by H 1 and H 2 ). For triangular TQD however, each dot is connected to its two tunnel partners and therefore there is no additional smallness in H 3 tunneling compared to H 1 and H 2 .
Summarizing, the five gates acting on the TQD (three for the quantum wells and two for the barriers) create eight random Gaussian fields (two single component random field for diagonal processes and three two-component random fields for the off-diagonal processes). The uniform noise gate acting simultaneously on all levels (not shown on the Fig. 5) can be added straightforwardly by an extra random field component.
The most general form of the single-particle Hamiltonian is: where h is the eight component Gaussian random potential (synthetic magnetic field) andˆ M is the eight-component vector consisting of eight Gell-Mann matrices. These arguments can be repeated for the serial complex Quantum dots consisting of M gated quantum wells separated by M −1 fluctuating barriers. Using the basis of SU (M ) group one straightforwardly obtains the Hamiltonian describing a particle in (M 2 − 1) -component random synthetic magnetic field.
InÉtudes we analyse the disorder-averaged Green's functionḠ( ) for Keldysh models with a multi-component Gaussian random potential. Since all models discussed in this manuscript are 0 + 1 dimensional with zero stands for spatial dimensions, DOS ν( ) is related to the single particle Green's functions as follows 4 : The single particle T-matrix is expressed in terms of the self-energy Σ( ) = G −1 0 −Ḡ −1 ( ): Knowing the exact form ofḠ R gives access to all physical single particle characteristics of QDs.
3Étude No 1. Keldysh model with one-component noise We start our firstÉtude by considering the original Keldysh model in the time domain [5]. The model describes, for example, single electron states in a single well potential (Single Quantum Dot) under the action of a fluctuating back gate (see the the Intermezzo). The "scalar" single component classical Gaussian random field (noise) is associated with the slow fluctuation of the confining potential. For the purposes of our notes we assume that the characteristic time associated with the change of the confining potential is greater than any other characteristic time of the system (e.g. Zener time, dwell time etc). Since we are dealing with single particle physics, statistics do not play any role. We will consider Fermi and Bose many-body models in a separate publication. One of the goals of ourÉtudes is to demonstrate the properties of the model and the response functions. To achieve this we will allow ourselves to show sketches of derivation and providing some minimal details of our calculations. The single particle Green's function (GF) is defined as follows 5 : Here the index (subscript) "1" refers to the one-component model. The inverse bare GF is G −1 0 ( ) = and Σ 1 ( ) is the self-energy part. As already mentioned, the classical Gaussian noise correlator is is quasi-static as the characteristic noise time t noise → ∞. Parameter W 2 is the variance of the Gaussian distribution. The self-energy (SE) part is defined by the following equation The first and second order Feynman diagrams for the self-energy are shown on Figure 6 a-c. SE is connected to the "triangular" vertex Γ 1 (Fig. 6 d-e). The Ward Identity (WI) establishes an explicit relation between Γ 1 and Σ 1 as the "transmitted" through the wavy line frequency is zero.
Combining the Dyson equation (31 -33) with the Ward Identity (34) we get a closed set of linear first order ordinary differential equations for GF (35) and a non-linear first order Riccati-type differential equation (36) for SE. 6 Boundary condition determines the large frequency behaviour of GF. It is convenient to introduce the following dimensionless variables and functions: As a result, the dimensionless equations for GF and SE are accompanied by the boundary conditions: The non-linear Riccati equation (40) can be transferred to a linear second order ordinary differential equation by substitution Corresponding boundary condition for the new function Φ 1 (z) is The equation for Φ 1 (z) reads as follows: or, equivalently, d dz Using the boundary condition Eq.(43) we re-write the equations in the following form: and therefore Φ 1 (z) ≡ Ψ 1 (z). Finally, the Ward Identity casts the form or, equivalently, Obviously, these equations can be obtained directly from the definition of σ 1 (z), the Ward identity and Eq.(39).

Analytic properties
This subsection is devoted to a discussion of the analytic properties of the single particle Green's function. The retarded dimensionless Green's function Ψ R 1 (z) can be obtained with the help of the Weierstrass transformation: Substituting it to the Dyson equation we get: The imaginary part ImΨ R 1 (z) can be easily computed: The real part ReΨ R 1 (z) can be determined by the Kramers -Kronig relations Here P.V. denotes Cauchy principle value of the integral. The above transformation is connected to the Hilbert transform (HT) of the Gaussian function H[e −u 2 ](z) = H 0 (z) [15] as shown in Appendix B: Evaluating Eq. (53) at z = 0 we obtain the obvious property ReΨ R 1 (z = 0) = 0 which can be used as a new initial condition for the linear ordinary differential equation for the function The solution of this equation is obtained by the method of indefinite coefficients and finally reads as: Implementing the initial condition ψ 1 (0) = 0 we obtain C = 0. Finally, ReΨ R 1 (z) can be expressed in terms of the Dawson function D + (z) = F (z) [16]: where erfi(z) is imaginary error function [17,18]. 7 As, the Green's function Ψ1 satisfies a linear differential equation, both real and imaginary parts of Ψ1 also satisfy linear differential equation. The equation for φ1(z) = ImΨ R 1 (z) reads as One can represent the Dawson function as the sin-Fourier-Laplace transform of the Gaussian function: The Dawson function can be expanded as follows for small values of the argument z → 0: and in the limit z → ∞ the expansion reads: Keeping first few terms of the Dawson function for small/large values of its argument, we get: Obviously, since GF is purely imaginary in the zero frequency limit Ψ R We plot the frequency dependence of the ReΨ 1 and ImΨ 1 in the left panel of Fig.7 Substituting the Taylor expansion for Ψ R 1 (z) into the equation for the self-energy we obtain the asymptotic behaviour of the self-energy at z → 0 In the z → ∞ limit The frequency dependencies of Reσ 1 and Imσ 1 are shown in the right panel of Fig.7. The DOS for the single component Keldysh model is represented by a single Gaussian-shaped peak centred at z = 0.
8 The real-time GF (dimensionless time s = t · W ) is defined through the Fourier transformation: Here Θ(s) is a step (Heaviside) function. The real time GF satisfies the differential equation One can notice a duality of Dyson equations in the time and frequency domain. This duality is associated with the Gaussian form of the random field.

4Étude No 2. Keldysh model with two-component noise
The two-component time-dependent Keldysh model describes, for example, the single electron states in a double well potential under action of slow fluctuations of the barrier separating the wells (see the Intermezzo) [5]. The complex tunnel matrix element has two independent components representing two Gaussian random fields (noise) in the x-and y-pseudospin directions.
The single particle Green's function of the two-component time dependent Keldysh model is defined as follows where G −1 0 ( ) = is the inverse barer GF. The two-component Gaussian random field correlator describes the two-component slow noise. For simplicity, we assume that two independent amplitudes of the noise coincide. Feynman diagrams for the Green's function before Gaussian averaging are shown on Fig. 8. At each vertex of the diagram the "pseudo-spin" is flipped and therefore black and white vertices form the staggered pattern. Dashed lines denote the Gaussian random potential. The self-energy part is proportional to the product of the bold GF G 2 and the triangular vertex Several first irreducible diagrams for the self-energy part are shown in Fig. 9. Each "black" vertex is connected only to a "white" vertex and vice-versa. As a result, there exists no second order irreducible diagram 9 Before we derive the differential equation for G 2 we show that GF can be obtained by direct re-summation of the Feynman diagrams of the "formal" perturbation theory where the coefficients A (2) n = (2n)!! reflect the fact that all diagrams at a given order of the perturbation theory have the same value. To proceed with the re-summation we use the integral representation for the Euler's Gamma-function. As the next step, we compute explicitly the sum over n to transform the integral to the form As a result, the Green's function is expressed in terms of the double integrals: 9 The reducible second order diagram and other reducible diagrams in all orders result in replacement of the bare GF in diagram Fig. 9 to bold GF. Transforming this integral from teh "Cartesian" coordinates to the "polar" coordinate system one gets: Due to the SO(2) symmetry of the noise the integral is independent on the polar angle. Here we would like to note that in fact the Eq. (75) describes an averaging of GF with the known Rayleigh distribution. Finally, by differentiating Eq. (75) with respect to the frequency we obtain a linear differential equation for the single particle Green's function. 10 The same equation can be derived with the help of the Ward identity The corresponding Dyson equation for the self-energy functions reads as follows: The boundary condition reads: Similarly to the firstÉtude we introduce dimensionless variables and functions: As a result, the dimensionless equations cast the form: being accompanied by the corresponding boundary conditions: It can be seen that in contrast the the one-component Keldysh model, an additional 1/z "centrifugal" term appears in the Dyson equation. As we will see below, this term is responsible for level repulsion. The boundary condition for σ 2 follows from the Ward Identity and the boundary condition for Ψ 2 (z): Therefore we have The solution of Eq. 81 with the boundary condition Eq. (83) is given by The non-linear equation (82) for the SE is a Riccati-type equation. Substituting and taking into account the boundary condition Eq. (83) we get the corresponding boundary condition for the auxiliary function Φ 2 The second-order linear ordinary differential equation (ODE) for Φ 2 (z) reads as follows: Figure 9: Irreducible self-energy part diagram in the first a) and third b) orders of the perturbation theory. Second order diagram is topological impossible. c) first non-trivial diagram for the vertex.
which can be re-written in more compact form: Here we introduced a new function The corresponding boundary condition for R(z) reads as: The general solution of the equation (90) is R(z) = A/z and the boundary condition (92) defines the constant A = 1. Now we use the definition of R(z) to obtain a first order linear ODE for Φ 2 (z): we obtain the ODE for Θ 2 : with the boundary condition By comparing Eq. (94) with Eq. (81) we conclude that Θ 2 (z) ≡ Ψ 2 (z). In addition we obtain new differential identities: The connection between the functions σ 2 (z), Γ 2 (z) and Ψ 2 (z) is given by the equations identical to Eq. (48):

Analytic properties
The retarded dimensionless Green's function Ψ R 2 (z) is given by The imaginary part ImΨ R 2 (z) can be easily computed as: Figure 10: The identity ReΨ R 2 (z = 0) = 0 follows directly from the analytic properties of Ψ 2 . The function ψ 2 (z) = ReΨ R 2 (z) satisfies the ordinary linear differential equation: The homogeneous solution is given by: The particular solution is obtained by the method of indefinite coefficients which leads to another differential equation for the function C(z): Finally, the solution of the boundary problem is given by: Here Ei(u) is the exponential integral function [17,18]. The real part ReΨ R 2 (z) is expressed through the Hilbert transform H 1 [ue −u 2 ] (see details in Appendix B): where γ is the Euler constant, |Ψ 2 (z → 0)| 2 ∝ z 2 log 2 [z]. Note, the ImΨ R 2 (z → 0) ∝ |z| is not differentiable at z = 0. Furthermore, ReΨ R 2 (z → 0) ∝ z log |z| is not analytic in the limit z → 0 and have a branch cut. So GF cannot be decompose into a Taylor expansion in the limit z → 0. The frequency dependencies of the real and imaginary parts of GF Ψ 2 are shown in the left panel of Fig. 10 11 . 11 Note that the real-time GF obtained by sin -Fourier transform can be written in terms of the Dawson function: The self-energy is diverging and non-analytic in the limit z → 0: The asymptotic behaviour of GF in the limit z → ∞ is: and . The frequency dependencies of the real and imaginary parts of the self-energy σ 2 are shown on the right panel of Fig. 10.
The DOS for the two -component Keldysh model is represented by two Gaussian peaks separated by a pseudo-gap with linear energy dependence in the vicinity z = 0.

5Étude No 3. Keldysh model with three-component noise
We start thisÉtude by considering a Double Quantum Dot realization of the Keldysh model. We assume that there exists a random Gaussian field proportional to τ x , τ y Pauli matrices, which is associated with the barrier fluctuations, and an additional alternating diagonal random Gaussian field acting on the levels in each dot. The last field is proportional to τ z . A realization of the alternating diagonal noise is done by applying two separate back gates acting to each QD. All three generators of SU(2) group are contained in the Hamiltonian. Without losing generality, we consider the case of equal noise amplitudes similarly to previousÉtude. Breaking the rotational symmetry is discussed in V ariation N o. 1. Another realization of the three component random field can be implemented when barrier fluctuations are accompanied by a uniform random Gaussian field acting on the levels of each QD through a single back gate applied to both QD in the double well potential. This noise is proportional to the unit matrix τ 0 . We discuss it in V ariation N o. 2.
The general equations relating the dimensionless SE σ 3 (z) and the dimensionless triangular vertex Γ 3 (z) with the single particle GF Ψ 3 (z) are: These equations are obtained directly from the definition of σ 3 (z) and the Ward Identity The differential equation for the Ψ 3 (z) and σ 3 (z) and the boundary conditions are given by Figure 11: The solution of the linear ODE for Ψ 3 (z) casts the form: Here Ω 3 = π/2 is the normalizing coefficient. Equivalently, the solution can be written in terms of the Hilbert transform H[u 2 e −u 2 /2 ](z) (see Appendix B) as: Less trivial identities are obtained by solving a non-linear Riccati-type equation for σ 3 (z).

Analytic properties
The retarded dimensionless Green's function Ψ R 3 (z) is given by the equation The imaginary part ImΨ R 3 (z) is easily computed The real part of the analytic in the upper half-plane of z function ReΨ R 3 (z) can be determined by the Kramers -Kronig relations: It can be also expressed by the Hilbert transform H[u 2 e −u 2 ](z) = H 2 (z) [15] (see Appendix B) which, in turn, is re-written in terms of the Dawson function as follows: Substituting H 2 (z) to the definition of ReΨ R 3 (z) we get: Using the properties of the Dawson function for small values of its argument, we obtain |Ψ R 3 (z → 0)| 2 ∝ z 2 . The frequency dependence of the real and imaginary parts of the GF Ψ 3 is shown in the left panel of Fig. 11.
Substituting the Taylor expansion for Ψ R 3 (z) into the equation for the self-energy we obtain the real and imaginary parts of the self-energy: The frequency dependencies of the real and imaginary parts of the SE σ 3 are shown on the right panel of Fig. 10.
The DOS for the three-component Keldysh model is represented by two Gaussian peaks separated by a pseudo-gap with the energy-square dependence at the vicinity z = 0. V ariation No. 1: Anisotropic two + one Keldysh model: Breaking rotational symmetry.
In this V ariation we assume that the amplitudes of the random complex Gaussian field associated with the fluctuations of the inter-dot barrier (which we denote below by W ⊥ ) and the amplitude of the random scalar Gaussian field associated with the fluctuations of the level position in each well of the double well potential, denoted by W are different. In this case it is not possible any longer to introduce the dimensionless functions of the dimensionless argument with a unique re-scaling parameter. For convenience of the reader we will not rescale the variables in this Subsection. The simplest way to calculate the Green function of the anisotropic three parametric Keldysh model is to use the path integral representation by introducing a generating functional (see details in [6]). In the extreme non-Markovian case corresponding to the "infinite memory" limit [5], the Green function of the anisotropic three-component Keldysh model (which we denote G 2+1 (ε)) is expressed as follows: Re-writing the integral in the spherical coordinate system and performing the integration over the angles we obtain after some simplification: Here erf(z) is the error function. [17,18] Let us analyse a particular case of anisotropy, namely W ⊥ > W . 12 The constant energy surface for the three-component anisotropic Keldysh model in both cases is an ellipsoid which transforms into a sphere if W ⊥ = W and isotropy is restored. Two extreme cases of the two-component Keldysh model W → 0 of the easy-plane anisotropy and one-component Keldysh model W ⊥ → 0 of the easy axis anisotropy can be obtained straightforwardly. To proceed with remaining integration over ρ it is convenient to use the series expansion of the error function As a result, Eq. (128) is transformed into Here α = (W 2 ⊥ − W 2 )/2, and G 0 ( ) = ( + iδ) −1 . The coefficients are the combinatorial coefficients in the n-th order of the "three-color" diagrammatic expansion for the full anisotropic three-component Keldysh model. In the isotropic limit α = 0 and only the term m = 0 survives. The coefficient C n0 = A n = (2n + 1)!!.

V ariation No. 2: Combined one + two Keldysh model
Let us assume that in addition to the slow barrier fluctuations in the double QD there exists a uniform noise to both QDs, which is produced by back gate. For simplicity we suppose that the noise amplitudes for the one-component uniform diagonal Gaussian random potential and the two-component off-diagonal Gaussian random potential are the same. The dimensionless GF Ψ 1+2 (z) is defined as: which is, in fact, the Weierstrass transform of the two-component Keldysh model Green's function.
Green's function Ψ 1+2 (z) satisfies the second order ODE which can be obtained by differentiating Ψ 1+2 (z) given by Eq. (132) with respect to it's argument, integrating by part and using the differential equation for Ψ 2 (z) Eq. (81): The corresponding boundary conditions follow from fixing two coefficients in the z → ∞ asymptotic: In fact, these boundary conditions (134) correspond to two conditions: The analytic properties of Ψ 1+2 (z) are fully defined by the corresponding properties of Ψ 2 (z) (see previousÉtude).
These equations can be obtained directly from the definition of σ d (z) and the Ward Identity Here we assumed that the GF has only "radial" component in d-dimensional spherical coordinate system. The Dyson equations for the GF and SE are accompanied by the corresponding boundary conditions: Here is the normalizing coefficient and Γ(u) is the Euler's Gamma-function. The (d − 1)/z "centrifugal" term is responsible for the level repulsion.
The differential identities connecting GF and SE are given by: These identities follow from the solution of Riccati-type equation for σ d (z).

Large D-limit
Differential equations describing the single-particle Green's function of the Keldysh model with a few-component Gaussian random potential contain coefficients O(1). Keldysh models with a D + 1 1-component random potential contain a natural small parameter 1/D which can be used to find an approximate but reliable solution based on the 1/D expansion. To obtain the differential equation in the form convenient for the expansion we first proceed with re-scaling of the functions and the arguments: The differential equations for new GF and SE read: The Ward Identity acquires the form: It is straightforward now to find an approximate solution of these equations. In particular, one can see that the boundary condition gives a reliable approximation for both the GF and SE at z > D: The imaginary part of the Ψ R D (z) in this approximation consists of two δ-functions separated by ∆z = 2 √ D. To summarize, adding components of the random Gaussian field allows to engineer the gap in the DOS: few-component Keldysh models are characterized by a "soft" (pseudo) gap which becomes "harder" when the number of components increases.

7Étude No 5. Perturbative expansion for the bold Green's function
The three last advancedÉtudes are devoted to combinatorics of the Feynman diagrams for the bold Green's function and derivation of recurrence relations allowing to find a number of skeleton Feynman diagrams in a given order of the perturbation expansion. The "skeleton" means that we are looking for irreducible diagrams containing the bold GF (see e.g. Fig  12). As it has been mentioned before, the advantage of slow non-Markovian classical random potential models is that all diagrams in a given order of the perturbation theory have the same value. Thus, the coefficients of the expansion give the total number of the diagrams. 13É tude 6 is based on the works of Sadovskii -Kuchinskii [4,9] and Suslov [10] (SKS) who considered the enumeration of skeleton diagrams for the "orthodox" single-component Keldysh model. In thisÉtude we present a simple pedagogical derivation of the SKS recurrence equation using a method different from the original SKS works. This method proves to be very convinient for deriving a general equation giving access to enumeration of skeleton diagrams for classical random Gaussian field models with arbitrary number of the random field components. These equations will be derived and discussed inÉtude 7. We will show that the SKS equations is just a particular limit of the general theory. In this preparatoryÉtude we create a basis for answering questions about the combinatorics of Feynman diagrams. The first question is: "How many Feynman diagrams exist in n-th order of the perturbative expansion?". Or, put more simply "How many Feynman diagrams contain n-wavy lines for d-color random Gaussian field models?" This question is trivial as we already know the answers for one-, two-and three-component models. To answer this question we need to extend the general solution for the d-color GF Ψ d in terms of the bare GF in the limit z → ∞: Using equation for Ω d = 2 d 2 −1 Γ( d 2 ) from the previousÉtude we get This coefficient defines the number of Feynman diagrams in the n-th order of GF expansion.
To check that the general answer is correct, we write down explicitly the values of A (d) n for d = 1, 2, 3: Another important check is the number of diagrams in the lowest (first) order of the perturbative expansion: 13 The situation is completely different for fast (Markovian) random potentials. In this case, the selection rule says that the rainbow-type diagrams give much large contribution in comparison to the diagrams with self-crossing. Finally, for the limit d 1: The numbers given by Eq. (146) being plugged into the equations for the self-energies and the triangular vertex give the answer about the corresponding combinatorics in terms of bare GF.
8Étude No 6. Enumeration of "skeleton" Feynman diagrams in a single-component Keldysh model: square recursion Summation and re-summation methods of Feynman diagrams are widely used in analytical (quantum field theories) and numerical (bold diagrammatic Quantum Monte-Carlo methods, see, e.g. [27][28][29]) solutions in various problems of modern many-body physics. Counting the number of Feynman diagrams is another challenging problem, which is important, for example, for benchmarking numerical techniques based on the diagrammatic expansion. Hedin in his pioneering work [19] formulated a set of equations (Hedin's equations) connecting the single particle propagator G with the effective potential W , self-energy Σ, polarization Π and vertex Γ. Another additional equation follows from the functional-derivative connection between Γ, Σ and G. In zero dimension of the space-time the Hedin equations become algebraic. By searching the solution of the Hedin equations in powers of the bare interacting potential and taking into account the number of closed loops, one gets the number of Feynman diagrams from the expansion coefficients [20][21][22]. Further development of these methods made it possible to formulate an iterative algorithm for counting Feynman diagrams via many-body relations [23] . In thisÉtude we develop analytical methods for counting skeleton Feynman diagrams for self-energy, triangular vertex and single-particle T-matrix for models with d-component classical (no loops) Gaussian random field. We are dealing with 0 + 1 models and therefore will end up with some differential (not algebraic, contrast to the Hedin's theory) equations. The methods discussed inÉtudes 6 and 7 can also be generalized to enumerate diagrams for the two-particle Green's function and the full four-point correlator. In order to find the combinatorics for the self-energies and the triangular vertices of the single-component Keldysh model in terms of the bold GF we start with the general equations Eq.(141) formally defining σ 1 and Γ 1 as Taylor series in terms of the bold GF Ψ 1 : Thus, the goal is to find z[Ψ 1 ]. Since Ψ 1 (z) is expressed in terms of the Dawson function, we aim to construct the inverse Dawson function [16]. The coefficients a (1) n in the asymptotic series fully determine the number of irreducible Feynman diagrams for the self-energy and the vertex correspondingly.
Using equation (150) we write: We outline a simple and elegant way to construct the inverse Dawson function below. The function z[ which is obtained by inversion of Eq. (39). The boundary condition for this equation reads Substituting expansion for z[Ψ 1 ] into (153) we get: which after opening the brackets is re-written as follows: Equating the coefficients in front of the n-th power of the bold GF Ψ n 1 we obtain the recurrence relations for the coefficients a (1) n : The equation (162) can be written in more compact form 14 The self energy σ1 satisfies Abel's second kind non-linear differential equation [24] σ1 · dσ1 dψ = σ1 where ψ=Ψ1 is the bold GF and the initial condition is defined by behaviour of σ1 (67) and Ψ1 (62) at z → 0 [9]. Strictly speaking as both σ1 and Ψ1 are complex functions of a parameter z (seeÉtude 1), the Eq. (154) represents a system of two coupled non-linear partial differential equations. We consider Eq. (154) assuming ψ→x=Reψ and s= Reσ1 while Imψ→ 0 and Imσ1→0. It corresponds to parametric regime z 1 where both ImΨ1/(ReΨ1)→ 0 and Imσ1/(Reσ1) → 0. The Eq. (154) is transferred to a normal form: Introducing a new function u(ξ) = −1/s allows to convert the non-linear ODE Eq. (155) to Abel's first kind non-linear differential equation [26]: To see the connection with Sadovskii -Kuchinskii (SK) theory [9] we redefine σ1(ψ) = ψ · Q(ψ). As a result, substituting σ1(ψ) into Eq. (154) we obtain the non-linear differential equation for the function Q which was derived by a different method by Sadovskii and Kuchinskii and analysed in [9]: The solution of this equation formally written in terms of the Taylor series with zero convergence radius defines the recurrence relations (163).
The number of the irreducible diagrams in the first ten orders of the perturbation theory is presented in the Table 1.
The recurrence relations for the diagram combinatorics of the single-component Keldysh model have been obtained by Sadovskii and Kuchinskii [9] and Suslov [10] using a bit different technique. The large-n asymptotic of the coefficients a (1) n has been analysed by Sadovskii and Kuchinskii in [9]. In the limit n 1 SK transferred the recurrence equation (163) into an ordinary differential equation. Its solution gives the asymptotic: We plot in Fig. 13 (left panel) the normalized number of irreducible Feynman diagrams for the first n = 50 orders of the perturbation theory (dots). The red solid line shows asymptotic behaviour given by Eq.164 accounting for the 1/n correction. The SK method, while powerful enough to handle 1/n corrections, does not, however, allow the coefficient (1/e) before (2n + 1)!!! to be determined analytically. This coefficient was calculated numerically by SK and derived analytically by Suslov [10]. In the next and the lastÉtude we will derive the recurrence relations for the self-energy and the triangular vertex combinatorics considering the multi-component Keldysh model. We will derive general recurrence equations which give the number of the skeleton diagrams in any order of expansion. To illustrate how to get the asymptotic behaviour of the coefficients of d-component models at the limit n 1 we apply Sadovskii and Kuchinskii method. ThisÉtude begins with derivation of the recurrence equation for the two-component Keldysh model. Following the key steps of the derivation we present detailed discussion for computing the large -n asymptotic for the number of skeleton diagrams. After "learning" how the method works for a two-component random Gaussian field we proceed with the general derivation for the multi-component case. At the end of thisÉtude we briefly discuss how the combinatorics of skeleton diagrams for self-energy can be used for the combinatorics of the diagrams for the T -matrix.
Enumeration of "skeleton" diagrams for the two-component Keldysh model The differential equation for the inverted function z[Ψ 2 ] is obtained by inversion the differential equation for Ψ 2 (z): We substitute the Taylor expansion for  to get the recurrence relations for the coefficients a (2) n : This equation finally leads to the following recurrence relations: We notice an interesting property of the case d = 2, namely, that the coefficient a 1 = 0. It reflects the fact that the second order irreducible Feynman diagram vanishes due to topology of the two-component model (seeÉtude 2). The equation (167) can be further simplified: This recurrence equation generalizes the combinatorics of the Keldysh model to the case of the two-component random Gaussian field. The main difference from SKS equation is that in addition to the bilinear recurrence, a cubic term appears in the recurrence equation.
To find the asymptotic behaviour of a (2) n at n 1 we will follow the method developed by Sadovskii and Kuchinskii in [9]. We look for the asymptotic solution of (167) in the form a (2) n ≈ (αn + β) a (2) n−1 .
Substituting Eq. (167) into Eq. (167) we find that at the limit n 1 We can therefore re-parametrize Eq. (167) by introducing the new coefficients as follows: This re-parametrization gives the asymptotic behaviour for the number of the skeleton diagrams for the self-energy in the large-n limit. In order to find the 1/n correction we need to deal with the equations for the coefficients c n−m−p−2 .
Taking the limit of large n and retaining only the leading 1/n 2 terms in (170) under assumption c = c n ≈ c n−1 ≈ c n−2 we finally obtain a simple differential equation: The solution of this equation is given by As it has been discussed before, Sadovskii-Kuchinskii method does not provide information about the coefficient in front of the exponential function. To obtain correctly this coefficient we follow the logic of Suslov's work [10]. Let's begin with d = 1 case. For all d = 1, since the noise components are independent (no cross-correlations), we can repeat the same arguments for each component independently and hence the constant found for the d = 1 case will be replaced by the constant to the power d. First, we observe that the coefficient B n−1 has the same n-dependence as the coefficient A (1) n which defines the total number of the Feynman diagrams for GF 15 . Second, we see that the diagrams contributing to A n /2! and so on. As a result, 16 we get the relation: Similarly, for d = 2 we obtain Finally, using const = 1/e 2 we get: Note that for d = 2 each random Gaussian field line gives the factor 2 because the correlator τ + τ − = τ x τ x + τ y τ y and the correlations in the x-and y-pseudospin directions are statistically independent. The second column of the Table 1 contains coefficient a (2) N −1 computed for the first ten orders of the perturbation theory. The number of irreducible diagrams in the n-th order is given by the coefficient b N −1 /2 N . Note that the coefficient b (2) 2 = 0 since the second order irreducible Feynman diagram is absent (seeÉtude 2).
We plot in Fig. 13 (right panel) the normalized number of irreducible Feynman diagrams for the first n = 50 orders of the perturbation theory (dots). The red solid line shows the asymptotic behaviour given by Eq. (174) accounting for the 1/n correction. 15 It is easy to see that the number of skeleton diagrams in the n-th order of the perturbation theory is defined strictly speaking by the coefficients B we obtain recurrence relations between the coefficients a (d) Finally, following the steps outlined in the previousÉtude we get: The last equation can be re-written as follows: This equation is the central point of thisÉtude and the central result of the paper. The large-n asymptotes and the 1/n expansion for the coefficients a (d) n can be straightforwardly performed using the Sadovskii-Kuchinskii-Suslov [4,9,10] method which leads to the following final result: where the function f (d) is computed in the previous subsection for two particular cases: f (d = 1) = 5/4 and f (d = 2) = 1. We leave the problem of determining f (d) for arbitrary d as an exercise for readers.

Enumeration of "skeleton" diagrams for the T-matrix
The T-matrix is related to the self-energy by the following relation: T d · G 0 = Σ · G. Equivalently, for the dimensionless T-matrix function of the dimensionless variable z, Combinatorics of the diagrams for the T-matrix is defined through the equation: Therefore, there is no need to construct independent equations for the T-matrix combinatorics since it can be determined from the number given by the coefficients a (d) n : with Thus, the combinatorics of the coefficients a (d) n for the self-energy fully determines the combinatorics of the vertex part and the single particle T-matrix.

Coda. Concluding remarks
To summarize, let us discuss what we can learn from playing SevenÉtudes on dynamical Keldysh model.
• The single particle Green's function of the Keldysh models with a d-component non-Markovian classical random Gaussian field satisfies the first order ordinary linear differential equation. The Dyson equations can be obtained in the closed form due to the Ward Identities and solved exactly for an arbitrary number of the random field components.
• The exact solution for the single-particle Green's function corresponds to an ensemble Gaussian averaging of the single particle Green's function over the Gaussian realizations of the random potential.
• The imaginary part of the single particle GF (proportional to the Density of States) has the Gaussian shape.
• If the number of components is d > 1, there exists a "centrifugal" term in the Dyson equation responsible for level repulsion. The level repulsion becomes more pronounced when the number component of the random fields grows.
• All Keldysh models with even numbers of components contain intrinsic non-analyticity in the single particle correlators in the zero frequency limit due to the Rayleigh-type distribution. The single particle correlators of the Keldysh models with odd numbers of random field components are analytic at all frequencies due to the Gaussian distribution.
• The total number of skeleton Feynman diagrams in a given order of expansion with respect to a multi-component random potential satisfy some recurrence equations containing linear, bilinear and cubic terms. The single component Keldysh model has some additional degeneracy as both linear and cubic terms are absent.
• The combinatorics of the skeleton Feynman diagrams for the self-energy completely determines the corresponding combinatorics of the full vertices and the single-particle T-matrices.
Multi-component Keldysh models describe the behaviour of complex (double, triple-etc) quantum dots subject to slowly fluctuating electric potentials which determine the shape of the quantum wells and inter-dot barriers. The application of the Keldysh model to computing dynamical correlation functions of many-body Fermi and Bose systems, as well as some other realizations of the Keldysh model in theories with multi-component synthetic gauge fields [30][31][32][33][34][35], polaron problem [36][37][38] and Sachdev-Ye-Kitaev [39] model as well as applications to Weyl semimetals [40][41][42] opens interesting directions for future research.

Acknowledgements
The Authors are thankful to late Kostya Kikoin for many years of fruitful collaboration on numerous topics of modern condensed matter physics including but not limited by dynamical symmetries and the Keldysh model. We are grateful to Boris Altshuler for very fruitful discussions on disordered systems and Nikolay Prokof'ev and Boris Svistunov for teaching the basics of the bold diagrammatic Quantum Monte Carlo method. DE thanks DFG (project number 405940956) for the partial financial support. The work of MK is conducted within the framework of the Trieste Institute for Theoretical Quantum Technologies (TQT).

Gell-Mann matrices
The Gell-Mann matrices constitute the fundamental representation of the group SU (3). The eight generators satisfy the commutation relations: We refer to the conventional form of the Gell-Mann matrices as the λ-basis, Λ α = λ α /2.: The λ-matrices are normalized by the identity T r(λ α λ β ) = 2δ αβ .
The structure constants are It is convenient to use the rotated basis of the Gell-Mann matrices [12]. The rotation is performed to embed the representation S = 1 of the group SU (2) as first three generators of the group SU (3). We refer to the rotated representation of Gell-Mann matrices as the µbasis, M α = µ α /2: The µ-matrices are normalized by identity T r(µ α µ β ) = 2δ αβ .
The structure constants are