# Dynamic inhomogeneities and phase separation after quantum quenches in strongly correlated systems

## Abstract

We present a Gutzwiller von Neumann dynamics (GvND) method for simulating equilibrium and nonequilibrium phenomena in strongly correlated electron systems. Our approach is a real-space formulation of the time-dependent Gutzwiller approximation method. Applying the GvND method to simulate interaction quenches in the Peierls-Hubbard model, we demonstrate the amplification of initial inhomogeneities, which in turn results in the collapse of quench-induced synchronized oscillation. Moreover, we find a dynamical phase transition separating two quasi-stationary regimes with rather distinct spatial distributions of physical quantities after the collapsed oscillation. In particular, in the strong-coupling regime, the system exhibits a dynamic phase separation in the quasi-stationary state. Our results thus underscore the importance of spatial fluctuations in the nonequilibrium dynamics of strongly correlated systems.

The recent enormous theoretical interest in nonequilibrium dynamics of correlated materials is partly driven by the impressive progress in experimental techniques for controlling and probing such systems (1); (2); (3). In particular, developments of time-resolved measurement techniques in solids and cold atoms now allow one to study dynamical phase transitions far from equilibrium on the microscopic time scale of correlated systems (4); (5). Theoretically, the quantum quench setup provides an idealized platform to investigate intrinsic out-of-equilibrium dynamics related to strong electron correlation (6); (7); (12); (8); (13); (9); (11); (10); (14). One crucial question here is whether the system eventually thermalize and reach a new equilibrium (15). After a quench to a large interaction parameter, the system exhibits characteristic collapse-and-revival oscillations which eventually fade out in the long time. But in some cases, the system is trapped in a nonthermal quasi-stationary state for very long times (13); (12); (14).

While there has been considerable progress in our understanding of nonequilibrium dynamics in one dimension (16); (17), investigation of quantum quench in high dimensional systems has only begun recently. Here we are mainly interested in interaction quench in fermionic systems, of which the single-band Hubbard model is a canonical example. The nonequilibrium dynamics of a Fermi sea after a sudden switch-on of the Hubbard repulsion in infinite dimensions has been studied in a pioneering work (10) using the time-dependent dynamical mean-field theory (DMFT) (18); (19). The results clearly indicate the existence of a dynamical phase transition separating two distinct out-of-equilibrium regimes depending on the Hubbard parameter after quench. However, up to date, most theoretical studies of quench dynamics ignores the spatial inhomogeneity and fluctuations, which seem to be ubiquitous in nonequilibrium dynamics of many-body systems. A famous example is the formation of topological defects, described by the Kibble-Zurek mechanism (20); (21), when a system is driven across a continuous phase transition. The DMFT method, even with its cluster or real-space generalization (22); (23); (24), is still very limited in its treatment of complex spatial structures due to its heavy computational cost. It is thus highly desirable to develop an approximate yet efficient approach for large-scale real-space simulations of correlated systems. In this regard, the recently developed time-dependent Gutzwiller approximation (TDGA) (25); (26); (27); (28) provides such an efficient alternative to the more accurate nonequilibrium DMFT. More importantly, TDGA was shown to capture many nontrivial effects observed in DMFT (25), such as the existence of a dynamical phase transition, when applied to quantum quenches in the Hubbard model.

In this paper, we present a Gutzwiller-von Neumann dynamics (GvND) formulation for real-space dynamical simulations of correlated systems. In this approach, the time evolution of the Gutzwiller variational parameters is coupled to the von Neumann equation governing the dynamics of many-electron wavefunction. We apply the GvND method to investigate the interaction quench in a triangular-lattice Hubbard model with electron-phonon coupling and uncover several dramatic effects caused by the dynamical lattice degrees of freedom. First, the phase-locked oscillation of, e.g. double occupation, caused by the sudden turn-on of interaction is disrupted by the development of inhomogeneous lattice deformations. Moreover, we find two regimes of quasi-steady states with rather distinct spatial distributions of physical quantities after the collapsed oscillation. In particular, for large , the system exhibits a spontaneous phase separation after the quantum quench.

We consider the Hubbard model on a deformable triangular lattice described by the Hamiltonian

(1) | |||||

