# Nonequilibrium self-energy functional approach to the dynamical Mott transition

###### Abstract

The real-time dynamics of the Fermi-Hubbard model, driven out of equilibrium by quenching or ramping the interaction parameter, is studied within the framework of the nonequilibrium self-energy functional theory. A dynamical impurity approximation with a single auxiliary bath site is considered as a reference system and the time-dependent hybridization is optimized as prescribed by the variational principle. The dynamical two-site approximation turns out to be useful to study the real-time dynamics on short and intermediate time scales. Depending on the strength of the interaction in the final state, two qualitatively different response regimes are observed. For both weak and strong couplings, qualitative agreement with previous results of nonequilibrium dynamical mean-field theory is found. The two regimes are sharply separated by a critical point at which the low-energy bath degree of freedom decouples in the course of time. We trace the dependence of the critical interaction of the dynamical Mott transition on the duration of the interaction ramp from sudden quenches to adiabatic dynamics, and therewith link the dynamical to the equilibrium Mott transition.

###### pacs:

71.10.-w, 71.10.Fd, 71.15.Qe, 05.70.Ln## I Introduction

Systems of strongly correlated electrons on a lattice exhibit diverse emergent phenomena such as insulating behavior caused by strong local Coulomb repulsion. The Mott insulator is believed to hold the key for an understanding of the complex physics that is characteristic for several transition metals and compounds and has been at the focus of a vast number of studies in recent decades.(Imada et al., 1998) The preparation and the study of a Mott-insulating state by experiments done with ultracold fermionic atoms (Jördens et al., 2008; Schneider et al., 2008) opens an exciting new perspective on the underlying many-body problem and on the related idealized many-body models, such as the single-band Fermi-Hubbard model.(Giamarchi, 2004) The new aspect in those experiments is the high degree of control over the microscopic model parameters which can be used to steer the system between different phases in the equilibrium phase diagram but also to initiate and to manipulate nonequilibrium processes. (Bloch et al., 2008; Lewenstein et al., 2012) The last point in particular has attracted some interest recently, and time-dependent experiments with ultracold atoms in optical lattices (Greiner et al., 2002; Strohmaier et al., 2010) can complement the investigation of condensed-matter dynamics on femtosecond time scales with ultrafast pump-probe experiments. (Iwai et al., 2003; Perfetti et al., 2006; Wall et al., 2011)

The study and the understanding of strongly interacting lattice-fermion systems far from equilibrium requires a critical examination of standard concepts of quantum statistics regarding, e.g., the thermalization of isolated quantum systems,(Deutsch, 1991; Srednicki, 1994; Rigol et al., 2008) and it can bring about entirely new concepts such as dynamical phase transitions. (Dziarmaga, 2010; Polkovnikov et al., 2011) Apart from such fundamental theoretical questions, a further development and application of numerical methods is highly needed to study relevant problems such as the correlation-driven metal-insulator transition.

The Mott transition in the single-band Hubbard model is the paradigmatic field of application for the dynamical mean-field theory (DMFT),(Metzner and Vollhardt, 1989; Georges et al., 1996) which may be characterized as an internally consistent and nonperturbative mean-field approach controlled by the limit of infinite spatial dimensions. Although the feedback of nonlocal magnetic correlations is neglected, the DMFT phase diagram of the (paramagnetic) Mott transition (Georges et al., 1996) represents an instructive example of a phase transition which must be described by nonperturbative means. This has triggered the study of the real-time dynamics at the Mott transition using the extension of the standard DMFT to the nonequilibrium case.(Schmidt and Monien, 2002; Freericks et al., 2006; Aoki et al., 2014) The simplest protocol to initiate the dynamics is a sudden quench of the Hubbard-, from to different final values . For weak , it has been found that thermalization is delayed and that the system gets trapped in a so-called prethermal metastable state.(Moeckel and Kehrein, 2008; Eckstein et al., 2009) On the other hand, for strong , a fast relaxation is again impeded by collapse-and-revival oscillations that are reminiscent of the dynamics in the atomic limit. Both regimes are separated by a sharp transition for a certain final interaction strength , at which a fast relaxation to thermal equilibrium takes place. This picture of the “dynamical Mott transition” has emerged by applying the nonequilibrium DMFT to the Hubbard model on an infinite-dimensional Bethe lattice, (Eckstein et al., 2009) and has been corroborated by different subsequent studies and using different methods.(Schiró and Fabrizio, 2010, 2011; Sandri et al., 2012; Hamerla and Uhrig, 2013, 2014)

Within the DMFT the original lattice-fermion system is mapped onto an impurity model where a single correlated site is embedded in a noninteracting dynamical mean field (“bath”) which must be determined self-consistently. This effective impurity model poses a demanding many-body problem, particularly for a general nonequilibrium situation. While continuous-time quantum Monte-Carlo methods represent highly efficient “solvers” for the equilibrium case, only short propagation times can be accessed due to a severe sign (or phase) problem showing up in the real-time domain. (Werner et al., 2009; Gull et al., 2011) Perturbative approaches are much more successful and have been used extensively in the strong- (Eckstein and Werner, 2010) and weak-coupling regime (Eckstein and Werner, 2011; Tsuji et al., 2013) but are clearly of limited use to address the Mott transition which takes place at intermediate coupling strengths. As concerns the single-band Hubbard model, an exact-diagonalization (ED) solver is an efficient method for the equilibrium case at zero temperature,(Caffarel and Krauth, 1994) though the determination of the one-particle parameters of the impurity Hamiltonian is essentially done in an ad hoc way. For a nonequilibrium problem, the necessary Hamiltonian representation of the effective mean field poses an even more severe complication which has so far been solved on a short time scale only (Gramsch et al., 2013) since in the course of time more and more bath degrees of freedom have to be coupled to the impurity. Still this has allowed ED approaches to operate, such as Krylov-space methods, (Gramsch et al., 2013) the multiconfiguration time-dependent Hartree method (Balzer et al., 2015) or density-matrix renormalization techniques based on matrix-product states.(Wolf et al., 2014)

