Phases of translation-invariant systems out of equilibrium: Iterative Green's function techniques and renormalization group approaches

We introduce a method to evaluate the steady-state non-equilibrium Keldysh-Schwinger Green's functions for infinite systems subject to both an electric field and a coupling to reservoirs. The method we present exploits a physical quasi-translation invariance, where a shift by one unit cell leaves the physics invariant if all electronic energies are simultaneously shifted by the magnitude of the electric field. Our framework is straightaway applicable to diagrammatic many-body methods. We discuss two flagship applications, mean-field theories as well as a sophisticated second-order functional renormalization group approach. The latter allows us to push the renormalization-group characterization of phase transitions for lattice fermions into the out-of-equilibrium realm. We exemplify this by studying a model of spinless fermions, which in equilibrium exhibits a Berezinskii-Kosterlitz-Thouless phase transition.


I. INTRODUCTION
Unconventional phases of matter play an integral role in condensed matter research and beyond. 1 Understanding the conditions under which systems harboring many particles conspire to give rise to these emergent, collective phenomena is crucial from a fundamental as well as a technological perspective. The description of such phases also poses a formidable theoretical challenge as they are usually driven by interactions and independent particle pictures fail spectacularly. To remedy this, powerful many-body techniques such as the renormalization group where developed. After years of research, much is known about the classification of phases of matter in thermal equilibrium as well as about the transitions between them. 1 As a second step, one might wonder about ways of controlling these phases beyond the possibilities offered by equilibrium means. 2 A particular non-equilibrium route that is routinely followed in experiments is to apply electric fields to solids. If the electric field is strong enough, the linear response regime is left, electrons are driven out of equilibrium, and non-linear effects become relevant. This so-called non-linear transport regime has attracted much interest in the last decades and many counter-intuitive effects were demonstrated. E.g., it was shown that a negative differential conductance 3-6 and oscillating currents (thyristor effect) 7 can arise and that this might have significant implications for highly efficient heat engines. 8,9 For very strong electric fields compared to the scattering rate of electrons, coherent Bloch oscillations are found, 10,11 which in the absence of scattering will not decay. 12 In a metallic condensed matter setup, these oscillations are challenging to observe experimentally because the scattering-induced relaxation is usually very fast on the time scale of the oscillation frequency and the steady-state current quickly relaxes to zero. [13][14][15] However, these oscillations can be accessed in semiconductors 16 or cold-atom systems. 17,18 For noninteracting electrons (i.e., in the absence of scattering), these oscillations can be understood as a gradient-field induced localization of the electron wave functions, 19 an effect known as Wannier-Stark localization. 20 Recently, it was shown that this localization might survive even when interactions are turned on, 21,22 yielding the concept of Stark many-body localization akin to many-body localization induced by quenched, quasi-periodic or programmable disorder. [23][24][25][26][27][28][29] These studies elevate closed electric-field driven quantum systems to the frontier of research concerning ergodicity breaking, and thus effects beyond the paradigm of statistical mechanics can be expected in these systems.
However, when considering interacting closed systems (as discussed above) under an external driving force and beyond the regime of many-body localization, the drive will continuously heat up the system until an infinite temperature state is reached by the growing deposition of energy. This state is not very interesting, but fortunately a more realistic model includes infinite baths 30 which can dissipate this additional energy. [31][32][33][34][35] In such a setup, an interesting non-equilibrium state (supporting, e.g., a finite value of the steady state current) is conceivable. In this context, negative differential conductance was reported if the electric field is increased at a constant coupling to the reservoirs. This negative differential conductance is a consequence of the current being suppressed for increasing fields, because the amount of energy per unit time dissipated by the bath remains constant, and thus the effective temperature increases, which in turn decreases the current. 31,32,36 Here we want to address the question of what happens to the electronic phases of matter as an increas-ingly strong field is driving the system out of its equilibrium state [37][38][39] from a microscopic model perspective. To this end, we extend a renormalization group approach, 40 which was successfully applied to characterize phases of matter in microscopic models in equilibrium, 41,42 to the non-equilibrium realm. This allows us to address how phases of matter can be controlled using non-equilibrium means via mechanism such as the dielectric breakdown of insulators. This mechanism is only one example out of the broader class of non-equilibrium control avenues and describes that a correlation driven Mott insulating state can be turned metallic after the field strength has surpassed a certain threshold value where the metallization occurs via the production of doublon-hole pairs. [43][44][45][46][47][48][49] While we apply the developed methodology to the case of an infinite one-dimensional nearest neighbor chain of spinless fermions, the general framework we derive allows to study non-equilibrium control of phases of matter in general tight-binding models, while keeping track of all of the microscopic details.
The rest of this paper is structured as follows: In section II, we introduce the class of models that can be treated using our methods. After briefly recapitulating the Keldysh Green's function formalism in Sec. III, we present a novel iterative algorithm to compute the Green's functions of an infinite system (Sec. IV). This algorithm is a central result of this paper and is widely applicable within all diagrammatic techniques such as (dynamical) mean-field theory. We develop a full secondorder implementation of the Keldysh functional renormalization group (that accounts for inelastic scattering) for infinite, open systems subject to an electric field in Sec. V. As an example, we then apply this methodology to an interacting tight-binding chain coupled to reservoirs in an electric field (Sec. VI). We thoroughly discuss numerical details, and we investigate the survival of the charge-density wave transition when the reservoir couplings and/or an electric field are switched on.