Here is the creation operator of electron with spin at site-, is the electron number operator, is the Hubbard repulsion parameter, is a elastic constant, is the momentum operator, and is the mass of the atom. denotes the displacement vector of the -th site, i.e. , and is a unit vector pointing from site- to . We consider the following dependence of hopping integral on the displacements

(2) |

where is the bare hopping constant, and is the electron-phonon coupling. The Hamiltonian (1) with given by Eq. (2) is also called the Peierls-Hubbard (PH) model (29); (30). The 1D version of the PH model is the famous Su-Schrieffer-Heeger model (31) with Hubbard interaction. Earlier interest on the 2D square-lattice PH model is motivated by the interest in high- superconductors (32); (33); (34). The PH model has served as the basic platform for investigating the interplay of Peierls instability and electron correlation. Since our main interest here is correlation related nonequilibrium phenomena without spontaneous symmetry-breaking, we focus on half-filled triangular-lattice PH model in which the lack of Fermi surface nesting prevents the system from weak-coupling instability either in the elastic or electronic subsystems. The phonons here are mainly used as a way to introduce dynamic inhomogeneities to the system.

We first briefly review the TDGA method (25); (26), which can be derived from the Dirac-Frenkel variational principle (35); (36). The Gutzwiller wavefunction is expressed as , where is a time-dependent Slater determinant constructed from a renormalized Hamiltonian. is a time-varying Gutzwiller operator, which can be parameterized by a set of variational matrices in the basis of local Hilbert space. This GA formulation is intimately related to the slave-boson approach (37); (38). In a sense, the matrix elements can be viewed as amplitudes of the slave boson coherent state (39); (40). Importantly, expectation value of local operator is now expressed as a trace: (39). The time evolution of is governed by the Schrödinger equation , while the variational matrix obeys the equation of motion: , where is the on-site Coulomb interaction expressed in the local basis of , and the electron binding energy is . Here we have introduced the reduced electron density matrix , and is the Gutzwiller renormalization factor (39)

(3) |

In our real-space formulation, the dynamical equation for the slave bosons reads

(4) | |||||

where and are the number and double-occupation operators, respectively, in the local basis, and are effective chemical potentials determined from the energy conservation condition. The above Eq. (4) indicates that the evolution of depends on the electron density matrix, whose dynamics is governed by the von Neumann equation , or explicitly:

(5) | |||

Here we have included an diagonal on-site potential to the Hubbard Hamiltonian. Eqs. (4) and (5) comprise a complete set of differential equations for the dynamics of correlated electrons in the real-space TDGA framework. We note in passing that the above formulation can be directly generalized to multi-orbital Hubbard-like models by treating , , as the spin-orbital indices.

In the presence of dynamical lattice degrees of freedom, these two equations are further coupled to the equation of motion of the displacement field:

(6) | |||

This equation is basically the Newton equation for the atoms, and the above formulation can be viewed as a non-adiabatic generalization of the Gutzwiller molecular dynamics (MD) method discussed in Ref. (41). While the motion of atoms is confined in the vicinity of the lattice points here, Ref. (41) considers full MD dynamics for strongly correlated systems in the liquid phase in the adiabatic limit, which assumes that electrons quickly relax to the equilibrium state of the instantaneous renormalized Hamiltonian. Consequently, both and are determined self-consistently from the ionic configuration at each time step (41). In our GvND formulation here, both the slave bosons and the electrons follow their own dynamics in the more general non-adiabatic and non-equilibrium situation.

We next apply the above GvND dynamics method to simulate interaction quench in the triangular-lattice PH model. For single-band Hubbard system, the most general form of is a diagonal matrix with diagonal elements for empty site, and , for singly occupied state with a spin-up or down electron, and finally for a doubly occupied site. In this basis, both the double-occupation and number operators are also diagonal: , , and . Here we focus on non-magnetic electronic states, and assume and , and . In this special case, we have for the chemical potentials.

In our interaction quench simulations, the system is initially in the ground state of a hal-filled PH model with . At , the Hubbard repulsion is suddenly turned on to . Fig. 1 shows the time dependence of the various physical quantities from simulations with and , where is the bandwidth. The parameters used in these simulations are , , and mass ; energy and inverse time are measured in units of the nearest-neighbor hopping . After the Coulomb interaction is switched on, the spatial averaged double-occupation and the quasiparticle weight exhibits the characteristic oscillation in short time scales, as shown in the insets of Fig. 1(a) and (d). Such oscillations have been reported previously in DMFT (10) as well as TDGA (25) simulations of the Hubbard model.