The (nonequilibrium) self-energy functional theory (SFT),(Potthoff, 2003a; Hofmann et al., 2013) offers a different route for the application of ED methods. Here, a reference impurity model with a given finite (small) number of bath degrees of freedom is considered. Instead of imposing the DMFT self-consistency condition, the time-dependent parameters of the reference system are fixed by applying a general variational principle, stating that the grand potential should be stationary as when expressed as a functional of the (nonequilibrium) self-energy. While the DMFT is recovered by choosing a reference system with a continuum of bath sites, any reference system with a finite bath generates a consistent dynamical impurity approximation (DIA) or, when choosing, in the spirit of cluster-DMFT approaches, a finite cluster of correlated sites to better account for short-range correlations, a variational cluster approximation (VCA). For an overview we refer to Refs. Potthoff, 2012, 2014.

This idea has successfully been employed for the study of the (equilibrium) Mott transition. Remarkably, the DIA with only a single additional bath site has turned out to recover the DMFT picture of the Mott transition in a qualitatively correct way, (Potthoff, 2003b; Požgajčić, 2004; Koga et al., 2005; Inaba et al., 2005; Eckstein et al., 2007) and with a few more bath degrees of freedom (Požgajčić, 2004) the agreement is even quantitative.

The present paper reports on results obtained by applying the nonequilibrium extension of the dynamical impurity approximation to study the real-time dynamics after quenches and ramps of the Hubbard interaction. For the nonequilibrium case, the fundamental concepts of self-energy functional theory are somewhat different and require a completely new strategy for the numerical evaluation.(Hofmann et al., 2013) Different practical issues of the implementation have been discussed recently in Ref. Hofmann et al., 2016. For the first application of the nonequilibrium DIA we restrict ourselves to a two-site reference system. As will be discussed, this is indeed sufficient in many respects to cover the essentials of the dynamical Mott transition as compared to previous work.(Eckstein et al., 2009; Schiró and Fabrizio, 2010) We also discuss the significance of the method in view of recently proposed Hamiltonian-based impurity solvers.(Gramsch et al., 2013; Balzer et al., 2015; Wolf et al., 2014)

The paper is organized as follows: In Sec. II we briefly introduce the nonequilibrium many-body problem, review the self-energy functional approach, and discuss the cornerstones of its numerical implementation. The application to study the Mott transition for both the equilibrium and the nonequilibrium case is presented in Secs. III and IV, respectively. Section V provides a summary.

## Ii Model and Methods

Using standard notations, the Hamiltonian of the Fermi-Hubbard model at half-filling reads

(1) |

Here, (creates) annihilates a fermion at site and with spin projection , and the number operator is given by . Fermions can tunnel between neighboring sites with the hopping amplitude . Two fermions on the same site are repelled by the local Coulomb interaction . Nontrivial real-time dynamics can be stimulated by controlling the explicit time dependencies of the model parameters. Here, we investigate both sudden quenches and ramps of different duration of the interaction parameter .

Calculations have been performed within the two-site dynamical impurity approximation (DIA), which provides a local approximation to the self-energy. It relates the full lattice problem to a small reference system consisting of a single correlated site () sharing the same time-dependent interaction as the original model but with an additional uncorrelated “bath” site () coupled to it via the time-dependent hybridization . For an illustration, see Fig. 1. Here and in the following, primed quantities refer to the reference system.

The parameter is determined via a variational principle set up within the general framework of the nonequilibrium self-energy functional theory (SFT).(Hofmann et al., 2013) The latter exploits the fact that the initial state equilibrium grand potential of the original system, , at inverse temperature can be expressed as a functional of the nonequilibrium self-energy. The self-energy functional is stationary at the physical self-energy of the model, i.e., , where it equals the physical grand potential, namely . Note that functionals are indicated by a hat and that is defined on the Keldysh-Matsubara contour ,(Wagner, 1991; van Leeuwen et al., 2006; Rammer, 2007) i.e., has elements with complex contour times .

Assuming that the problem posed by the reference system can be solved, the self-energy functional of the original system can be evaluated exactly for a certain subclass of trial self-energies, namely for the exact self-energies of the reference system at different parameters . We have

(2) |

where the grand potential , the self-energy , and the Green’s function of the reference system are functionals of . Furthermore, is the inverse temperature of the initial equilibrium state. One must consider with on the upper and lower Keldysh branch as independent variables, since otherwise the functional dependence of on all real-time quantities would disappear. Finally, is the free Green’s function of the original system. The trace contains an implicit integration along , i.e., we defined , where is infinitesimally later on . The optimal value is determined at each instant of time according to the stationarity principle

(3) |

which is evaluated on the space of physical parameters . This Euler equation is the central equation of the (two-site) dynamical impurity approximation, which gives us access to the optimal local self-energy and therewith to the one-particle Green’s function .

Equation (3) is inherently causal, i.e., optimal parameters can be determined at some time without affecting results at earlier times . This allows for the implementation of a time-propagation scheme for . For reasons discussed in Ref. Hofmann et al., 2013, it is beneficial to carry out the functional derivative in Eq. (3) analytically, which turns the variational problem into the problem of finding the zeros of . The latter, however, proves to be highly unstable. (Hofmann et al., 2016) Fortunately, this problem can be bypassed by only fixing the initial condition (at ) via and for later times requiring the respective time derivative to vanish, i.e., . Thus, we finally have to solve