II. CLASS OF MODELS
First, we will outline the class of systems that can be treated using our method. We eventually aim at modelling infinitely extended, one-dimensional chains of charged fermions which are coupled to reservoirs and which are subject to an electric field. As a starting point, we consider a general fermionic Hamiltonian with an infinite number of degrees of freedom that features a kinetic energy as well as a two-particle interaction: where c ( †) denote the fermionic annihilation (creation) operator. We also refer to H sys as the chain. This system is assumed to be coupled to an infinite set of fermionic, non-interacting reservoirs: where a ( †) denote the fermionic annihilation (creation) operators within the reservoirs. The total Hamiltonian is given by: The initial state is assumed to be one where the reservoirs are decoupled (H ν coup = 0) and are by themselves in thermal equilibrium. The influence of the reservoirs can then be characterized by the following hybridization functions: where n ν (ω) is a thermal (Fermi) distribution function: We also assume that all degrees of freedom within the chain feature some decay channel into the reservoirs, guaranteeing a well-defined stationary state that is independent of the initial preparation of the chain itself. Guided by the picture of a chain in an electric field, we restrict ourselves to Hamiltonians which have a discrete translational shift symmetry. With a given L ∈ N defining a unit cell, we demand that The quantity E > 0 has the interpretation of an electric field in arbitrary units. The hybridization and distribution function are similarly required to fulfill which directly yields a similar relation for Γ ν,K ij . Finally, we assume that all terms in the Hamiltonian are strictly local, i.e., there exists an M ∈ N such that We stress that we do not impose any constraints on the size L of the unit cell.  Pictorial representation of the model defined in Eq. (9) and used in Sec. VI. A tight-binding chain with a nearest-neighbor hopping t and a nearest-neighbor interaction U is subject to an in-plane electrical field E and is coupled to wide-band reservoirs with a hybridization Γ. In addition, we consider a finite staggered on-site potential of strength s (not depicted).
In Sec. VI, we will discuss the specific example of an interacting tight-binding chain coupled to zerotemperature wide-band reservoirs (L = 1, M = 2; see All other components are uniquely defined by the symmetries of the system. This model will serve as the main testbed for our method. In the limit Γ = E = 0, the phase diagram can be computed analytically using the Bethe ansatz: 50 The system is a gapless Luttinger liquid for U ≤ 2t and a Mott insulator with a spontaneouslybroken translational symmetry for U > 2t, respectively. In the latter case, the ground-state is two-fold degenerate and features a charge-density wave (CDW). In order to break a potential ground-state degeneracy within a numerical method, we introduce a staggered potential s that increases (decreases) the on-site energies of odd (even) sites within the chain by s and therefore breaks translational symmetry. This increases the unit cell to L = 2.

III. GREEN'S FUNCTIONS
We will now introduce Keldysh Green's functions, which are a key ingredient to diagrammatic methods such as the FRG formalism. 51 The single-particle Green's functions in the stationary state can be written as The retarded component is given by and can be related to the non-interacting retarded Green's function g ret (ω) by virtue of the Dyson equation: where the self-energy Σ ret is associated with the twoparticle interaction v ijkl . The Keldysh Green's function is defined as and the corresponding Dyson equation reads where we have used that All quantities in Eqs. (11) and (14) are matrices defined by two single-particle indices. To simplify the notation, we will frequently employ multi-indices 1 = (i 1 , α 1 ) that include both this single-particle index i 1 as well as the Keldysh index α 1 ∈ {1, 2}. The frequency-dependence will still be written out explicitly. If the entire system is in an equilibrium configuration described by the (Fermi) distribution function n(ω), the Green's functions obey the fluctuation-dissipation theorem: The FRG approximation we introduce in Sec. V preserves this symmetry in the equilibrium limit, which is essential in order to avoid unphysical, anomalous heating effects. The symmetry described by Eqs. (6) and (7) translates directly to non-interacting Green's function g, where 1 + L denotes a shift of the single-particle index, 1 + L = (i 1 + L, α 1 ). This is a direct consequence of the Dyson equations (12) and (15). It follows from diagrammatic arguments (an expansion into an infinite perturbation series) that the exact self-energy Σ and thus also the full Green function G Λ (see, e.g., the Dyson equation) inherit this symmetry:

IV. COMPUTING GREEN'S FUNCTIONS IN AN INFINITE SYSTEM
In this section, we discuss how to compute the retarded and Keldysh Green's function of an infinite system under the assumption that the corresponding self-energies are known (e.g., from an FRG calculation). This cannot be done straightforwardly but requires an iterative algorithm, which we will now present. Our algorithm does not involve any additional approximations but is (numerically) exact. More importantly, it is not specifically tailored to the methods of this paper but is applicable in a completely general setting.
In the following, we will assume that i) our system fulfills the translation symmetry of Eqs. (6) and (7), that ii) the single-particle Hamiltonian h, the self-energies Σ ret,K , and the reservoir coupling Γ ν,ret ij are of limited range N where L evenly divides N , i.e., and that iii) the Green's functions G ret,K ij are only needed for |i − j| < N .
In our concrete example of the tight-binding chain [see Eq. (9)], we have h ij = Γ ν,ret ij = 0 for |i − j| ≥ 2. In Sec. V, we will show that our FRG approximation to the self-energy fulfills Eq. (19) and that only Green's function with |i − j| < 3M enter into the flow equations due to the approximations made in Eqs. (48) and (57). Thus, we would have N = 3M in this case.