Interestingly, we find that this oscillation only lasts up to a time , beyond which both the quasiparticle weight and double-occupation approach quasi-stationary value and , respectively, after a transient period. By investigating the behavior of individual sites, we find that both quantities continue to oscillate but with a site-dependent amplitude and frequency. This phenomenon can thus be viewed as the collapse of synchronized oscillations among individual sites. To understand what causes this collapse, we note that the collapse at coincides with the rapid rise of the standard deviations, and , of double-occupation and on-site density, indicating the development of a inhomogeneous configuration.

It is worth noting that the coherent oscillation after the quench is a unstable quasi-stationary state. We find that the collapse of the oscillation results from the amplification of initial disorder or inhomogeneity, no matter how small it is. In order to have a more controlled simulation, we introduce random displacements of the order of in the initial state. The electron contribution to the forces acting on an atom is given by the second term in Eq. (6). In the coherent oscillation regime, we can approximate the forces as , where is deviation from the density matrix of a uniform electronic state. In a perfectly ordered state, the atoms experience no forces coming from the electrons. However, in the presence of disorder, the rapid oscillation of the renormalization factors tends to amplify the effect of inhomogeneity , hence injecting energy to the phonons. The increased lattice distortion in turn enhances the inhomogeneity in density matrix. This positive feedback eventually leads to the collapse of the coherent oscillation. Indeed, detailed examination shows that after the quench the elastic energy grows exponentially until the time . After the collapse of the oscillation, both the potential and kinetic energy relax to quasi-stationary values, as demonstrated in Fig. 1(c) and (f).

An intriguing phenomenon in the interaction quenches of Hubbard model is the existence of dynamical phase transition that separates the weak and strong coupling regimes. This result was first reported from the DMFT simulation (10), and was later reproduced using the TDGA approach (25). Our GvND simulations also find a dynamical transition at . Fig. 2(a) and (b) show the time trace of the (spatial) averaged double occupation in the weak and strong coupling regimes, respectively. The period of the coherent oscillation (before the collapse), shown in Fig. 2(c) as a function of , diverges as approaches the critical value, consistent with previous result. Moreover, both and tend to zero when approaching the transition point, indicating an emergent dynamical Mott insulating state; see Fig. 2(d).

The dynamical transition at not only separates two distinct temporal behaviors as already suggested in previous studies (10); (25). Our GvND simulations reveal rather distinct spatial patterns in the quasi-stationary states after the collapse of the coherent oscillation. Fig. 3 shows the spatial profile of double-occupation and displacement field at a time for varying . In the weak-coupling regime (), the distribution of double-occupation is relatively uniform and is dominated by short-range fluctuations. Indeed, its probability distribution , shown in Fig. 4(a), exhibits a single peak with a width that increases as approaches the critical value. Moreover, the displacements clearly shows high-frequency spatial modulation; see Fig. 3(d). At , both and show fractal-like structures with various length scales, resembling the spatial distribution of order parameters in the vicinity of an equilibrium phase transition. The corresponding exhibits a pronounced peak at , consistent with the fact that at the dynamical transition point.

On the other hand, a rather distinct distribution is obtained in the strong-coupling regime as demonstrated in Fig. 3(c) and (f). In this regime where , domains with large double-occupation and on-site density are interspersed with Mott droplets, or regions with vanishing and . This picture of phase separation can also be clearly seen from the double-peak structure in the probability distribution of double-occupation in this regime shown in Fig. 4(b). Moreover, we find that the Mott droplets are characterized by reduced electron hoppings. This is consistent with the fact that large displacements are found to coincide with these insulating regions, where the hopping amplitude is reduced due to the Peierls coupling; see Eq. (2). This dynamic phase separation in the regime of large is reminiscent of the site-selective Mott transition observed in highly disordered Hubbard model (42); (43). It is found that the insulating phase of the Anderson-Hubbard model with strong on-site disorder is a highly inhomogeneous state with Mott insulating droplets interlaced with regions containing Anderson-localized quasiparticles (43).