(4a) | ||||

(4b) |

For precise details on the SFT framework and its numerical implementation we refer to Refs. Hofmann et al., 2013, 2016.

In principle, approximations within the self-energy functional theory can be constructed such that they respect the macroscopic conservation laws and the respective continuity equations for particle number, spin, and energy. Conservation laws for one-particle quantities in fact hold for any choice of the reference system, but obeying energy conservation requires a continuum of variational degrees of freedom.(Hofmann et al., 2013) Hence, for small reference systems, energy conservation is weakly violated, but one can expect to gradually improve on this by adding further variational degrees of freedom. It is noteworthy, that by providing a continuous bath (i.e., ) one formally recovers the dynamical mean field theory (DMFT),(Georges et al., 1996; Freericks et al., 2006; Schmidt and Monien, 2002; Aoki et al., 2014) (for or its cluster extension for ) which in fact is a fully conserving approximation.

For our calculations we consider the half-filled Hubbard model on a one-dimensional lattice of 40 sites with periodic boundary conditions, which is sufficient to ensure numerically converged results.(convergence, ) Choosing a one-dimensional system is convenient for numerical reasons. It is important to note, however, that the lattice dimension and geometry enters the DIA only via the free density of states (DOS). Moreover, we expect that results essentially depend on the variance of the DOS only.(Potthoff, 2003b) For the one-dimensional lattice this is . Energy (and time) units are fixed by setting . Calculations have been performed for different inverse temperatures , which set the length of the Matsubara branch in Eqs. (2)–(4). All integrations over imaginary time have been carried out using accurate high-order numerical integration schemes with step sizes varying from for larger to for smaller . For the real-time propagation and integration along both Keldysh branches, however, we are limited to the trapezoidal rule.(convergence, ) Sufficiently converged time propagations are obtained for time steps for maximum times up to .

## Iii Equilibrium Mott transition

Before the real-time dynamics of the Hubbard model can be analyzed within the two-site DIA, a proper initial state has to prepared, that is, the equilibrium variational problem [Eq. (3) at ] has to be solved. To the best of our knowledge, all previous SFT studies evaluated the grand potential [cf. Eq. 2] (and possibly its derivatives(Požgajčić, 2004)) directly to search for stationary points. Here we instead determine the equilibrium solution by evaluating its derivative analytically and look for the roots of by solving Eq. (4a). Thus, as a benchmark and to introduce our method, in this section we reproduce and discuss the known equilibrium results for a two-site reference system, as depicted in Fig. 1.

Results for the optimal hybridization parameter are shown in Fig. 2(a). In general, the on-site energies of the correlated and of the bath site have to be used as additional variational parameters, but in the present case (half filling) their value is fixed due to particle-hole symmetry. Starting from high temperature (low ) and weak interaction , one can easily perform a global search to obtain a solution of Eq. (4a). The full – phase diagram is then explored using a local search based on Broyden’s method,(Kelley, 1987; Broyden, 1965) starting from the solution at a nearby point in the phase diagram. By lowering the temperature , i.e., increasing , we find three solutions for certain values of : There is a metallic solution with large , which is adiabatically connected to the metallic phase at weak , an insulating one with small , connected to the strong- limit, and a third solution with intermediate , which is thermodynamically unstable (see below). This indicates a first-order phase transition with coexisting metallic and insulating phases.

In Fig. 2(b) we additionally present the double occupancies for the respective optimal solutions. In the coexistence region the double occupancy, like the optimal hybridization, has three branches. The branch for which increases with increasing corresponds to the thermodynamical unstable solution, because for this phase one has , which violates a thermodynamical stability condition. In fact, the system undergoes a first-order phase transition at a critical interaction, the value of which can be inferred from the Maxwell construction,(Strand et al., 2011) as shown in the inset of Fig. 2(b): The double occupancy jumps at the critical interaction , determined by requiring that the shaded areas on both sides of the jump be equal. In addition, the lower and upper boundaries of the coexistence region, and , can be read off at the spinodal points of the curve [see arrows in Fig. 2(b)].

Results for different temperatures are collected in the phase diagram shown in Fig. 3. Metallic and insulating solutions coexist in a triangular-shaped region, bounded by the curves and . Within the coexistence region, there is a line of first-order transitions terminating in a second-order critical point at the temperature . For temperatures above the Mott metal-insulator transition becomes a smooth crossover. Extrapolating our data to , we find and , both of which fall within a range of results obtained earlier for other lattices (where the variance of the density of states has been used as the energy unit).(Potthoff, 2003b) The value of obtained within the DIA for the Bethe lattice is in remarkably good agreement with DMFT+NRG(Bulla, 1999) (). The value obtained for the critical temperature, , is more sensitive to the lattice geometry: For the semi-elliptical DOS one finds within the two-site DIA.(Potthoff, 2003b) This value underestimates the DMFT result by 50%, but quantitative agreement is obtained already by adding only three bath sites.(Požgajčić, 2004)

As there is no Mott transition at a finite Hubbard- in the one-dimensional model, (Giamarchi, 2004) let us point out again that the DIA is a mean-field approach. It is therefore not really sensitive to the lattice dimension, and the results rather depend on the lattice geometry via the variance of the DOS only. The one-dimensional case is studied here for purely technical reasons, and the results should be seen as representative for the model on higher-dimensional lattices.