A. Notation
For the rest of this section, we introduce the following notation for the single-particle indices: with the implicit understanding that, e.g., G ret CL (ω) refers to the retarded Green's function G ret (i∈C)(j∈L) (ω), which is a matrix of size N × N . The same convention is used for the self-energy as well as for all other quantities carrying two single-particle indices. We employ an Einstein convention for summations, e.g., We finally note that using our notation, the translation symmetry in Eq. (18) takes the form and likewise for the self-energy.

B. Retarded Green's function
We first discuss how one can obtain the retarded part of the Green's function. We need to invert a matrix which by construction has the following form: where a block structure is defined by T and D as follows: D CC , T CL , T LC , T CR , and T RC are matrices of size N ×N . Note that in general T LC = T † CL due to the inclusion of the self-energy. From now on, we will often omit the frequency dependence to improve readability.
A crucial ingredient is that the inverse of a block matrix is given by or equivalently By successively applying Eq. (25) and (26), one can prove that where in the first step we identified four blocks as indicated on the lhs. Per our assumption, the retarded Green's function G ret ij is only needed for |i − j| < N ; it is thus sufficient to determine G ret CC . If we apply Eq. (27) to Eq. (23) and use that T C(L\L) = 0, we obtain In order to solve Eq. (28), we need to determine the objects [D −1 ] LL as well as [D −1 ] RR , i.e., we need to calculate the first and last block of the inverse of the two matrices DLL and DRR associated with semi-infinite systems. This can be achieved iteratively by exploiting translationinvariance, which we will now discuss.

Iterative algorithm for the auxiliary Green's function
For notational simplicity, we define an auxiliary Green's function as the inverse of the matrix in Eq. (23) for T LC = T CL = T RC = T CR = 0 (which becomes block diagonal in this case). For later use, we introduce similar objects G ret, ¡ L and G ret, ¡ R as the inverse of Eq. (23) where only T LC = T CL = 0 and T RC = T CR = 0, respectively. The advanced components are defined as G adv, We note that while in the presence of finite interactions Σ ret = 0, these are no longer physical Green's functions of the underlying Hamiltonian, they inherit all of its symmetries such as translation-invariance: The first line is a direct consequence of Eqs. (6), (7), and (18) combined with the fact that the last N × N block (L) of an isolatedL system is, up to a shift in energy, identical to the last N × N block (C) of a system wherẽ R is removed (the second line follows similarly). If we use Eq. (30), we can now set up a recursion relation to where we have applied Eq. (25) to Eq. (23) with T CR = T RC = 0. A similar expression can be derived for (31) is local in ω and easily solved using a self-consistency loop. At finite electric field, however, this equation couples Green's functions at different frequencies. One can solve it by using lim ω→±∞ G ret, ¡ R (ω) = 0 as an initial condition; in practice, it is sufficient to set G ret, ¡ R (±Ω) = 0, where Ω far exceeds all other energy scales. Eq. (31) can then be used to successively calculate the auxiliary Green's function on a discrete grid of frequencies. How to do this in practice is outlined in Appendix A.

C. Keldysh Green's function
Next, we illustrate how to compute the Keldysh Green's function G K CC (ω). The Dyson equation (14) takes the form If we employ the lower-left component of Eq. (25), the Dyson equation can be simplified as follows: where we have defined as well as the auxiliary Green's functions The latter are the only unknown quantities in Eq. (34); have already been calculated in the previous section. In Eq. (34), we have employed that Σ K LR =Σ K C(L\L) =Σ K C(R\R) = 0 holds per our assumption. We will now discuss how G K, can be computed iteratively by exploiting translation-invariance.

Iterative algorithm for the auxiliary Green's function
In order to compute G K, ¡ L ¡ R LL and G K, ¡ L ¡ R RR , we will use the relation which is analogous to Eq. (33) and follows by applying Eq. (25) to Eq. (23) with T CR = T RC = 0. Moreover, we exploit translation-invariance forΣ K [see Eq. (18)] as well as for G ret,adv : which can be derived in analogy to Eq. (30). This yields In the last line, all quantities carry a frequency argument ω, which we have omitted to improve readability. G K, ¡ L ¡ R RR follows similarly. This equation has the same form as Eq. (31) and can be solved either self-consistently (if E = 0) or successively by utilizing G K, The functional renormalization group is an implementation of the RG idea on the level of correlation functions. 40 It sets up flow equations for the self-energy as well as for higher-order vertex functions with respect to a flow parameter Λ introduced as an infrared cutoff within the non-interacting Green's functions g ret,Λ (ω) and g K,Λ (ω). A detailed description to this method can be found in Refs. 40 and 52.
From now on, we focus solely on the case that H ν res describes zero-temperature wide-band reservoirs with a frequency independent hybridization that are coupled uniformly to the chain: In this case, it is convenient to use the single scale Γ as the flow parameter Λ = Γ. 34 The advantage of this approach is its physical interpretation -intermediate results during the solution of the flow equations can simply be viewed as physical results at a stronger coupling. We will now set up a full-fledged second-order Keldysh FRG implementation for a one-dimensional chain that is translation-invariant up to shifts in energy. A variety of different first-order Keldysh FRG calculations can be found in the literature, 53-56 but second-order FRG schemes have so far been developed exclusively for singleimpurity models 57,58 or for chains which are in thermal equilibrium. 41,42,59,60 The only key approximation in our scheme is the socalled channel decomposition of the vertex flow equation (see Sec. V B), which has been widely applied in equilibrium FRG calculations. 41,42,57,60 Moreover, we will assume that all vertex functions have a limited support of range M (see Sec. V E), which will serve as our key numerical control parameter. The original second-order flow equations (which only assume the channel decomposition) are recovered in the limit M → ∞, and we will demonstrate that convergence in M can be reached in all practical applications.