To summarize, we have presented a Gutzwiller von Neumann equation method for simulating real-space nonequilibrium dynamics in strongly correlated electron systems. Applying the GvND method to study the interaction quench in Peierls-Hubbard model, we find that the quench-induced coherent oscillation is intrinsically unstable against any initial disorder. In fact, this instability against inhomogeneity exists even without the dynamical phonons. A similar collapse of the synchronized oscillation can be induced by a tiny on-site disorder as in the Anderson-Hubbard model (44). In the presence of dynamical lattice degrees of freedom, we further observe dynamically generated phase separation when the system is quenched to a strong-coupling regime. Such amplification of inhomogeneity has already been reported previously in a similar quench dynamics study of the Luttinger liquid (45). However, in the 1D case, the amplification of the initial inhomogeneity is attributed to fractionalized quasiparticles (45), which is a special feature of 1D phase. On the other hand, similar dynamical inhomogeneity has also been found in the interaction quenches of high dimensional superfluid (46). In our case, the highly inhomogeneous nonequilibrium state certainly results from the nontrivial interplay of correlations and electron-phonon couplings.

The fact that the elastic part of the Peierls-Hubbard model is a linear system might lead to prolonged post-quench quasi-stationary states since there is no energy exchange between the phonons. Although different phonon modes can talk to each other indirectly through the electrons, we expect inclusion of elastic nonlinearity could introduce additional relaxations and help with the thermalization. Such issues can be addressed with molecular dynamics simulations based on our GvND method, which will be left for future studies.

###### Acknowledgements.

The author would like to thank K. Barros, C. Batista, G. Kotliar, and H. Suwa for collaborations on related works and numerous insightful discussions. This work is partly supported by the Center for Materials Theory as a part of the Computational Materials Science (CMS) program, funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.### References