We conclude that the present implementation of the two-site DIA reproduces the known results for the Mott transition obtained earlier where the phase diagram has been constructed from the explicit calculation of the grand potential. The agreement between the results of the two different numerical approaches is fully quantitative. As compared to the full DMFT solution, the two-site approximation qualitatively captures the correct topology of the equilibrium phase diagram. Quantitatively, is predicted quite accurately while is over- and is underestimated. For the present study, we will nevertheless restrict ourselves to the two-site approximation since the computational effort is considerably higher for the nonequilibrium case. More importantly, nonequilibrium calculations with more than a single bath site are not easily stabilized numerically with the present implementation. (Hofmann et al., 2016)

## Iv Dynamical Mott transition

In the following we will discuss the real-time dynamics of the Hubbard model induced by sudden quenches or ramps of the interaction from a “free” initial state to arbitrary final interactions . Within the SFT the optimal parameters of the reference system are undefined for a free system due to the vanishing self-energy. For practical reasons we therefore consider initial states with .

We furthermore choose the initial inverse temperature . Essentially, this corresponds to a zero-temperature initial state: As concerns the reference system, there is hardly any change in the optimal hybridization parameter with temperature in the limit , as can be seen from Fig. 2(a). The remaining temperature dependence via the noninteracting Green’s function of the original system [see Eq. (2)] is very weak for lower temperatures. Consequently, there is hardly any temperature dependence seen in the nonequilibrium results. This has been verified numerically (up to ).

The interaction is switched from to via by either quenching,

(5) |

or conducting cosine-shaped ramps of different duration , i.e.,

(6) |

Both cases will be discussed successively in the next two subsections.

### iv.1 Interaction quenches

Following the time evolution after a quench, we find two qualitatively different response patterns for weak and strong final interactions, which are well separated by a sharp transition point at a critical interaction . Results are presented in Fig. 4, where we show the time dependence of the optimal hybridization parameter, the double occupancy, and the total energy. Moreover, in Fig. 5 we show for each time-dependent quantity the time average

(7) |

and the fluctuations

(8) |

Let us first focus on weak quenches, i.e., . For the optimal hybridization parameter we observe a quick drop to smaller values within approximately one inverse hopping, followed by moderate oscillations around some constant value, see Fig. 4 (top left). For final interactions , the long-time average of the optimal hybridization slightly decreases with increasing (Fig. 5, top). On the same short time scale, the double occupancy decays from its noninteracting initial value, i.e., the Coulomb repulsion quickly suppresses doubly occupied sites. For final interactions we find a strong initial drop and pronounced periodic recurrences. However, these recurrences shift to later and later times upon increasing . In addition, small regular oscillations around some value close to zero become apparent, see Fig. 4 (middle left).

The exact value of the total energy right after the quench (at ) is given by the expectation value of the Hamiltonian in the noninteracting initial state, , which increases linearly with the final interaction. For weak quenches we find that this value is relatively well conserved, apart from a small drop of the time-averaged value (of less than ), and some moderate oscillations of about or less for increased quench size (up to ) and when compared to the respective long-time average (Fig. 5, bottom). By comparison with a thermal ensemble for the same interaction , i.e., by comparing with equilibrium two-site DIA calculations, we can thus ascribe an effective temperature to the long-time averages by demanding that . The effective temperature increases from for to for . The corresponding thermal value for the double occupancy roughly agrees with the respective time-averaged value (see Fig. 5, middle). This is in agreement with the prethermalization scenario (Moeckel and Kehrein, 2008; Eckstein et al., 2009) observed in DMFT calculations.

We now turn to strong quenches, i.e., . Here the time-dependent behavior of the system drastically differs from that in the regime of weak quenches. For we find a quite regular oscillatory behavior for all relevant quantities. A Fourier analysis of the oscillations in the double occupancy reveals that oscillations occur with frequencies approximately given by , as shown in Fig. 6, i.e., the characteristic frequency for collapse-and-revival oscillations in the atomic limit. On top of this, there are slow beatings, which probably should be ascribed to finite-size effects and which appear to be independent of the interaction. In the long-time limit, the double occupancy slowly increases with . However, it does not reach its free value (i.e., ) again, as perturbative arguments would suggest,(Eckstein et al., 2009) i.e., the two-site approximation seems to underestimate the actual double occupancy in the strong-coupling limit (see Fig. 5, middle).

The optimal hybridization parameter strongly oscillates around zero, see Fig. 4 (top right). In equilibrium and for strong interactions the quasiparticle weight (for a two-site system) is given by ,(Lange, 1998) so that strong collapses and revivals of the square of the optimal parameter would correspond to an oscillatory behavior of the Fermi-surface discontinuity, as has been observed in DMFT calculations. (Eckstein et al., 2009) In Fig. 5 (top) we therefore show the long-time behavior of . With increasing interactions we find that both its average and its fluctuations quickly saturate.

For strong final interactions, conservation of total energy becomes rather poor and, compared to the weak-coupling case, the time-averaged value differs more significantly from the exact value. Nevertheless, we may again compare the long-time average to an appropriate thermal value to extract the effective temperature . As compared to the weak-coupling regime, the effective temperatures are roughly an order of magnitude higher and, apart from an offset, scale linearly with . Interestingly, from this we find thermal values of the double occupancy which in fact coincide with the overall minima of its time-dependent oscillations, see Fig. 5 (middle). This is again in line with the DMFT results.(Eckstein et al., 2009)

We now focus on the dynamics close to the critical point , see Fig. 7. In this regime the behavior of for and becomes very similar: Within two inverse hoppings, the optimal hybridization strength decays almost to zero, but then revives to positive values for and shows slow oscillations with relatively large amplitude. The same dynamics, but with opposite sign of the optimal parameter at long times, is observed for . This is accompanied by a decay of the double occupancy down to almost zero (), followed by strong revivals which are in phase with . As discussed for the weak-quench regime, these oscillations shift to later and later times for quenches closer and closer to the critical value. Finally, right at the critical point, no revivals are observed, i.e., the bath dynamically decouples, and remains zero up to the longest simulated times. For the double occupancy merely shows weak oscillates around its long-time average.