A. Flow equations
The flow equation for the self-energy reads The single-scale propagator is given by where ∂ * Λ indicates a derivative that acts only on the explicit Λ-dependence of the cutoff (but not on Σ Λ ). The quantity γ denotes the one-particle irreducible twoparticle vertex function; it preserves energy conservation due to the time-translation invariance of the system, and its frequency-dependence can thus be parametrized via with the coordinates Using this notation, the flow equation for γ takes the form where we already truncated the otherwise infinite hierarchy of differential equations by neglecting the flow of the three-particle vertex. This approximation is controlled in a perturbative sense and all terms neglected are at least of O U 3 . The flow equations (41) and (45) need to be complemented by an initial condition. When the coupling to the reservoirs is large (Λ → ∞), the vertex functions can be obtained analytically: where we introduced the Keldysh-space version of the two-particle interaction The initial value of the retarded self-energy is frequencyindependent and can therefore be absorbed into the noninteracting Hamiltonian h.

B. Channel decomposition
The vertex flow equation (45) depends on three independent frequencies and is thus difficult to tackle numerically. Hence, we need to resort to an additional approximation, the so-called channel decomposition. 57 We make the following ansatz for γ: γ Λ 1 2 12 (Π, X, ∆) =v 1 2 12 + γ p,Λ 1 2 12 (Π) + γ x,Λ 1 2 12 (X) + γ d,Λ 1 2 12 (∆), and assume that (i) the flow equation for γ p,Λ , γ x,Λ , and γ d,Λ is given by the first, second, and third term of Eq. (45), respectively, and that (ii) each channel is only fed back into its own flow equation. This yields The initial condition reads γ α,Λ→∞ 1 2 12 = 0. The self-energy flow equation (41) now takes the form The channel decomposition makes the vertex flow equations manageable by numerics (each term γ α,Λ depends only on a single frequency) but still include all terms of O U 2 . In addition to decoupling the frequency structure, the channel decomposition also simplifies the dependence on the spatial indices of the vertex functions. Since γ α,Λ→∞ = 0, one trivially finds that in this limit: where |1 − 2| := |i 1 − i 2 | refers to the distance of the single particle-indices within the multi-indices 1, 2 = (i 1,2 , α 1,2 ). Per the assumption in Eq. (8), the same holds true for the single-particle structure of the initial vertex v. One can easily see that the flow equations (49) preserve Eq. (51), which thus remains true throughout the flow; e.g., the indices 1 and 2 (2 and 1) appear as the first and last (second and third) argument on the rhs of the flow equation for γ x,Λ 1 2 12 . We emphasize that Eq. (51) is a direct consequence of the channel decomposition and does not constitute an additional approximation.
C. Making use of the system's symmetry The symmetries in Eq. (18) are self-consistently preserved within our approximation scheme (i.e., after truncation). If we assume that Eq. (18) holds for a given Σ Λ (and thus also for G Λ as well as S Λ ) and exploit that the initial, frequency-independent vertex fulfills Eq. (6), we can use the flow equation (49) to show that If we now plug Eq. (52) into the self-energy flow equation (50), it follows immediately that the symmetry relations in Eq. (18) are preserved. In a nutshell, Eqs. (18) and (52) imply that one of the spatial indices (say i 1 ) of the self-energy as well as of the two-particle vertex can be restricted to {0, . . . , L − 1} when solving the flow equations.

D. Integrations as convolutions
The flow equations (49) and (50) can all be rewritten in terms of convolutions: If we define the shorthand notatioñ and split up the self-energy flow equation (50) into three terms, Σ Λ = α∈{p,x,d} Σ α,Λ , we find No frequency-dependence is generated in the last term. Similarly, the flow equations (49) for the vertex can be recast as ∂ Λ γ p,Λ 1 2 12 (Π) = i 2π 33 44 γ p,Λ 1 2 34 (Π) G Λ 44 * S Λ 33 (Π)γ p,Λ 3 4 12 (Π) , ∂ Λ γ x,Λ 1 2 12 (X) = i 2π 33 44 This shows that a numerically-efficient implementation of the flow equations can be based on an efficient implementation of convolutions, which in turn can be achieved by employing fast Fourier transforms to perform all integrations (see Appendix B). While at T = 0 some components of the Green's functions and single-scale propagators are discontinuous, this is not true for the vertex functions, which one can understand as follows: The rhs of the flow equations (56) is governed by a convolution of two functions G Λ and S Λ that decay sufficiently quickly for ω → ±∞; this yields a continuous function.
Overall, we are left with a maximum of L(2M − 1) 3 independent, frequency-dependent components of the twoparticle vertex functions. In many cases (e.g., for a nearest neighbor interactionv), this number can be further reduced by imposing additional constraints on Eq. (51).