- A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
- J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nature Phys. 11, 124 (2015).
- T. Langen, R. Geiger, and J. Schmiedmayer, Ultracold atoms out of equilibrium, Annu. Rev. Condens. Matter 6, 201 (2015).
- L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf, Time Evolution of the Electronic Structure of -TaS through the Insulator-Metal Transition, Phys. Rev. Lett. 97, 067402 (2006).
- D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K. Haule, Electrodynamics of correlated electron materials, Rev. Mod. Phys. 83, 471 (2011).
- M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons, Phys. Rev. Lett. 98, 050405 (2007).
- M. A. Cazalilla, Effect of Suddenly Turning on Interactions in the Luttinger Model, Phys. Rev. Lett. 97, 156403 (2006).
- M. Kollar and M. Eckstein, Relaxation of a one-dimensional Mott insulator after an interaction quench, Phys. Rev. A 78, 013626 (2008).
- D. Rossini, A. Silva, G. Mussardo, and G. E. Santoro, Effective Thermal Dynamics Following a Quantum Quench in a Spin Chain, Phys. Rev. Lett. 102, 127204 (2009).
- M. Eckstein, M. Kollar, and P. Werner, Thermalization after an interaction quench in the Hubbard model, Phys. Rev. Lett. 103, 056403 (2009).
- M. Moeckel and S. Kehrein, Interaction Quench in the Hubbard Model, Phys. Rev. Lett. 100, 175702 (2008).
- C. Kollath, A. M. Läuchli, and E. Altman, Quench Dynamics and Nonequilibrium Phase Diagram of the Bose-Hubbard Model, Phys. Rev. Lett. 98, 180601 (2007).
- S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Strongly Correlated Fermions after a Quantum Quench, Phys. Rev. Lett. 98, 210405 (2007).
- N. Tsuji, M. Eckstein, and P. Werner, Nonthermal Antiferromagnetic Order and Nonequilibrium Criticality in the Hubbard Model, Phys. Rev. Lett. 110, 136404 (2013).
- M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
- S. R. White and A. E. Feiguin, Real-Time Evolution Using the Density Matrix Renormalization Group, Phys. Rev. Lett. 93, 076401 (2004).
- A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Time-dependent density-matrix renormalization-group using adaptive effective Hilbert spaces, J. Stat. Mech. P04005 (2004).
- J. K. Freericks, V. M. Turkowski, and V. Zlatić, Nonequilibrium Dynamical Mean-Field Theory, Phys. Rev. Lett. 97, 266408 (2006).
- H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
- T. W. B. Kibble, Topology of cosmic domains and strings, J. Phys. A: Math. Gen. 9, 1387 (1976).
- W. H. Zurek, Cosmological experiments in superfluid helium ?, Nature 317, 505 (1985).
- A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
- T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Quantum cluster theories, Rev. Mod. Phys. 77, 1027 (2005).
- G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
- M. Schiró and M. Fabrizio, Time-dependent mean-field theory for quench dynamics in correlated electron systems, Phys. Rev. Lett. 105, 076401 (2010).
- M. Schiró and M. Fabrizio, Quantum quenches in the Hubbard model: Time-dependent mean-field theory and the role of quantum fluctuations, Phys. Rev. B 83, 165105 (2011).
- M. Sandri and M. Fabrizio, Nonequilibrium dynamics in the antiferromagnetic Hubbard model, Phys. Rev. B 88, 165113 (2013).
- G. Seibold and J. Lorenzana, Time-Dependent Gutzwiller Approximation for the Hubbard Model, Phys. Rev. Lett. 86, 2605 (2005).
- S. Mazumdar and S. N. Dixit, Coulomb Effects on One-Dimensional Peierls Instability: The Peierls-Hubbard Model, Phys. Rev. Lett. 51, 292 (1983).
- J. E. Hirsch, Effect of Coulomb Interactions on the Peierls Instability, Phys. Rev. Lett. 51, 296 (1983).
- W. P. Su, J. R. Schrieffer, and A. J. Heeger, Soliton excitations in polyacetylene, Phys. Rev. B 22, 2099 (1980).
- S. Tang and J. E. Hirsch, Peierls instability in the two-dimensional half-filled Hubbard model, Phys. Rev. B 37, 9546 (1988).
- H. Fehske, M. Deeg, and H. Büttner, Two-dimensional Peierls-Hubbard model within the slave-boson approach, Phys. Rev. B 46, 3713 (1992).
- Q. Yuan and T. Kopp, Coexistence of the bond-order wave and antiferromagnetism in a two-dimensional hall-filled Peierls-Hubbard model, Phys. Rev. B 65, 085102 (2002).
- P. A. M. Dirac, Note on Exchange Phenomena in the Thomas Atom, Proc. Cambridge Philos. Soc. 26, 376 (1930).
- J. Frenkel, Wave Mechanics: Advanced General Theory (Clarendon, Oxford, 1934).
- G. Kotliar and A. E. Ruckenstein, New Functional Integral Approach to Strongly Correlated Fermi Systems: The Gutzwiller Approximation as a Saddle Point, Phys. Rev. Lett. 57, 1362 (1986).
- T. Li, P. Wölfle, and P.J. Hirschfeld, Spin-Rotation-Invariant Slave-Boson Approach to the Hubbard Model, Phys. Rev. B 40, 6817 (1989).
- N. Lanatá, Y. Yao, X. Deng, V. Dobrosavljević, and G. Kotliar, Slave Boson Theory of Orbital Differentiation with Crystal Field Effects: Application to UO, Phys. Rev. Lett. 118, 126401 (2017).
- M. Behrmann, A. I. Lichtenstein, M. I. Katsnelson, and F. Lechermann, Versatile approach to spin dynamics in correlated electron systems, Phys. Rev. B 94, 165120 (2016).
- G.-W. Chern, K. Barros, C. D. Batista, J. Kress, and G. Kotliar, Mott transition in a metallic liquid – Gutzwiller molecular dynamics simulations, Phys. Rev. Lett. 118, 226401 (2017).
- K. Byczuk, W. Hofstetter, and D. Vollhardt, Mott-Hubbard Transition versus Anderson Localization in Correlated Electron Systems with Disorder, Phys. Rev. Lett. 94, 056404 (2005).
- M. C. O. Aguiar, V. Dobrosavljević, E. Abrahams, and G. Kotliar, Critical Behavior at the Mott-Anderson Transition: A Typical-Medium Theory Perspective, Phys. Rev. Lett. 102, 156402 (2009).
- G.-W. Chern, unpublished.
- M. S. Foster, E. A. Yuzbashyan, and B. L. Altshuler, Quantum Quench in One Dimension: Coherent Inhomogeneity Amplification and “Supersolitons”, Phys. Rev. Lett. 105, 135701 (2010).
- M. Dzero, E. A. Yuzbashyan, and B. L. Altshuler, Cooper pair turbulence in atomic Fermi gases, Europhys. Lett. 85, 20004 (2009).