We conclude that, within the two-site DIA, the dynamical Mott transition is described as a sharp transition characterized by critical behavior in the -dependence of the quantities shown in Fig. 5. One may speculate that in calculations with more bath degrees of freedom in the DIA, some bath sites which represent low energy degrees of freedom would decouple whereas others would remain connected to the correlated impurity. Nevertheless, even on the level of the two-site approximation, there is a surprisingly good agreement of the critical interaction with results from the DMFT(Eckstein et al., 2009) () and the Gutzwiller ansatz(Schiró and Fabrizio, 2011) () when comparing with the value rescaled by the variance of the one-dimensional DOS, i.e., with .

Within the DMFT (Eckstein et al., 2009) a rapid thermalization is found at , and the thermalized state is characterized as a bad metal. Opposed to this, within the two-site DIA, a complete decoupling of the bath site right at the critical point implies that the final state is described on the level of the Hubbard-I approximation. (Hubbard, 1963) Note that therewith the dynamical Mott transition is very similar to the equilibrium Mott transition at zero temperature which is also characterized by a vanishing hybridization to the bath site. In both cases, the Hubbard-I approximation must be seen as a comparatively crude description of the bad metal or Mott insulator, respectively, and one cannot expect a fully consistent picture on this level. In the nonequilibrium case, for example, the determination of the effective temperature by comparison with equilibrium two-site DIA calculations via yields . The resulting thermal double occupancy of , however, turns out too large as compared with the time average . A better agreement is found when estimating by comparing with the Hubbard-I solution, where the bath site is decoupled. This yields and . However, at finite temperatures, the Hubbard-I solution is only metastable. One may summarize that for the final state at , our findings more resemble the predictions of the Gutzwiller approach (Schiró and Fabrizio, 2011) rather than those of the DMFT.

It is also instructive to interpret our results in the context of a recent proposal (Schiró and Fabrizio, 2011; Žitko and Fabrizio, 2015) stating that the metallic phase of the paramagnetic and particle-hole symmetric Hubbard model in infinite dimensions can be seen as a phase with a spontaneously broken local gauge symmetry. At the zero-temperature transition to the Mott insulating phase, the symmetry is restored. This picture of the Mott transition is nicely captured by the mean-field theory of a slave-spin model, (de’Medici et al., 2005; Rüegg et al., 2010) which is essentially equivalent to the Gutzwiller approach. Using the numerical renormalization group, it has been demonstrated (Žitko and Fabrizio, 2015) that the symmetry breaking at the Mott critical point takes place in the full model and is not an artifact of the mean-field approach. The slave-spin magnetization can serve as an order parameter for the metal-insulator transition.

Within the DIA, the optimal hybridization with the auxiliary bath site at zero one-particle energy can be seen as an analogous order parameter. There is, however, an important difference since the corresponding symmetry group is U(1): Physical observables depend on only and are therefore invariant under a phase change of . The fact that in the metallic phase can thus be interpreted as a spontaneous breaking of a local U(1) gauge symmetry which is restored, i.e., , in the Mott insulator at . The choice for the metallic phase just fixes the gauge.

Our numerical results for the quench dynamics at indicate that there is a transition from a symmetry broken, state at time to a symmetric state, for , i.e., a time-dependent Mott transition – similar to that in the Gutzwiller approach. (Schiró and Fabrizio, 2011) In the Gutzwiller approach and in the two-site DIA as well, the final state that is reached for is not the thermal state. Namely, a full two-site DIA equilibrium calculation for would give at any finite temperature, and the corresponding thermal state has a lower grand potential than that of the Hubbard-I-like state which is obtained in the equilibrium calculation by ad hoc setting (which is always a stationary point of the grand potential). Hence, within the two-site DIA restoring the local U(1) gauge symmetry in the time-dependent Mott transition necessarily implies that the final state is nonthermal. Within the DMFT, in contrast, the final state, which is (rapidly) reached after a quench to , is a thermal state with a high temperature (more than an order of magnitude higher than ), (Eckstein et al., 2009) which does not break the symmetry. (Žitko and Fabrizio, 2015) One may speculate that within the DIA the hybridization of the zero-energy mode with the impurity also vanishes in the thermal state at if more and more bath sites are added, so that a rapid decoupling of this mode at would not be at odds with rapid thermalization.

### iv.2 Ramps of the interaction

The previous discussion provokes the question whether the dynamical Mott transition and the conventional equilibrium Mott transition can be smoothly connected to each other. Since the dynamical transition occurs at a much weaker interaction, it is not at all obvious whether the two phenomena are related at all. One route to study this question is to consider a ramp with a finite duration rather than an instantaneous quench of the interaction as has been done using the Gutzwiller approach in Ref. Sandri et al., 2012. Ramping the interaction in a short time from to will make contact to the results found for a sudden quench. On the other hand, for ramps with infinite duration, i.e., if the interaction is changed adiabatically rather than suddenly, the system evolves along paths within the equilibrium phase diagram and will cross the line of equilibrium transitions (see Fig. 3). In fact, assuming that there is a critical interaction for any ramp time at all, one should expect that, with increasing , the critical interaction crosses over from (sudden quench) to (), since starting from a zero-temperature initial state, an adiabatic process will result in a zero-temperature final state.

To test our expectation we therefore consider a sequence of (cosine-shaped) ramps with different duration , see Eq. (6). Here we are limited to finite propagation times, , for practical reasons. Nevertheless, this allows us to study the relevant critical behavior for ramp times up to .