F. Single-scale propagators
In Sec. IV, we discussed how the Green's functions G Λ of an infinite system can be computed iteratively. Since the rhs of the FRG flow equations also contains the single-scale propagator S Λ , we will now illustrate how this quantity can be computed along the lines of Sec. IV. To improve readability, we will frequently refrain from writing out the Λ-dependence as well as frequencyarguments of various quantities such as G Λ (ω) and S Λ (ω) throughout this section.

Retarded single-scale propagator
We begin with the retarded part of the single-scale propagator; the advanced part follows from [S ret ] † = S adv . Since our cutoff stems from wide-band reservoirs coupled to each site, it only enters on the diagonal of [G ret ] −1 and thus This restriction reduces the number of terms in the following expressions, but a generalization to a more involved cutoff scheme is straightforward. The retarded part of the single-scale propagator can be computed by taking the derivative of Eq. (28): The only unknown quantity in this equation is S ret, ¡ L ¡ R := ∂ * Λ G ret, ¡ L ¡ R , which can be obtained via the derivative of Eq. (31): and similarly for S ret, ¡ L ¡ R RR . Yet again, this equation can either be solved self-consistently (at E = 0) or successively (for E = 0) using the boundary conditions S ret, ¡ L ¡ R (±∞) = 0 as outlined in Appendix A.

Keldysh single-scale propagator
We now proceed with the Keldysh component of the single-scale propagator. The cutoff only enters into the diagonal ofΣ K , which leads to [see Eqs. (4), (35), and (40)] Taking the derivative of Eq. (34) yields (all quantities carry a frequency argument ω)

G. Frequency discretization
For a numerical treatment, it is necessary to discretize the frequency space. In order to faithfully represent the physical system at hand, one must choose a grid that accounts for all of its relevant energy scales. For simplicity, we evaluate both the vertex functions as well as the Green's functions and single-scale propagators on the same set of frequencies.
The Green's functions decay on the scale of the system's bandwidth and are broadened by inelastic scattering. Furthermore, the Keldysh Green's function is linked to the distribution function within the reservoirs via Eq. (14), making the temperature of the reservoirs a relevant energy scale. This motivates the use of an equidistant grid whose width scales with the bandwidth and which includes additional points around the chemical potentials of the reservoirs. For the case of zerotemperature reservoirs, it is most convenient to choose with δ ω . At non-zero temperature, a more sophisticated choice of Ω extra is required. In Appendix B, we show that using such a grid allows for an efficient implementation of convolutions using fast Fourier transforms if the number of points in Ω extra is small.
When the coupling to the reservoir dominates all other energy scales, we require that ω max Λ δ ω . In the opposite limit, the bandwidth B = 4t of the closed system determines the width of the grid, ω max B, and δ −1 ω should be chosen much larger than the coherence time of the closed system. If not otherwise stated, we use Since ω max depends on Λ, it is necessary to adapt the grid during the solution of the flow equations. This is done after every step of the differential equation solver; in order to obtain the vertex functions on the new grid, we use linear interpolation (note that extrapolation is never required as ω max is only decreased).

VI. APPLICATION
We will now apply our iterative Green's function algorithm as well as our novel non-equilibrium FRG approach to the tight-binding chain introduced in Sec. II. A pictorial representation of this model is shown in Fig. 1.
First, we will benchmark the iterative Green's function algorithm introduced in Sec. IV against analytical results available for the non-interacting Hamiltonian U = 0. In this limit, all single-particle eigenfunctions of the closed system (Γ = 0) are exponentially localized for any nonzero E, leading to Wannier-Stark insulating behavior. The conductivity becomes finite for any Γ > 0.
Secondly, we will explore the capabilities of our novel non-equilibrium functional RG approach in describing finite interactions U > 0. There are two reference results that we can compare against. First, we will set up a mean-field treatment to compute the phase diagram in the large-U limit; this serves as a highly non-trivial test for the FRG, which is perturbative w.r.t. U . Secondly, it is known that the closed system with zero electric field (Γ = E = 0) undergoes a Berezinskii-Kosterlitz-Thouless quantum phase transition from a gapless Tomonaga-Luttinger liquid (U < 2t) to a charge density wave (U > 2t). The order parameter of the charge density wave (CDW) phase can be defined as the occupation difference between even and odd sites: The corresponding susceptibility diverges in a CDW phase but is finite otherwise. We will explicitly compute these quantities using the FRG for arbitrary U , Γ, and E.
A. Non-interacting properties

Density of states
We now discuss the physics for U = 0 in more detail and use this limit as a testing ground for our iterative Green's function algorithm. In the absence of an electric field, the local density of states (LDOS) is known analytically: In Fig. 3(a), we demonstrate how this exact result is recovered by the iterative algorithm of Sec. IV. For finite fields E > 0 and Γ = 0, the system is a Wannier-Stark insulator, and all single-particle eigenstates are exponentially localized. At small Γ t, the local density of states becomes sharply peaked on the scale of Γ. In Fig. 3(b), we demonstrate how our iterative algorithm can be used to compute the LDOS in a numerically-exact fashion. (One should note that for finite U , additional inelastic processes lead to smoother Green's functions, which simplifies computations.)

Current
Wannier-Stark localization has a profound impact on the current I flowing through the system in the presence of a finite electric field. In particular, the hybridizationdependence I(Γ) was determined analytically for small fields E t (note that the bandwidth is given by B = 4t): 34 The behavior in these three regimes can be understood qualitatively. For small Γ, Wannier-Stark localization leads to a vanishing current; the number of particles entering the chain from the reservoirs scales with Γ. Once in the system, each fermion is repeatedly reflected by the electric field in the form of Bloch oscillations. The typical distance traveled before eventually leaving the chain scales as t/E, resulting in a total current that scales as Γt/E. If E Γ t, fermions tunnel coherently between far-apart reservoirs. To analyze this in more detail, consider the defining equation for an element of the noninteracting Green's function: At ω = 0, E Γ t, and i j or j i, this equation is approximately solved by where I νν can be understood as the current between reservoirs ν and ν . For a given distance ∆ between pairs of reservoirs, there is O(∆) pairs that are connected by the selected bond. The function n ν (ω) − n ν (ω) has a width ∆E, and the Green's function can be assumed to be constant within such an interval. Therefore, the total current through that bond is obtained as At large coupling Γ t, the scattering into the reservoirs dominates and the correlation becomes small [compare Eq. (12)]: The current obtained using our numerical algorithm is shown in Fig. 4 as a function of Γ for various values of E < t. The three distinct regimes of Eq. (70) can be clearly identified; this serves as an additional benchmark for our iterative approach.

B. Large interaction limit
In the previous section, we have tested the (numerically exact) iterative algorithm introduced in Sec. IV by comparing with analytical results in the non-interacting limit U = 0. We now turn to benchmarking our FRG approximation scheme for U = 0. To the best of our knowledge, no reference data is available except in the limit E = Γ = 0. We therefore resort to a mean-field treatment, which is expected to give reasonable results for U → ∞. This is a highly non-trivial testing ground for the FRG, which is perturbative w.r.t. the interaction strength.

Mean-field approach
If the interaction dominates the bare, decoupled Hamiltonian (i.e., for U t), the system can be reduced to a chain of disconnected sites. One can easily show that for t = 0, the effects of a finite field E = 0 can be eliminated by virtue of a gauge transformation analogous to the Peierls substitution: 63 where τ denotes time (within the action). This transformation effectively shifts all energies on site j and the adjacent reservoir by −jE. Hence, we can restrict ourselves to E = 0. For Γ = 0, the ground state is a perfect charge density wave, while for large couplings Γ U , the ground state is not spontaneously ordered. To our best knowledge, the critical coupling strength which characterizes the transition between these two regimes is not known.
It is reasonable to treat the limit t = 0 using mean-field theory. Since the Green's functions become diagonal, one needs to solve the following self-consistency equation: Using Eq. (67) as well as n even = 1− n odd , this reduces to ∆n = 2 π arctan (s/Γ + U ∆n/Γ) .
At s = 0, this equation has nontrivial solutions beyond a critical interaction strength of In the disordered phase (U/Γ < π/2), the susceptibility for small symmetry breaking (s U, Γ) scales as For reasons of completeness, we will now present technical details of how to solve the mean-field equations for t = 0; results will be discussed elsewhere. In thermal equilibrium, a common technique to finding solutions of the mean-field equations is to use a self-consistency loop based on an initial guess for the field. If multiple solutions are found, one picks the one which minimizes the free energy; this solution corresponds to a stable fixed point of the mean-field equations.
For finite fields E = 0, the above procedure needs to be modified. First, a self-consistency loop is insufficient to identify all solutions of the mean-field equations as not all fixed points can be obtained as the limit of a self-consistency loop, regardless of the initial guess. One therefore needs to resort to a version of Newton's method, which also converges to fixed points where a self-consistency loop fails. Secondly, the free energy can no longer be used to determine a unique, stable solution. One can gauge the stability from the Lipschitz constant, which is estimated from the behavior of the self-consistency loop under a small perturbation. Out of equilibrium, one can thus only report the existence of solutions of the self-consistent equations and their stability with respect to perturbation. This will be discussed in a separate publication.

Benchmarking the FRG
The mean-field results in the limit U → ∞ can be used to benchmark our FRG approach (which is perturbative w.r.t. U ). Since h is diagonal, the same holds true for all single-particle propagators g (which can thus be computed straightforwardly without resorting to the iterative algorithm of Sec. IV) as well as for the self-energy. This allows us to greatly simplify the FRG flow equations. Note that the two-particle vertex functions remain infinitely extended as terms with arbitrary dist(1 , 2 ) and dist(1, 2) can be generated by the flow. However, γ Λ 1 212 = 0 if dist(1 , 1) > 0, and the self-energy itself is strictly local.
The main FRG results are presented in Fig. 5(a). For large reservoir couplings Γ, the CDW order parameter ∆n scales linearly in the initial symmetry breaking s, which corresponds to a finite susceptibility χ. Below a critical coupling, however, ∆n takes a finite value which does not decrease with decreasing s; the response of the system diverges, and translation invariance is spontaneously broken. One can identify the critical point Γ crit of this transition as the point where the curves at different s intersect. In Fig. 5(b), we show Γ crit as a function of M , which is the control parameter in the numerical solution of the FRG flow equations. For large M , we obtain which is in good agreement with the mean-field prediction of Eq. (79). For reasons of completeness, we show  Fig. 5(c), which again confirms that convergence can be reached. At the critical point, the susceptibility appears to diverge as a power law, and the FRG result for the critical exponent reads: This is illustrated in Fig. 5(d).