We begin the discussion for ramps of different duration to the same final interaction , starting from the same initial state that has also been considered for the quenches discussed in Sec. IV.1. As a measure of adiabaticity of the process we compare the total energy as a function of the interaction during the ramp with the corresponding equilibrium result. This is shown in Fig. 8 where the time-dependent value of the total energy during the ramp is plotted against the instantaneous value of the interaction. The resulting function can be compared with the equilibrium total energy (dashed line) as well as with the total energy after a sudden quench (straight line). With increasing ramp duration the curves converge to the equilibrium result, i.e., to the -dependence of the ground-state energy. For ramp times , the process is almost perfectly adiabatic.

As is discussed in the context of the (quantum) Kibble-Zurek mechanism, (Dziarmaga, 2010; Francuz et al., 2016) the dependence of the energy difference on the ramp time should asymptotically follow an inverse power law. Figure 8 demonstrates that the excitation energy does decrease with increasing . Here, one would expect the mean-field exponent, i.e., (see Ref. Sandri et al., 2012 for a discussion), and, roughly, our data are in fact consistent with this expectation. To extract a reliable value for the exponent, however, calculations for much longer ramp times would be necessary.

Independent of the ramp duration, we find essentially the same distinction between a weak- and a strong-coupling regime that has been discussed for the case of a sudden quench. The two regimes are sharply separated by a critical interaction which depends on the ramp duration (see discussion below). With increasing ramp time the dynamics becomes more well behaved in the sense that energy conservation becomes almost perfect for weak final interactions and is strongly improved in the strong-coupling case. As an example, in Fig. 9, we show the time dependencies of the optimal hybridization, of the double occupancy and of the total energy for different values of which are located below, right at, and above the dynamical critical interaction for a ramp with . Note that the process is clearly nonadiabatic. For we find that there are almost no oscillations of the optimal hybridization parameter after the ramp is completed. The same holds true for the double occupancy. Its (time-averaged) value after the ramp only slightly increases with increasing final interaction after having reached a minimum (close to zero) right at the critical point. This already indicates proximity to an adiabatic process where the double occupancy would just follow its equilibrium value, i.e., where it would monotonically decrease with increasing final interactions [see Fig. 2(b)].

We also obtain reasonable results for the time-dependent momentum distribution and, contrary to the study of quench dynamics, can therefore more comprehensively focus on the question of thermalization. In Fig. 10, we show three different examples for the final-state dynamics of , exemplary for the weak- [Fig. (a)a] and for the strong-coupling case [Fig. (b)b] as well as for [Fig. (c)c]. We again consider ramps with . The initial state is characterized by a sharp Fermi-surface discontinuity, which is slightly washed out by the finite temperature (). The final state that is reached in the long-time limit either shows a sharp jump of at the Fermi surface (in case of ) or collapse-and-revival oscillations (). This is very similar to the DMFT results in the quench case.(Eckstein et al., 2009) Right at the critical point () we also find fast thermalization toward a hot thermal distribution immediately after the ramp is completed. Comparing with Hubbard-I equilibrium calculations, we find an effective temperature of , see Figs. (c)c and (d)d. Note that this is somewhat lower than the effective temperature that had been obtained for the quench ().

Let us finally come back to the original motivation to study ramps of the interaction. We consider ramps with various durations bridging the limit of an instantaneous quench and the adiabatic limit . For each ramp time, we have performed a series of calculations with different to extract the respective value of the dynamical critical interaction . The latter is indeed well defined in the whole regime. Its dependence on the ramp time for is shown in Fig. 11.

monotonically increases with and seems to approach the value of the critical interaction for the zero-temperature Mott transition, as obtained by the two-site DIA (cf. Sec. III). However, the convergence turns out to be very slow. We also cannot fully exclude that the low but nonzero initial-state temperature has some effect on the result expected for and that, even for a perfectly adiabatic process, the final-state effective temperature becomes nonzero which would imply that converges to a somewhat lower value. Nevertheless, the results indeed clearly indicate that the dynamical Mott transition and the equilibrium Mott transition are related phenomena which are smoothly connected – at least within the two-site DIA. The same conclusion can be drawn from the results of the Gutzwiller calculations (Sandri et al., 2012) which, however, show additional oscillations of the critical interaction when increasing the ramp duration. This effect is absent in the two-site DIA.

### iv.3 Discussion of the method

To conclude this section, let us contrast our approach with Hamiltonian-based methods which strive to solve the effective impurity model of nonequilibrium DMFT exactly by mapping it onto a single-impurity Anderson model with a finite number of bath sites,(Gramsch et al., 2013) which can then be treated numerically. The number of bath orbitals is systematically increased until the properties of the DMFT bath are accurately represented. In the current implementation, the number of bath sites needed for an accurate representation of the DMFT bath scales roughly linearly with the maximum simulation time (and it also weakly depends on the parameter regime). This limits simulations to short times. Different implementations have been put forward to solve the finite impurity model, using exact-diagonalization techniques, (Gramsch et al., 2013) the multiconfiguration time-dependent Hartree method, (Balzer et al., 2015) as well as an approach based on the matrix-product state representation.(Wolf et al., 2014) With exact-diagonalization methods(Gramsch et al., 2013) propagation up to inverse hoppings has been possible by providing bath sites at weak interactions, whereas using matrix-product states, (Wolf et al., 2014) () could be reached with () sites at strong (weak) interactions, i.e., the Hamiltonian based solvers are currently aiming at a numerical exact solution at short times.