C. Phase diagram at intermediate interaction
At intermediate interaction, the ordering tendencies driven by the interaction compete with the kinetic energy, resulting in a non-trivial phase-diagram. In the absence of an electric field, mean field theory predicts an ordered phase at any interaction U when the coupling Γ is small enough (U crit /t → 0 for Γ → 0). In equilibrium and for Γ → 0 this is known to be an artifact of the mean-field approximation. The exact Bethe-ansatz solution predicts that the CDW order is destabilized by quantum fluctuations and that a finite critical U crit /t = 2 1 0.1 0.01 is needed to drive the phase into an ordered state by a Berezinskii-Kosterlitz-Thouless mechanism. 64,65 In contrast to the mean-field approach and in accord with the exact solution, we do not observe symmetry breaking at small interactions within the FRG calculation [see Fig. 6(a)]. While a phase transition at lower, inaccessible reservoir couplings Γ can in principle not be ruled out, an alternative, Matsubara FRG scheme carried out directly at Γ = 0 predicts a finite value of U crit /t. 41 Beyond a critical interaction of U crit /t ≈ 1.4, the initial symmetry breaking yields a response that does not vanish linearly for s → 0, and the system enters a CDW phase [see Fig. 6 After pinpointing the position of the phase transition, we can analyze the critical behavior itself. Mean-field theory predicts a divergence as independent of the interaction strength. In contrast, the susceptibility obtained using the FRG diverges as a power-law with an interaction-dependent exponent [see Fig. 7(a)]. The exponent increases upon lowering the interaction towards the critical point, beyond which no phase transition can be identified at any Γ > 0 [see Fig. 7(b)]. The critical interaction strength extracted from the analysis is approximately U crit /t ≈ 1.4 in the small-Γ limit (the exact solution yields U crit /t = 2).

D. Phase diagram at non-zero electric field
We now turn to the phase diagram of the system driven out of equilibrium by a finite electric field E > 0. Results for the susceptibility as well as for the CDW order parameter are shown in Fig. 8 and 9 for a constant electric field E/t = 0.2 and a constant interaction U/t = 5, respectively.
For a constant, small electric field of E/t = 0.2, the system is in a disordered phase for U/t 3 at any value of Γ. However, the susceptibility does not decrease monotonically with Γ but features novel structures which reflect the emergence of the multiple energy scales away from equilibrium [see Fig. 8(a)]. The critical interaction beyond which we observe a transition into a CDW phase is drastically enhanced (here by roughly a factor of 2) compared to the case of E = 0; the electric field drives a current through the system, and the tendencies to form charge order are suppressed.
In Fig. 9, we show results for constant U/t = 5 and various E. Small fields E/t 2 induce a current and thus reduce the tendency to form charge order. However, we observe a re-entrance into the CDW phase as the strength of the field is increased: For large E, the current is suppressed due to localization and can no longer inhibit CDW order.

VII. CONCLUSION
We have introduced a new method to treat an infinite quantum system which is driven out of thermal equilibrium in a translation-invariant fashion. In particular, we have set up recursion relations for Keldysh-Schwinger Green's functions and demonstrated how they can be solved efficiently in a numerically-exact way. The algorithm is directly applicable to any diagrammatic Green's function based method such as (dynamical) mean field theories or the functional renormalization group. Furthermore, the presented iterative Green's functions scheme has direct relevance to periodically driven systems described by a Floquet Green's functions approach. 20,66,67 We applied this general machinery within a novel, second-order Keldysh formulation of the functional renormalization group; significant effort was devoted to efficiently implementing the FRG flow equation via fast Fourier transforms. As a physically relevant example, we studied a tight-binding chain of interacting spinless fermions coupled to reservoirs and driven out of equilibrium by an electric field. It is known that the closed system features an equilibrium Berezinskii-Kosterlitz-Thouless phase transition into a CDW phase beyond a critical interaction strength U crit . In contrast to meanfield approaches, the FRG reproduces this result and can be used to determine U crit in the presence of a finite reservoir coupling. A small electric field induces a current and thus suppresses tendencies to form charge order. For large field, however, Wannier-Stark localization leads inhibits currents and leads to a re-entrance into the CDW phase.
where G ∈ G ret, . At vanishing E = 0, such equations can be solved with a self-consistency loop. For E > 0, Eq. (A1) is non-local in the frequency, and we will now discuss to to solve such an equation efficiently. In all cases, the following boundary condition holds: With the initial assumption G(±k max N E) = 0, Eq. (A1) can be used to successively calculate the Green's function on the entire grid.
b. Arbitrary grids The method described above becomes inefficient if N E is much smaller than the required grid spacing. Alternatively, we can consider an arbitrary frequency discretization {ω 1 , . . . , ω kmax } and assume G(ω kmax ) = G(ω 1 ) = 0. For simplicity, we restrict us to the case of negative sign in Eq. (A1). In that case one uses the following recursive algorithm to successively obtain G(ω) on the entire grid. One assumes that G(ω) has already been calculated for all frequencies ω k+1 , . . . , ω kmax . In order to compute ω k , we find the largest frequency ω i which fulfills Note that one always has i ≥ k. If we assume that G is well approximated by a linear interpolation between grid points, we obtain: where Whenever k = i, Eq. (A5) is solved using a selfconsistency loop, otherwise G(ω i ) has previously been

Appendix B: Efficient convolution
For a numerical implementation of the FRG algorithm discussed in this work, it is essential to perform the integrals appearing on the rhs of the flow equations efficiently. As mentioned above, this can be achieved by treating all integrals as convolutions and by utilizing an efficient algorithm to perform these. To that end, we discretize all frequencies on a grid as discussed in Sec. V G. In the following, N denotes the total number of frequency points.
We define the convolution of two functions as If this integral is carried out naively, this requires O (N ) operations, and the effort to obtain the rhs of all of the flow equations thus scales as O N 2 . In this Appendix, we will discuss how such convolutions can be obtained efficiently on equidistant grids, on grids with an arbitrary spacing, and on mixed grids (which are required within the FRG).

Equidistant grids
For equidistant grids, we can rewrite the convolution in terms of a discrete Fourier transform, which will allow us to perform significantly faster computations.

a. Discretization
We employ an equidistant grid with N ∈ 2N + 1 points and approximate f and g by piecewise constant functions (compare left panel of Fig. 10): with a vector (h i ) i={− N −1 2 ... N −1 2 } , and h i = h (iw) and h = f, g. A convolution of f, g then reduces to a convolution of vectors: with the summation bounded appropriately. The summation runs over O (N ) elements, and it thus takes O N 2 operations to compute the convolution on the entire original grid (which has N points). This can be improved by employing a Fourier transform.

b. Fourier transform
Using the discrete Fourier transform one rewrites the discrete convolution as Nf k1 e −i(n−j)k2 2π which is the n-th element of the back-transform of the product off kḡk . By employing the fast Fourier transform (FFT) algorithm, the discrete Fourier transform can be computed in O (N log(N )) operations, 68 which allows us to obtain all components of the convolution on the same grid in O (N log(N )) operations.
While this algorithm is very efficient, it is specifically designed for equidistant grids and can not easily be generalized to more general discretizations. For our application with its vastly different energy scales, one is forced to use an overly dense grid, which diminishes the advantage of this method.

Arbitrary grids
Alternatively, one can work with an entirely arbitrary grid defined via {x 1 , . . . , x N } with x 1 < · · · < x N . We assume that the functions f, g are approximated piecewise linearly (compare center panel of Fig. 10), for h = f, g. The convolution of the two functions can then be written as Depending on the chosen discretization, the support S kk of the integrand at a given k might be non-zero for more than one k ; however, the remaining convolutions of linear functions can be easily evaluated.
Note that for an efficient algorithm it is of crucial importance that the grid of x i is sorted and that appropriate algorithms to identify the contributing k are employed. On an arbitrary grid with N points, this algorithm requires O N 2 log N operations to obtain the convolutions on the entire grid.
While this algorithm allows for grids that take the vastly different energy scales of the problem into account, it is (compared with the case of equidistant grids) slow and is the bottleneck of such an implementation. with − N −1 2 w < x 1 < x 2 < · · · < x k < N −1 2 w. We first analyze sections of the equidistant subgrid.
Then there is exactly one decomposition f (x) ≈ a n + f n lin (x) fulfilling the following three conditions: Firstly, we require f (x) = a n + f n lin (x) ∀x ∈ {nw} ∪ S n , where f lin is a piecewise linear function with support only in nw − w 2 , nw + w 2 . Secondly, we require df n lin (x) dx → 0 for x → nw ± w 2 . (B14) Lastly, we impose the condition with f 0 (x) = a n |x − nw| < w 2 0 otherwise.
Note that only a small number of f n lin (x) ≡ 0 (in fact, k at most), since these linear functions are only required if one of the additional points x i falls into the interval nw − w 2 , nw + w 2 . A visual example of such a decomposition is shown in Fig. 11, a concrete example is presented in the next section.
We now evaluate this decomposition on the grid and have to distinguish between two cases: a. For all y in the equidistant grid the first term can be computed as outlined in Sec. B 1 in O (N log(N )) operations. Due to our choice dxf i lin (x) = 0 = dxg i lin (x), the support of these functions as well as the fact that f 0 and g 0 are piecewise constant, the second and third term vanish. The last term can be explicitly computed in O k 2 log(N ) operations: the summands are only non-zero for O k 2 values of (i, j), and for each such combination O (log(N )) operations are required to identify the (small) support of f i lin * g j lin . b. For each y outside the equidistant grid one can obtain the convolution explicitly using the algorithm for arbitrary grids discussed in Sec. B 2; as there are only k such points, this results in an algorithm of O (kN log(N )) operations.