In contrast to this, the self-energy functional approach maps the original lattice-fermion problem onto an auxiliary model with a fixed, small number of bath sites and, in the case of the dynamical impurity approximation, a single correlated site. Rather than aiming at an exact solution of the nonequilibrium DMFT equations, the SFT provides an independent variational scheme to determine the time-dependent one-particle parameters of the reference system which only in the limit of an infinite number of bath sites recovers the DMFT. Formally, a qualitatively correct time evolution on much longer propagation times is thus possible with a very small number of bath sites. The present study has in fact shown that even with the most simple reference system (with a single bath site only) one can make close contact with full DMFT results. The agreement between the two-site DIA and full DMFT is qualitatively satisfying and close to the critical point for the (dynamical) Mott transition even quantitative. This demonstrates that much of the essential physics can be captured with a single time-dependent variational parameter.

The general framework of the SFT ensures that variational approximations are conserving with respect to the particle number and spin. The possible violation of energy conservation, however, must be seen as a major drawback of the present implementation of the DIA. Ways to overcome this problem have been discussed in Ref. Hofmann et al., 2013. Here, we could show that energy conservation is in fact violated but that, on the other hand, this violation is moderate in the weak-coupling limit after a quantum quench and even for strong interactions does not generally invalidate the results which still agree qualitatively with DMFT. Furthermore, if the dynamics is initiated by ramping the interaction, energy conservation is respected to a much higher degree.

## V Summary

The nonequilibrium extension of the self-energy functional theory (SFT) has been applied to study the real-time dynamics of the Fermi-Hubbard model initiated by sudden quenches and by ramps of the interaction parameter, starting from the noninteracting limit to different final values . As the simplest nontrivial approximation within the SFT, which provides a local trial self-energy, we have employed the two-site dynamical impurity approximation (DIA). The dynamical Mott metal-insulator transition represents an ideal first test case for this method.

We have studied the dynamical Mott transition by systematically tracing the time evolutions of the double occupancy, the total energy, the momentum distribution, and the optimal hybridization parameter of the reference system for different . All quantities are found to exhibit distinct response behavior in the weak- and in the strong-coupling regime. Within the two-site DIA these regimes are separated by a sharp critical interaction at which the low-energy bath site decouples from the correlated site in the course of time. By analyzing long-time averages and comparing these with thermal results, we have found fast thermalization for quenches to and clear indications for prethermalization in both, the weak- and the strong-coupling regime. In all relevant aspects, this is in surprisingly good qualitative agreement with a previous nonequilibrium DMFT study. (Eckstein et al., 2009) This also holds for the numerical value of the dynamical critical interaction which turns out to be roughly a factor of two smaller than the critical value for the equilibrium Mott transition at zero temperature.

Comparing results for ramps of different duration, we could trace the critical behavior in the whole range from a sudden quench to the limit of an adiabatic quasistatic process. We found that there is a well-defined critical interaction in all cases which monotonically increases with the ramp time and which converges to the zero-temperature critical point in the adiabatic limit. Qualitatively, this agrees well with the predictions of the Gutzwiller approach. (Schiró and Fabrizio, 2011; Sandri et al., 2012)

In view of the comparatively simple but successful two-site dynamical impurity approximation it appears very promising to improve the approach by considering more complex reference systems. An improved study of the mean-field dynamics using a reference system with more bath degrees of freedom suggests itself. Furthermore, cluster approximations generating non-local trial self-energies appear highly interesting. Both routes are computationally demanding (Hofmann et al., 2016) concerning both complexity but also stability, but on the other hand also very promising. Furthermore, the nonequilibrium self-energy functional theory is a completely general approach and can be applied to more complex systems beyond the single-band Hubbard model at half filling – similar to the equilibrium case. (Potthoff, 2012, 2014) We expect that the method can serve as a highly useful tool to uncover and understand intriguing phenomena of strongly-correlated many-body systems out of equilibrium.

###### Acknowledgements.

We would like to thank Christian Gramsch for numerous helpful discussions. Support of this work by the Deutsche Forschungsgemeinschaft (Germany) within the Sonderforschungsbereich 925 (projects B5 and B4) is gratefully acknowledged.## References

- Imada et al. (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998)
- Jördens et al. (2008) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature 455, 204 (2008)
- Schneider et al. (2008) U. Schneider, L. Hackermüller, S. Will, T. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008)
- Giamarchi (2004) T. Giamarchi, Quantum Physics in One Dimension (Oxford University Press, 2004)
- Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008)
- Lewenstein et al. (2012) M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems (Oxford University Press, Oxford, 2012)
- Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415, 39 (2002)
- Strohmaier et al. (2010) N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Altman, and E. Demler, Phys. Rev. Lett. 104, 080401 (2010)
- Iwai et al. (2003) S. Iwai, M. Ono, A. Maeda, H. Matsuzaki, H. Kishida, H. Okamoto, and Y. Tokura, Phys. Rev. Lett. 91, 057401 (2003)
- Perfetti et al. (2006) L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf, Phys. Rev. Lett. 97, 067402 (2006)
- Wall et al. (2011) S. Wall, D. Brida, S. R. Clark, H. P. Ehrke, D. Jaksch, A. Ardavan, S. Bonora, H. Uemura, Y. Takahashi, T. Hasegawa, H. Okamoto, G. Cerullo, and A. Cavalleri, Nat. Phys. 7, 114 (2011)
- Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991)
- Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994)
- Rigol et al. (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008)
- Dziarmaga (2010) J. Dziarmaga, Adv. Phys. 59, 1063 (2010)
- Polkovnikov et al. (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011)
- Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989)
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996)
- Schmidt and Monien (2002) P. Schmidt and H. Monien, cond-mat/0202046, (2002)
- Freericks et al. (2006) J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006)
- Aoki et al. (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014)
- Moeckel and Kehrein (2008) M. Moeckel and S. Kehrein, Phys. Rev. Lett. 100, 175702 (2008)
- Eckstein et al. (2009) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. 103, 056403 (2009)
- Schiró and Fabrizio (2010) M. Schiró and M. Fabrizio, Phys. Rev. Lett. 105, 076401 (2010)
- Schiró and Fabrizio (2011) M. Schiró and M. Fabrizio, Phys. Rev. B 83, 165105 (2011)
- Sandri et al. (2012) M. Sandri, M. Schiró, and M. Fabrizio, Phys. Rev. B 86, 075122 (2012)
- Hamerla and Uhrig (2013) S. A. Hamerla and G. S. Uhrig, Phys. Rev. B 87, 064304 (2013)
- Hamerla and Uhrig (2014) S. A. Hamerla and G. S. Uhrig, Phys. Rev. B 89, 104301 (2014)
- Werner et al. (2009) P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 035320 (2009)
- Gull et al. (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011)
- Eckstein and Werner (2010) M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010)
- Eckstein and Werner (2011) M. Eckstein and P. Werner, Phys. Rev. Lett. 107, 186406 (2011)
- Tsuji et al. (2013) N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. 110, 136404 (2013)
- Caffarel and Krauth (1994) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994)
- Gramsch et al. (2013) C. Gramsch, K. Balzer, M. Eckstein, and M. Kollar, Phys. Rev. B 88, 235106 (2013)
- Balzer et al. (2015) K. Balzer, Z. Li, O. Vendrell, and M. Eckstein, Phys. Rev. B 91, 045136 (2015)
- Wolf et al. (2014) F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Phys. Rev. B 90, 235131 (2014)
- Potthoff (2003a) M. Potthoff, Eur. Phys. J. B 32, 429 (2003a)
- Hofmann et al. (2013) F. Hofmann, M. Eckstein, E. Arrigoni, and M. Potthoff, Phys. Rev. B 88, 165124 (2013)
- Potthoff (2012) M. Potthoff, in Strongly Correlated Systems: Theoretical Methods, Springer Series in Solid-State Sciences, Vol. 171, edited by A. Avella and F. Mancini (Springer, Berlin, 2012)
- Potthoff (2014) M. Potthoff, in DMFT at 25: Infinite Dimensions, Modeling and Simulation, Vol. 4, edited by D. V. Eva Pavarini, Erik Koch and A. Lichtenstein (Forschungszentrum Jülich, 2014)
- Potthoff (2003b) M. Potthoff, Eur. Phys. J. B 36, 335 (2003b)
- Požgajčić (2004) K. Požgajčić, cond-mat/0407172, (2004)
- Koga et al. (2005) A. Koga, K. Inaba, and N. Kawakami, Progress of Theoretical Physics Supplement 160, 253 (2005)
- Inaba et al. (2005) K. Inaba, A. Koga, S. I. Suga, and N. Kawakami, Phys. Rev. B 72, 085112 (2005)
- Eckstein et al. (2007) M. Eckstein, M. Kollar, M. Potthoff, and D. Vollhardt, Phys. Rev. B 75, 125103 (2007)
- Hofmann et al. (2016) F. Hofmann, M. Eckstein, and M. Potthoff, Journal of Physics: Conference Series 696, 012002 (2016)
- Wagner (1991) M. Wagner, Phys. Rev. B 44, 6104 (1991)
- van Leeuwen et al. (2006) R. van Leeuwen, N. E. Dahlen, G. Stefanucci, C. O. Almbladh, and U. von Barth, in Time-Dependent Density Functional Theory, Vol. 706, edited by M. A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K. Burke, and E. K. U. Gross (Springer, Berlin Heidelberg, 2006)
- Rammer (2007) J. Rammer, Quantum field theory of nonequilibrium states (Cambridge University Press, 2007)
- () Concerning the system size, the equilibrium results are converged up to the sixth digit already for 40 sites. The corresponding error is three orders of magnitude smaller than the error resulting from the finite time-stepping (see below). Consequently, there is in fact no recognizable impact on the real-time dynamics when enlarging the system. Applying the trapezoidal rule for all integrations involved on the Keldysh branch, we find an error scaling quadratically with the size of the time step (the error is defined as , i.e., the maximum-norm difference of the optimal parameters obtained from calculations with time steps and ). In particular, the error of all our calculations is of the order of , i.e., plots for results with a smaller time step could hardly be distinguished from the results shown.
- Kelley (1987) C. T. Kelley, Solving nonlinear equations with Newtonâs method, Fundamentals of algorithms (Society for Industrial and Applied Mathematics, 1987)
- Broyden (1965) C. G. Broyden, Math. Comp. 19, 577 (1965)
- Strand et al. (2011) H. U. R. Strand, A. Sabashvili, M. Granath, B. Hellsing, and S. Östlund, Phys. Rev. B 83, 205136 (2011)
- Bulla (1999) R. Bulla, Phys. Rev. Lett. 83, 136 (1999)
- Lange (1998) E. Lange, Modern Physics Letters B 12, 915 (1998)
- Hubbard (1963) J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963)
- Žitko and Fabrizio (2015) R. Žitko and M. Fabrizio, Phys. Rev. B 91, 245130 (2015)
- de’Medici et al. (2005) L. de’Medici, A. Georges, and S. Biermann, Phys. Rev. B 72, 205124 (2005)
- Rüegg et al. (2010) A. Rüegg, S. D. Huber, and M. Sigrist, Phys. Rev. B 81, 155118 (2010)
- Francuz et al. (2016) A. Francuz, J. Dziarmaga, B. Gardas, and W. H. Zurek, Phys. Rev. B 93, 075134 (2016)