---
abstract: 'We present the first-ever microscopic dynamical simulation of the temperature-controlled Mott metal-insulator transition in the Hubbard model. By combining the efficient Gutzwiller method with molecular dynamics simulations, we demonstrate that the transformation from the correlated metal to the Mott insulator proceeds via the nucleation and growth of the Mott droplets. Moreover, the time evolution of the Mott volume fraction is found to follow a universal transformation kinetics. We show that after an initial incubation period, the early stage of the phase transformation is characterized by a constant nucleation rate and an interface-controlled cluster growth mechanism, consistent with the classical theory developed by Kolmogorov, Johnson, Mehl, and Avrami. This is followed by a novel intermediate stage of accelerated phase transformation that is significantly different from the prediction of the classical theory. Morevoer, the cluster-growth dynamics in this stage exhibits an unexpected avalanche behavior, similar to the Barkhausen noise in magnetization dynamics, even in the absence of quenched disorder. Detailed structural characterization further uncovers a universal correlation function for the transient mixed-phase states of the Mott transition. The implications of our findings for the recent nano-imaging experiments on metal-insulator transition of correlated materials are also discussed.'
author:
- 'Gia-Wei Chern'
title: Kinetics of thermal Mott transitions in the Hubbard model
---

Introduction {#sec:intro}
============

Mott metal-insulator transition (MIT) continues to be an important subject in modern condensed matter physics even after more than sixty years of study [@mott90; @imada98; @dobrosavljevic12]. In particular, the unusual properties of the metallic state near the MIT, sometimes also called the bad metal state [@emery95], have attracted considerable attention both experimentally and theoretically [@limelette03; @bartosch10; @itou17; @pustogow18; @terletska11; @sordi12; @vucicevic15; @jawad15; @lenz16; @imada05]. Conventionally, the effort has been focused on the thermodynamic, optical, magnetic, and transport properties in the vicinity of the Mott transition. These are usually achieved experimentally through bulk or macroscopic measurements. Coupled with the advances in many-body techniques such as the quantum Monte Carlo simulations [@leblanc15; @motta17; @kent18] and the dynamical mean-field theory (DMFT) [@georges96; @kotliar06; @gull11], significant progress has been made toward understanding the microscopic mechanisms and the quantum criticality of the MIT.

An important aspect that has received attention only recently is the [*complexity*]{} of MITs in real materials [@dagotto05]. For example, several near-field spectroscopy and nano-imaging experiments [@bonnell12; @atkin12; @liu17] have revealed highly inhomogeneous electronic states with complex nano-scale textures during the first-order metal-insulator transformation in correlated electron systems [@qazilbash07; @liu13; @madan15; @lupi10; @qzailbash11; @mcleod16; @stinson18; @mattoni16; @ronchi18; @preziosi18; @singer18]. In fact, MIT in most functional materials, such as vanadium dioxide, is [*discontinuous*]{} at ambient conditions. Understanding the structural and dynamical properties of the mesoscopic patterns during the first-order MIT is not only of fundamental interest, but also has important technological implications.

In the past decade, considerable experimental effort has been devoted to understand the properties of the heterogeneous electronic states in the vicinity of the Mott transition. For instance, the nano-scale phase separation in several temperature-driven MITs was shown to be consistent with the nucleation and percolation scenario [@qazilbash07; @lupi10; @mcleod16; @liu13], similar to that observed in the colossal magnetoresistant materials (CMR) [@uehara99; @fath99; @zhang02]. On the theoretical side, however, a [*microscopic*]{} modeling of these emergent dynamical textures is still lacking. The difficulty is partly due to the multiscale and multi-faceted nature of the MIT kinetics. On one hand, large-scale dynamical simulation is required in order to describe the complex structures of the mixed-phase states in the Mott transition. On the other hand, accurate modeling of the electron correlation effects requires advanced many-body techniques [@leblanc15; @motta17; @kent18; @georges96; @kotliar06; @gull11], most of which are computationally too expensive to be directly combined with large-scale simulations. One exception is the Gutzwiller approach [@gutzwiller63; @gutzwiller64; @gutzwiller65], which allows for computational efficiency similar to the Hartree-Fock methods or the density functional theory (DFT). Importantly, the Gutzwiller method captures essential correlation effects such as the bandwidth renormalization and electron localization [@brinkman70], and offers perhaps the only feasible way for large-scale real-space simulations of MIT.

The mesoscale electronic inhomogeneities associated with the MIT are reminiscent of the complex patterns observed in numerous first-order phase transitions of, e.g. binary alloys, polymer mixtures, and many other physical systems [@langer80; @gunton83; @puri09]. The study of the first-order transition kinetics has a long history [@binder87; @onuki02]. Microscopically, the kinetic Ising model with, e.g. Glauber-type dynamics, has been extensively used to simulate the phase transformation dynamics [@marro79; @phani80; @sadiq83; @stauffer82; @shneidman99; @shneidman99b; @wonczak00; @soisson00; @brendel05; @ryu10; @binder16]. The Ising, or the liquid-gas, transition is also often invoked to model the Mott MIT by, e.g. associating the gas and liquid phases with the metallic phase and Mott insulator, respectively [@limelette03; @bartosch10; @itou17; @liu16; @bar18; @papanikolaou08]. Several works have indicated that Mott transitions likely belong to the Ising universality class. One objective of this work is to investigate whether the kinetic Ising model can also account for the nucleation dynamics of the Mott transition.

The phase transformation dynamics is conventionally formulated in terms of the appropriate order parameters that characterize the phase transition [@binder87; @onuki02]. However, since the MIT is a transition between different transport behaviors and is not characterized by broken symmetries, there is no clear choice of the order parameter. For Mott transitions, however, making use of the analogy with the liquid-gas transition, one can introduce an order parameter based on, e.g. the density of doublons, or doubly occupied sites [@lee17]. Since the Gutzwiller method can be viewed as a mean-field theory of the Mott transition within the Kotliar-Ruckenstein slave-boson framework [@kotliar86], it offers a systematic construction of the relevant order parameters. Specifically, they correspond to the coherent state $\Phi$ of the slave bosons [@lanata12; @lanata15]. For example, $\Phi$ in the single-band Hubbard model is represented by the amplitudes $\Phi = (e, p_{ \uparrow}, p_{ \downarrow}, d)$, whose square give the probability of the empty, singly occupied (up and down spin), and doubly occupied states, respectively [@kotliar86]. It thus encodes information about the local electron density $n_\alpha = |p_\alpha|^2$ as well as the double occupancy $D = |d|^2$. This formulation will be particularly useful for MIT in multi-orbital systems, where a large number of local electron states and potential order-parameters have to be considered [@bunemann98; @lanata15; @lanata17].

Phenomenologically, the dynamics of the order-parameter field depends on the conservation laws as well as the relevance of other hydrodynamic variables [@hohenberg77]. The simplest dynamics for a non-conserved order parameter $\Phi(\mathbf r, t)$, as in the case of Mott MIT, is the time-dependent Ginzburg-Landau (TDGL) equation [@binder73; @valls90; @iwamatsu08] $$\begin{aligned}
    \label{eq:TDGL}
    \frac{\partial \Phi}{\partial t} = -\Gamma \frac{\partial \mathcal{E}}{\partial \Phi^\dagger} + \eta(\mathbf r, t).\end{aligned}$$ Here $\mathcal{E}[\Phi]$ is an effective free-energy functional, $\Gamma$ denotes the inverse damping coefficient, and $\eta(\mathbf r, t)$ is a Gaussian white noise. In the case of the kinetic Ising model with Glauber dynamics, the TDGL corresponds to the coarse-grained approximation of the mean-field master equation [@puri09], and is often used in simulations of liquid-gas or magnetic transitions with non-conserved kinetics [@binder87; @onuki02]. The stochastic $\eta$ term, which plays an important role in the nucleation process, comes from the thermal fluctuations of the microscopic variables that are “integrated” out in the above formulation. Although the TDGL has also been used in modeling the hysteresis dynamics of the MIT [@bar18], its microscopic foundation has yet to be carefully examined. Indeed, a formal derivation of the dynamical equation for the order parameter of the Mott transition requires integrating out the microscopic electron and lattice degrees of freedom. This is generally very difficult if not impossible, and phenomenological expressions are often used in the kinetics simulations.

In this paper, we perform the first-ever microscopic simulation of the thermal Mott transition in the Hubbard model by keeping both lattice and electron degrees of freedom. We employ the Langevin molecular dynamics (MD) method to simulate the lattice fluctuations, which provide the relaxation mechanism and the source of stochastic force for the order parameter field. To explicitly include the microscopic correlation effects, the electron degrees of freedom are [*dynamically*]{} integrated out during the MD simulation. This is achieved by solving the Hubbard model using the efficient Gutzwiller method at every MD time-step. We show that the temperature-driven Mott transition is described by the nucleation and growth scenario, consistent with several recent nano-imaging experiments. We further characterize the different stages of the transformation dynamics and the structure of the intermediate mixed-phase states.

The rest of the paper is organized as follows. Section \[sec:GMD\] presents the new quantum MD scheme based on the Gutzwiller method, and the Hubbard model on a deformable lattice to be studied in this work. Section \[sec:thermal-mott\] discusses the finite-temperature phase diagram of the Hubbard model and outlines the Gutzwiller MD simulations of the thermal Mott transition driven by temperature quenches. In Section \[sec:nucleation\], the evolution of the cluster-size distribution from our MD simulations is discussed within the classical theory of nucleation-and-growth kinetics [@becker35; @zeldovich43; @frenkel46]. Section \[sec:interface\] presents the numerical evidence of an interface-controlled cluster-growth mechanism at the early stage of the phase transformation. Comparison of the transition kinetics with the classical Kolmogorov, Johnson, Mehl, and Avrami model [@kolmogorov37; @avrami39; @avrami40; @johnson39] is also discussed. In Section \[sec:avalanche\], we discuss the novel intermediate stage of the transformation and the associated avalanche cluster-growth behavior. Section \[sec:structure\] presents the numerical structure factor and two-point correlation function of the mixed-phase states during the Mott transition. And finally, Section \[sec:outlook\] concludes the paper with an outlook of future directions.

Gutzwiller molecular dynamics {#sec:GMD}
=============================

We consider the half-filled Hubbard model on a compressible triangular lattice $$\begin{aligned}
    \label{eq:H}
    \mathcal{H} = \sum_{\langle ij \rangle, \,\alpha } t(\mathbf u_i - \mathbf u_j) \left(c^\dagger_{i,\alpha} c^{\;}_{j, \alpha}  +  c^\dagger_{j,\alpha} c^{\;}_{i, \alpha} \right) 
    + U \sum_i n_{i,\uparrow} n_{i,\downarrow} \nonumber \\
    + \frac{K_0}{2} \sum_i |\mathbf u_i|^2 + \frac{K_1}{2} \sum_{\langle ij \rangle} \left[\hat{\mathbf e}_{ij} \cdot (\mathbf u_j - \mathbf u_i) \right]^2 + \sum_i \frac{|\mathbf p_i|^2}{2m}. \quad\end{aligned}$$ Here $c^\dagger_{i, \alpha}$ is the creation operator of electron with spin $\alpha = \uparrow, \downarrow$ at site-$i$, $n_{i, \alpha} = c^\dagger_{i,\alpha} c^{\;}_{i, \alpha}$ is the electron number operator, $U$ is the Hubbard repulsion parameter, $K$ is a elastic constant, $\mathbf p_i$ is the momentum operator, $m$ is the mass of the atom, $\mathbf u_i$ denotes the displacement vector of the $i$-th site, i.e. $\mathbf r_i = \mathbf r^{(0)}_i + \mathbf u_i$, and $\hat{\mathbf e}_{ij}$ is a unit vector pointing from site-$i$ to $j$. We assume the following dependence of hopping integral on the displacements $$\begin{aligned}
    \label{eq:t_ij}
    t_{ij} = t(\mathbf u_i - \mathbf u_j) = t^{(0)}_{ij} \left[1 + g \, \hat{\mathbf e}_{ij} \cdot (\mathbf u_j - \mathbf u_i) \right],\end{aligned}$$ where $t^{(0)}_{ij}$ is the bare hopping constant, and $g$ is the electron-phonon coupling. The Hamiltonian (\[eq:H\]) with $t_{ij}$ given by Eq. (\[eq:t\_ij\]) is called the Peierls-Hubbard model [@mazumdar83; @hirsch83], which has served as the basic platform for studying the interplay between the Peierls instability and electron correlation [@tang88; @fehske92; @yuan02]. Here we choose the frustrated triangular lattice to avoid the unnecessary complexity due to spontaneous symmetry-breaking or structural transition. As discussed in Sec. \[sec:intro\], the lattice degrees of freedom here serve to introduce the stochastic dynamics for the order parameter of the Mott transition.

The experimental time scale of temperature-driven Mott transition, which is also the time scale $\tau_\Phi$ for the order-parameter dynamics, is usually much longer than the characteristic time of lattice fluctuations. On the other hand, the relaxation time of electrons $\tau_e$, which ranges from a few tens to hundred femtoseconds, is short compared with the lattice time scale $\tau_L$ of the picosecond order. This separation of time scales allows us to apply the Born-Oppenheimer (BO) approximation in our MD simulations. Indeed, BO approximation is often employed in [*ab inito*]{} or quantum MD method [@marx09; @plasienka17]. The state-of-the-art DFT-based quantum MD, however, cannot describe the Mott transition. In order to capture the electron correlation effects, here we employ a new quantum MD scheme in which the atomic forces are computed on-the-fly from the Gutzwiller solution of Hubbard-type models [@chern17]. Applying the Gutzwiller MD to a Hubbard liquid model, we have recently demonstrated for the first time a Mott MIT in an atomic liquid [@chern17].

The atomic displacements $\mathbf u_i$ in the finite-temperature MD simulations are governed by the Newton equation of motion with the Langevin thermostat [@allen89; @rapaport11] $$\begin{aligned}
    \label{eq:langevin}
    m \ddot{\mathbf u}_i = \mathbf f^l_{i} + \mathbf f^e_{i} -\gamma \dot{\mathbf u}_i + \bm \xi_i(t),\end{aligned}$$ where $\gamma$ is the damping coeficient, $\bm\xi_i(t)$ are stochastic forces with statistical properties $\langle \xi^a_i \rangle = 0$ and $\langle \xi^a_i(t) \xi^b_j(t') \rangle = 2 \gamma k_B T \delta_{ij} \delta_{ab} \delta(t - t')$, $a, b = x, y, z$ denote the cartesean components, and $T$ is the temperature of the reservoir. Standard velocity-Verlet algorithm [@allen89; @rapaport11] is used to integrate the above Langevin equation in our simulations. The atomic forces within the BO-MD framework is computed using the finite temperature generalization of the Hellmann-Feynman theorem $\mathbf f_{i} = - {\rm Tr}(\rho_G \, \partial \mathcal{H} / \partial \mathbf u_i)$. There are two contributions to the atomic forces. The elastic restoring force is $$\begin{aligned}
    \mathbf f^l_i = -K_0 \mathbf u_i -K_1\,\sum_j \!^{'} \, \hat{\mathbf e}_{ij}  \cdot  (\mathbf u_i - \mathbf u_j),\end{aligned}$$ where the summation is restricted to the nearest neighbors. The electronic force originates from the $\mathbf u$-dependence of the transfer integral $t_{ij}$, and is given by $$\begin{aligned}
    \label{eq:f_e}
    \mathbf f^e_{i} = g \sum_{j, \alpha} t_{ij}^{(0)} \hat{\mathbf e}_{ij} \left[ {\rm Tr}(\rho_G \, c^\dagger_{i,\alpha}c^{\;}_{j\alpha} ) + \mbox{c.c.}\right],\end{aligned}$$ where $\rho_G = \mathcal{P}_G \, \rho^{\rm qp} \, \mathcal{P}_G$ is the Gutzwiller electron density matrix [@sandri13], $\rho^{\rm qp}$ is the density matrix of quasi-particles, and $\mathcal{P} = \prod_i \mathcal{P}_i$ is the Gutzwiller projector. As in Gutzwiller’s original approach, the local projectors $\mathcal{P}_i$ are introduced to reduce the on-site double-occupancy variationally [@gutzwiller63; @gutzwiller64; @gutzwiller65]. In the modern formulation of the Gutzwiller approximation (GA), these projectors can be conveniently parameterized using the so-called $\Phi$-matrix [@lanata12], $\mathcal{P}_i \equiv \sum_{\alpha, \beta} \Phi_{i, \alpha\beta} /(\Pi_{i, \beta})^{-1/2} |\alpha\rangle \langle \beta|$, where $\Pi_{i, \beta}$ is the occupation probability of local atomic state $|\beta\rangle$ for the uncorrelated state. For single-band Hubbard model, the $\Phi$-matrix is diagonal $$\begin{aligned}
    \Phi_i = \left[ \begin{array}{cccc}
    e_i & 0 & 0 & 0 \\
    0 & p_{i, \uparrow} & 0 & 0 \\
    0 & 0 & p_{i, \downarrow} & 0 \\
    0 & 0 & 0 & d_i
    \end{array} \right],\end{aligned}$$ and its elements represent the empty, singly occupied (spin up and down), and doubly occupied electron state of the $i$-th atom [@lanata15; @lanata17]. As discussed in Sec. \[sec:intro\], they can also be viewed as the coherent-state amplitude of the slave bosons in the Kotliar-Ruckenstein theory, whose saddle-point solution corresponds to the GA.

To obtain the Gutzwiller density matrix $\rho_G$, one needs to self-consistently solve the quasi-particle density matrix $\rho^{\rm qp}$ and the SB amplitudes $\Phi_i$. For a given set of atomic displacements$\{\mathbf u_i \}$, which defines the tight-binding $t_{ij}$ through Eq. (\[eq:t\_ij\]), the Gutzwiller solution is obtained by minimizing the electron free energy $$\begin{aligned}
    \label{eq:free-energy}
    & & \mathcal{F} = {\rm Tr}(\rho^{\rm qp} \mathcal{H}^{\rm qp}) + U \sum_i {\rm Tr}(\Phi^\dagger_i n_{i, \uparrow} n_{i, \downarrow} \Phi_i) \nonumber \\ 
    & &  \quad +  \, T \, {\rm Tr}(\rho^{\rm qp} \ln \rho^{\rm qp}) + T \sum_i {\rm Tr}\left[ \Phi^\dagger_i \Phi^{\;}_i \ln\left(\Phi^\dagger_i \Phi^{\;}_i / \Pi_i\right) \right] \nonumber \\
     & & \quad + \, \sum_{i, \alpha} \mu_{i,\alpha} \left[ {\rm Tr}(\Phi^\dagger_i \, \hat{n}_{i, \alpha} \Phi^{\;}_i)  - {\rm Tr}(\rho^{\rm qp} \hat{n}_{i, \alpha})\right].\end{aligned}$$ The first line above denotes the total energy within the Gutzwiller approximation. The two terms in the second line correspond to quasi-particle entropy and the relative entropy between correlated distribution $\Phi^\dagger_i \Phi_i$ and the corresponding uncorrelated $\Pi_i$ [@wang10; @sandri13]. The atom-dependent Lagrangian multipliers $\mu_{i, \alpha}$ are introduced to enforce the Gutzwiller constraints in the third line, and $\hat{n}_{i\alpha} \equiv c^\dagger_{i\alpha} c^{\;}_{i\alpha}$ is the electron number operator.

The quasi-particle density matrix is given by the Boltzmann distribution $\rho^{\rm qp} = e^{-\beta \mathcal{H}^{\rm qp}} / Z^{\rm qp}$, where $\mathcal{H}^{\rm qp}$ is the renormalized tight-binding Hamiltonian $$\begin{aligned}
    \label{eq:H_renorm}
    \mathcal{H}^{\rm qp} &=&  \sum_{\langle ij \rangle, \alpha} \left( \mathcal{R}^{\;}_{i,\alpha} \mathcal{R}^*_{j,\alpha}\, t_{ij} c^\dagger_{i,\alpha} c^{\;}_{j, \alpha} + \mbox{h.c.} \right) \nonumber \\
     & & \qquad - \sum_{i, \alpha} \mu_{i, \alpha}\, c^\dagger_{i, \alpha} c^{\;}_{i, \alpha}, \end{aligned}$$ where $t_{ij}$ depends on the displacements through Eq. (\[eq:t\_ij\]), and the renormalization factors are [@lanata12] $$\begin{aligned}
    \label{eq:R}
    \mathcal{R}_{i, \alpha} = \frac{ {\rm Tr}\left(\Phi^\dagger_i c^\dagger_{i, \alpha} \Phi^{\;}_i c^{\;}_{i, \alpha} \right) }{  \sqrt{n_{i, \alpha} (1 - n_{i, \alpha})} }.\end{aligned}$$ and $n_{i,\alpha} = {\rm Tr}(\Phi^\dagger_i \, \hat{n}_{i, \alpha} \Phi^{\;}_i)$ is the quasi-particle number at the $i$-th atom.

The minimization of $\mathcal{F}$ with respect to $\Phi_i$, subject to constraint ${\rm Tr}(\Phi^\dagger_i \Phi^{\;}_i) = 1$ can be recast into a nonlinear eigenvalue problem $\mathcal{H}^{\rm SB}_i \Phi_i = \epsilon_i \Phi_i$ for each atom [@lanata12]. Since $\Phi_i$ is diagonal in the case of single-band Hubbard model, we use the same notation $\Phi_i = (e_i, p_{i, \uparrow}, p_{i, \downarrow}, d_i)$ to denote the wavefunction or the coherent state of the slave bosons. The embedded Hamiltonian has the form $$\begin{aligned}
    \label{eq:H_SB}
    \mathcal{H}^{\rm SB}_i &=& \sum_{\alpha}  \frac{  \left(\Delta_{i, \alpha} \mathcal{M}_{\alpha} + \Delta^*_{i, \alpha} \mathcal{M}^{\dagger}_{\alpha} \right)}{\sqrt{n_{i\alpha} (1-n_{i\alpha})}} + U \,\mathcal{U} \nonumber \\
      & & + \sum_\alpha  \left( \Lambda_{i\alpha} + \mu_{i\alpha} \right) \mathcal{N}_\alpha + T \Theta_i,\end{aligned}$$ where $\Delta_{i, \alpha}$ encodes the local electron bonding $$\begin{aligned}
    \label{eq:Delta}
    \Delta_{i, \alpha} = \sum_j t_{ij} \mathcal{R}^*_{j,\alpha} {\rm Tr}(\rho^{\rm qp}  c^\dagger_{i\alpha} c^{\;}_{j \alpha} ).\end{aligned}$$ Explicit expressions for the various matrices $\mathcal{M}_\alpha$, $\mathcal{N}_\alpha$, $\mathcal{U}$ and $\Theta_\alpha$, as well as expression for the coefficient $\Lambda_{i, \alpha}$ are given in Appendix \[sec:GA\]. The Hamiltonian $\mathcal{H}^{\rm SB}$ essentially describes the competition between electron bonding favored by the first $\mathcal{M}$ term, and the Coulomb repulsion represented by the second $\mathcal{U}$ term. The Lagrangian multiplier $\mu_{i\alpha}$ appears in the third term acts as a chemical potential in order to satisfy the self-consistent Gutzwiller constraint.

The above two steps, namely solving the density matrix $\rho^{\rm qp}$ and variational matrices $\Phi_i$, have to be iterated until convergence is reached at every MD time-step. In particular, the Gutzwiller constraint $$\begin{aligned}
    {\rm Tr}(\Phi^\dagger_i \Phi^{\;}_i \hat{n}_{i, \sigma})  = {\rm Tr}(\rho^{\rm qp} \hat{n}_{i, \alpha}) = n_{i,\alpha},\end{aligned}$$ has to be satisfied. The bottleneck of the MD simulation here is solving the renormalized tight-binding Hamiltonian, which is required at every Gutzwiller iteration in every time-step. In our simulations, we adopt the gradient-based probing technique combined with the kernel polynomial method (KPM) [@barros13; @wang18; @chern18]. This allows for an efficient estimation of the electron density matrix $\rho^{\rm qp}$ at a computational cost that scales linearly with system size. All simulations discussed below were performed on triangular lattice of unprecedented size of $N = 150^2$ sites.

In this work we focus on the intrinsic Mott transition in the absence of magnetic ordering. We consider only non-magnetic Gutzwiller solutions, $p_{i, \uparrow} = p_{i, \downarrow} = p_i$, and $n_{i, \uparrow} = n_{i, \downarrow} = n_i$, in all simulations discussed below. The following model parameters are used. $t^{(0)}_{\rm nn1} = -1$, $g = 0.55$, $K_0 = K_1 = 0.35$, $m = 2$, $\gamma = 0.2$, and the time step $\Delta t = 0.02$. The number of the Chebyshev polynomials $M = 2000$. In order to assist transitions between the local minima of the local order parameter $\Phi_i$, a small noise of the order of $\delta\Phi_i = 0.1$ is added to the slave-boson amplitudes obtained from the previous time-step at the beginning of the new Gutzwiller iteration.

Thermal Mott transition {#sec:thermal-mott}
=======================

![(Color online) \[fig:phase-diagram\] (a) Schematic phase diagram of the half-filled Peierls-Hubbard model on the triangular lattice. The phase boundary is estimated from the Gutzwiller MD simulations. The nearest-neighbor hopping $t^{(0)}_{1} = 1$ serves as the units for energy and temperature. The critical line $U_c(T)$ ends at $T = 0$ and $T_c \approx 0.9$. (b) The spatially averaged double-occupancy $\langle D \rangle = \sum_i^N D_i/N$ versus the temperature with the Hubbard repulsion set at $U = 13.1$. The shaded area marked the coexistence region of the phase transition. ](phase-diagram.pdf){width="0.99\columnwidth"}

A schematic phase diagram of the half-filled Peierls-Hubbard model is shown in Fig. \[fig:phase-diagram\](a) based on the Gutzwiller MD simulations. Similar phase diagrams are obtained from the finite-temperature Gutzwiller calculation [@wang10; @mezio17] as well as DMFT [@georges96; @bulla01; @liebsch09]. The first-order MIT line $U_c(T)$ ends at two critical points at $T = 0$ and a finite $T_c > 0$. The first analytical theory of the zero-temperature MIT was proposed by Brinkman and Rice using the Gutzwiller method [@brinkman70]. In the Brinkman-Rice scenario, the double-occupancy decreases linearly with increasing $U$ according to: $D(U) = \frac{1}{4}(1 - U/U_c)$, where the critical $U_c = 8 |E_{\rm band}|$ and $E_{\rm band}$ is the electron band energy. It is worth noting that the quasi-particle behavior in the Mott transition is well captured, at least qualitatively, by the Brinkman-Rice theory, although it cannot describe the incoherent electronic excitations and the Hubbard bands.

At finite temperatures, the Mott transition becomes first-order with the double-occupancy showing a discontinuity $\Delta D$ at the transition point. This is demonstrated in Fig. \[fig:phase-diagram\](b) which shows the temperature dependence of the spatially averaged double-occupancy $\langle D \rangle$ serving as an order parameter of the Mott transition. In addition to a clear discontinuity in $\langle D \rangle$, the first-order nature of the MIT also manifests itself in the hysteresis loop from the heating and cooling curves.

Here we consider the evolution of the Hubbard model subject to a temperature quench. For simplicity, we assume the temperature of the system is controlled by a reservoir which only interacts with the lattice through the stochastic forces $\bm\xi_i(t)$ in the Langevin equation (\[eq:langevin\]). The system is initially equilibrated at a temperature $T_i$ in the metallic phase. At $t = 0$, the reservoir temperature is suddenly increased to a higher temperature $T_f$ across the phase boundary. Interestingly, due to the negative slope of the phase boundary $dU_c/dT < 0$ [@note-dUdT], such a thermal quench drives the system from the correlated metal to a Mott insulator. Within the BO approximation, the electrons are assumed to quickly reach equilibrium at the instantaneous lattice temperature, i.e. $k_B T_e(t) = m \langle \dot{\mathbf u}(t)^2 \rangle /3$. The electron temperature thus increases continuously from $T_i$ to $T_f$ as the lattice is heated up by the reservoir.

![(Color online) \[fig:trace\] (a) The various energy densities versus time in a MIT driven by a temperature quench. Here $E_e$, $E_l$, and $E_k$ denote the electronic, lattice (elastic), and atomic kinetic energies, respectively. The initial and final temperatures are $T_i = 0.09$ and $T_f = 0.13$, respectively. The (red) horizontal dashed line corresponds to $k_B T_f$. (b) The averaged double occupation $\langle D \rangle$ and the volume fraction of the Mott phase $\phi_M$ versus time. The shaded area bounded between $t_1$ and $t_2$ indicates the region where most of the phase transformation takes place. ](MIT-trace.pdf){width="0.99\columnwidth"}

The time dependence of the various energy densities is shown in Fig. \[fig:trace\](a) for a quench from $T_i = 0.09$ to $T_f = 0.13$. The electronic energy $E_{\rm elec}$ includes both the hopping and Hubbard term. The lattice energy $E_l$ comes from the elastic $K$ term in the Hamiltonian. The atomic kinetic energy density is $E_k = m \sum_i |\dot{\mathbf u}|^2 / 2N$. The phase transformation mostly takes place in the temporal regime $t_1 < t < t_2$, marked by the shaded area. Immediately after the quench, the lattice temperature gradually increases toward the final temperature $T_f$. The enhanced thermal fluctuation of displacements in this initial period ($t < t_1$) result in an increase in elastic energy $E_l$. The enlarged fluctuations of $\langle |\mathbf u_i|^2 \rangle$ also lead to an enhanced variance of the nearest-neighbor hopping $t_{ij}$, which in turn gives rise to a slightly reduced electron energy. However, this relatively short transient phase does not play a significant role in terms of phase transformation.

![(Color online) \[fig:snapshots\] Snapshots of the on-site double occupancy $D(\mathbf r_i) = \langle n_{i, \uparrow} n_{i, \downarrow} \rangle$ at four different Mott volume fractions, corresponding to $t = 30$, 50, 58, and 65, during the thermal driven MIT. See Fig. \[fig:trace\](b) for the corresponding time dependence of the Mott volume fraction and the average double-occupancy. ](snapshots.pdf){width="0.99\columnwidth"}

As discussed in Sec. \[sec:intro\], the double occupancy is one of the order-parameters for characterizing the Mott transitions. Another important quantity is the fraction of Mott phase $\phi_M$, which is often used in the study of phase transformation kinetics. The time dependence of these two quantities during the thermal Mott transition is shown in Fig. \[fig:trace\](b); here $\langle D \rangle = \sum_i D_i / N$ is the spatial average of the local double occupation. Numerically, we label a given site as the Mott insulator if its double-occupancy $D_i < 10^{-8}$. It is worth noting that the transformed volume $\phi_M$ is negligible in the initial period when $t < t_1$. In fact, the average $\langle D \rangle$ even slightly increases due to the lowered electron energy.

After this initial period, the system is trapped in a meta-stable metallic state with a temperature $T \approx T_f$. The transformation to the Mott insulator mainly takes place in the temporal window $[t_1, t_2]$, in which the Mott volume fraction $\phi_M$ gradually increases from zero to its equilibrium value at $T_f$ while the average $\langle D \rangle$ slowly decays to zero, see Fig. \[fig:trace\](b). Once the system reaches the new equilibrium at $T_f$ for both the lattice and electrons, the electron binding energy is renormalized to almost zero. This also indicates a vanishing electronic force $\mathbf f_i^e \approx 0$; see Eq. (\[eq:f\_e\]). Consequently, the lattice system becomes essentially a collection of free harmonic oscillator, and both energy densities $E_l$ and $E_k$ approach $k_B T_f$ consistent with the equipartition theorem.

Nucleation and growth {#sec:nucleation}
=====================

To further investigate the nature of the phase transformation, Fig. \[fig:snapshots\] shows the spatial profile of the local double-occupancy $D(\mathbf r_i) = \langle n_{i, \uparrow} n_{i, \downarrow} \rangle$ at four different times during the relaxation. These snapshots clearly show that the transition to the new phase is initiated by the nucleation of Mott insulating puddles, i.e. regions with vanishing double-occupancy, in the metallic matrix. Importantly, the double-occupation in the metallic regions remains roughly the same. This observation also means that the gradual decrease of the average $\langle D \rangle$ in Fig. \[fig:trace\](b) results mostly from the increasing volume fraction of the Mott phase in a heterogeneous electronic state, instead of the decrease of a a homogeneously distributed double-occupancy.

This scenario is further confirmed by the histogram or probability distribution function $h(D)$ of the double-occupancy, shown in FIg. \[fig:hist\_db\](a), at various Mott volume fractions. The distribution function consists mainly of a delta-function at $D = 0$ corresponding to the Mott clusters, and a broad hump centered at $D \approx 0.055$, which remains roughly the same throughout the transition. This result is consistent with picture of two coexisting phases in the intermediate state. As the phase transformation progresses, the delta-peak keeps increasing at the expense of the finite-width distribution at finite $D$. Interestingly, the distribution of the local electron density rather uniform, as demonstrated by its histogram $h(n)$ shown in Fig. \[fig:hist\_db\](b). The on-site electron density $n_i$ is very narrowly distributed around the half-filling with a width $\delta n \lesssim 0.01$, and is dominated by a huge delta-peak at $n = 0.5$, with contributions from both the metallic and Mott regions.

![(Color online) \[fig:hist\_db\] Computed histogram of (a) double-occupancy and (b) on-site electron density at varying Mott volume fraction for the temperature-quench simulations of Fig. \[fig:trace\]. ](hist_db.pdf){width="0.99\columnwidth"}

As the phase transformation progresses, both the number and size of the Mott clusters increase with time. At the late stage of the phase transformation, the Mott clusters start to merge through a process akin to the percolation transition. Of central interest here is the evolution of the cluster-size distribution $n(s, t)$, which is shown in Fig. \[fig:histogram\] for different stages of the phase transformation. For small clusters, it is expected that their distribution is mainly determined by the interface energy which can be expressed as $\varepsilon(A) = \sigma A$, where $\sigma$ is the surface tension and $A \sim s^{(d-1)/d}$ is the surface area of the nuclei with size $s$. In two dimensions their distribution, given by the Boltzmann probability, is then $$\begin{aligned}
    \label{eq:ns_exp}
    n(s, t) \sim \exp(- \beta \sigma_{\rm eff} \sqrt{s}),  \end{aligned}$$ where $\sigma_{\rm eff}$ is an effective surface tension. As shown in Fig. \[fig:histogram\], the expression in Eq. (\[eq:ns\_exp\]) indeed fits the numerical distributions very well for small $s$, particularly at the initial stage of the phase transformation. The effective surface tension obtained from the fitting is $\sigma_{\rm eff} \approx 0.1$.

![(Color online) \[fig:histogram\] The probability distribution function $n(s)$ of the size of Mott clusters at different stages, characterized by the volume fraction of the Mott clusters, of the phase transformation. The curves are shifted vertically for clarity. The corresponding curve of volume fraction $\phi_M$ versus time is shown in Fig. \[fig:trace\](b). Each curve is obtained by average over a small time interval $[t-\Delta t, t+\Delta t]$ with $\Delta t = 2.5$, and over 15 independent simulations. The solid lines are the fit using Eq. (\[eq:ns\_exp\]). The dashed line corresponds to power law $s^{-\tau'}$ with $\tau' = 1.42$. The estimated error of the exponent is $\pm 0.06$. ](histogram.pdf){width="0.95\columnwidth"}

Significant deviation from Eq. (\[eq:ns\_exp\]) starts to occur at cluster size $s_c \approx 100$, which provides an estimate for the so-called critical size. The concept of critical nuclei was first discussed in the classical droplet model. In this picture, the decay of a meta-stable state is initiated by a collection of independent fluctuating droplets. Their quasi-equilibrium distribution is again given by the Boltzmann probability $n_0(s) \propto \exp[- \beta \varepsilon(s)]$, where $\varepsilon(s)$ is the minimum work required to form a cluster of size $s$. In general $d$ dimensions, this minimum work is $$\begin{aligned}
    \varepsilon(s) = \sigma s^{(d-1)/d} + \Delta f\, s.\end{aligned}$$ The first term represents the interface energy as discussed above, while the second term denotes the volume energy of the droplets with $\Delta f$ being the free-energy difference between the two competing phases. For $\Delta f > 0$, the droplets correspond to thermodynamically unfavorable phase, and $n_0(s)$ describes the fluctuations of these small domains. On the other hand, for $\Delta f < 0$, as in a quench into meta-stable states, $\varepsilon(s)$ describes the competition between the surface tension and the volume energy gain. In this case, the distribution $n_0(s)$ decreases as a function of $s$ for small clusters until it reaches the minimum at the critical size $s_c \sim (\sigma/ |\Delta f|)^d$. For $s > s_c$, the probability $n_0(s)$ then grows exponentially with increasing cluster size. These super-critical clusters are thus seeds of the nucleation process. The maximum of the work function $\varepsilon_{c} = \varepsilon(s_c) \sim \sigma^d/|\Delta f|^{d-1}$ represents an energy barrier to nucleate a critical cluster.

![(Color online) \[fig:Nc\] The total number of nuclei $\mathcal{N}_c$ as a function of time for various $\Delta T = T_f - T_i$ from the Gutzwiller MD simulations. The initial temperature $T_i = 0.09$. The dashed lines are fit to the linear equation in Eq. (\[eq:linear\]). ](num_vs_t.pdf){width="0.99\columnwidth"}

Once the phase transformation is initiated by these seed clusters, the quasi-equilibrium $n_0(s)$ no longer describes the time-dependent cluster distribution. Numerical methods are often required to compute the evolution of the $n(s, t)$. As Fig. \[fig:histogram\] shows, the distribution function extends beyond the estimated critical size with an increasingly longer tail. Interestingly, as the volume fraction approaches $\phi_M^* \approx 0.42$, the cluster-size distribution exhibits a power-law tail $$\begin{aligned}
    \label{eq:tau}
    n(s, t^*) \sim 1/s^{\tau'}, \end{aligned}$$ with an exponent $\tau' \approx 1.42$. In percolation transitions, the size of the percolating clusters follows a similar power-law distribution, $n(s; p \to p_c) \sim s^{-\tau}$, when the percolation probability $p$ approaches the critical value $p_c$. Here $\tau$ is called the Fisher’s exponent, and its numerical value is $\tau = 187/91 \approx 2.055$ for 2D percolation [@stauffer94]. Despite the similarity, it is worth noting that the nature of these two power-law behaviors is very different. The distribution function $n(s, p)$ in the case of percolation is stationary and describes the static structure of the percolating clusters. The emergence of the power-law characterized by Fisher’s exponent $\tau$ also indicates a divergent length scale $\xi \sim (p - p_c)^{-\nu}$ related to the underlying critical point [@stauffer94]. On the other hand, the power-law in Eq. (\[eq:tau\]) describes a [*transient*]{} distribution of the Mott clusters during a first-order phase transformation. It does not necessarily imply a divergent length scale, also see discussion in Sec. \[sec:structure\]. On the other hand, the integrated distribution $\mathcal{D}(s) = \int_0^{t^*} n(s, t) dt \sim 1/s^{\tau}$ also exhibits a power-law tail with $\tau \approx 2.03$, although it is less pronounced due to the limitation of system sizes. This power-law tail might be related to the Ising critical end point of the Mott transition, where the distribution of the geometrical spin clusters is described by a similar power-law distribution with $\tau =  379/187 \approx 2.027$ [@janke05]. A more detailed study of the phase coexistence regime is required to clarify this connection.

We next consider the time dependence of the total number of nuclei, which is related to the distribution function via the integral $\mathcal{N}_c(t) = \int_0^\infty n(s, t) \, ds$. Fig. \[fig:Nc\] shows the numerically computed $\mathcal{N}_c$ versus time for varying depths of the temperature quench $\Delta T = T_f - T_i$. Here we only show the early stages of the nucleation, in which two distinct kinetic behaviors can be seen from the figure: (i) the incubation period with negligible number of clusters, and (ii) the constant nucleation regime in which the total number of nuclei increases linearly with time. To understand these behaviors, we briefly discuss the classical nucleation theory first developed by Becker and Döring [@becker35], Zedolvich [@zeldovich43], and Frenkel [@frenkel46]. The distribution function $n(s, t)$ is governed by the Fokker-Planck equation [@gunton83] $$\begin{aligned}
    \label{eq:FK}
    \frac{\partial n}{\partial t }= -\frac{\partial J}{\partial s} =  \frac{\partial }{\partial s}\left[ R \, n_0\, \frac{\partial }{\partial s} \left(\frac{n}{n_0} \right) \right],\end{aligned}$$ where $J$ is the flux of nuclei in the cluster-size space, with its definition given in the second equality, and $R = R(s)$ is the diffusion coefficient. This theory essentially describes the random work of the nuclei in the cluster-size space. The flux can also be expressed as $J = -R \,\partial n / \partial s + A n$, where $A$ plays the role of the drift velocity. The coefficient $A$ is related to $R$ by the Einstein relation which originates from the condition that $J = 0$ for the equilibrium distribution $n(s, t) = n_0(s)$, leading to the flux expression in Eq. (\[eq:FK\]).

The regime with a constant nucleation rate, on the other hand, corresponds to a nonequilibrium steady-state solution with $J(s, t) = \mathcal{I} = $ constant, subject to the so-called “sink and source" boundary conditions: $n(s \to \infty, t) = 0$ and $n(s, t) \to n_0(s)$ for $s \ll s_c$ [@becker35]. In this case, the total number of clusters satisfies the simple differential equation $d\mathcal{N}_c/dt = -\int_0^\infty (\partial J/\partial s) \, ds = \mathcal{I}$, which can be readily integrated to yield $$\begin{aligned}
    \label{eq:linear}
    \mathcal{N}_c(t) =  \mathcal{I} (t - t^*),\end{aligned}$$ where we have introduced the so-called induction time $t^*$ as the constant of integration. Fig. \[fig:DT\_dep\] shows the numerically extracted nucleation rate $\mathcal{I}$ and induction time $t^*$ versus $\Delta T$. As expected, increasing the depth of the quench enhances the nucleation rate, while shortening the incubation time for nucleation. The constant flux $\mathcal{I}$ can be computed from Eq. (\[eq:FK\]), $$\begin{aligned}
    \mathcal{I} = \left( \int_0^\infty \frac{ds}{R\,n_0} \right)^{-1} \approx \mathcal{I}_0 \, \exp \! \left(\frac{-\varepsilon_c}{k_B T_f} \right).\end{aligned}$$ Here the saddle-point approximation has been used to compute the integral, $\varepsilon_c = \sigma^2/2 |\Delta f|$ is the nucleation barrier discussed previously, and the prefactor $\mathcal{I}_0$ is related to the curvature at the work function maximum [@gunton83]. In the inset of Fig. \[fig:DT\_dep\](a), we plot the negative logarithm of the nucleation rate as a function of $\Delta T$. Since $\varepsilon_c = - k_B T_f \ln (\mathcal{I}/\mathcal{I}_0)$, our result shows that the energy barrier decreases roughly linearly with the quench depth.

![(Color online) \[fig:DT\_dep\] (a) The nucleation rate $\mathcal{I}$ and (b) induction time $t^*$ as functions of the quench depth $\Delta T$ extracted from the curves in Fig. \[fig:Nc\]. The inset in panel (a) shows the negative logarithm of $\mathcal{I}$ versus $\Delta T$. The lines are guide for eye. ](DT_dep.pdf){width="0.99\columnwidth"}

Interface-controlled mechanism {#sec:interface}
==============================

After the initial nucleation period, the Mott clusters start to occupy contiguous areas and their growth is hindered. A quantity of central interest here is the volume fraction $\phi$ of the transformed phase. In our case, we are interested in the fraction $\phi_M$ of the Mott phase, which is already introduced in Fig. \[fig:trace\](b) and is shown in Fig. \[fig:phi\_M\](a) as a function of time for varying depths of the temperature quench. The standard model for the time dependence of the transformed volume fraction was developed by Kolmogorov, Johnson, Mehl, and Avrami (KJMA model) in the thirties [@kolmogorov37; @avrami39; @avrami40; @johnson39]. The KJMA model is a mean-field type theory for treating the impingement of clusters. To account for this effect, the KJMA introduced the concept of extended volume $V^e$, which is the volume occupied by the new phase [*if*]{} they can grow without being hindered by the presence of other nuclei. Let $dV_e = d\phi_e \mathcal{V}$ be the increase of the extended volume within a small time interval $dt$, and $\mathcal{V}$ is the fixed total volume of the system. The real volume $dV$ that is available for the growing nuclei is the fraction $ f = 1 - \phi$ of the extended volume, i.e. $dV/dV^e = d\phi/d\phi^e = f$. Integrating this equation leads to a relation between the fraction of the actual transformed volume and that of the extended one $$\begin{aligned}
    \label{eq:avrami}
    \phi(t) = 1 - \exp[-\phi^e(t)] = 1 - \exp\left( - K t^n \right).\end{aligned}$$ In the second expression, we have assumed a power-law growth of the extended volume and $K$ is a pre-factor proportional to the nucleation rate. This simple expression, also called the Avrami equation, has been widely used in the analysis of phase transformation experiments. It has been suggested that $n$ can be partitioned as $n = \alpha + d p$ [@christian75], where $a$ depends on the nucleation rate ($\alpha = 0$ for existing nuclei, and $\alpha = 1$ for constant nucleation), $d$ is the spatial dimensions, and $p$ is related to the growth mechanism. Notably, $p = 1$ for interface-controlled growth, and $p = 1/2$ for diffusion controlled process.

![(Color online) \[fig:phi\_M\] (a) The area fraction $\phi_M$ of the Mott phase versus time for varying quench depth $\Delta T = T_f - T_i$, here the initial temperature $T_i = 0.09$. (b) The log-log plot of $-\ln(1-\phi_M)$ versus the shifted time $t - t^*$. The dashed lines correspond to the Avrami equation with an exponent $n = 3$. ](phi_vs_t.pdf){width="0.99\columnwidth"}

From the Gutzwiller MD simulations, we can extract the exponent $n$ by plotting the effective extended fraction $ \phi^e_M = -\ln(1 - \phi_M)$ versus the simulation time, as shown in Fig. \[fig:phi\_M\](b). Since the Avrami equation is obtained under the assumption of constant cluster growth, the time here is shifted by the induction time $t^*$ obtained in Fig. \[fig:DT\_dep\](b). The log-log plot clearly shows a linear regime with slope $n = 3$ for small $t$. This result is consistent with the standard interface-controlled ($p = 1$) homogeneous nucleation ($a=1$), giving rise to $n = d + 1 = 3$ in our case. In general, the evolution of the extended volume fraction can be expressed as $\phi^e(t) \propto \int_0^t  \langle {s}(t; t') \rangle\, \mathcal{I}(t') \, dt'$, where $\mathcal{I}(t)$ is a time-dependent nucleation rate and $\langle {s}(t; t') \rangle$ is the size of nucleus at time $t$ that was nucleated at $t'$. For constant nucleation rate, as demonstrated in Sec. \[sec:nucleation\] for the early stage of the transformation, and assuming time-invariant cluster growth, we have $\phi^e(t) \propto \mathcal{I} \int_0^t \langle {s}(t - t') \rangle dt'$.

In the case of interface-controlled mechanism, the growth rate of a cluster of size $s$ is proportional to its surface “area" $A \sim s^{(d-1)/d}$, i.e. $ds/dt \sim s^{(d-1)/d}$, from which one obtains a growth law $$\begin{aligned}
    \langle {s} \rangle \sim t^d.\end{aligned}$$ The characteristic linear size of the cluster $\ell \sim \langle s \rangle^{1/d}$ thus increases [*linearly*]{} with time in interface-controlled mechanism. The corresponding extended volume fraction then obeys a power-law behavior $\phi^e(t) \sim t^{n}$ with the Avrami exponent $n = d+1$. Our results from Fig. \[fig:phi\_M\](b) indicates that the early stage of the phase transformation is consistent with the interface-controlled cluster growth, which is also directly verified in Fig. \[fig:size\] from our Gutzwiller MD simulations. By taking into account the incubation period, the time dependence of the average cluster size is well approximated by a parabolic function $\overline{s} = a t^2 + b$. In particular, a linear segment with unity slope can be seen in the log-log plot of the characteristic linear size, defined as $\ell = \langle {s} \rangle^{1/2}$, versus the shifted time $t - t^*$, where $t^*$ is the induction time from Fig. \[fig:DT\_dep\](b).

![(Color online) \[fig:size\] The average size of nuclei $\overline{s}$ as a function of time for varying depths of the temperature quench. The dashed lines are the fit with parabolic functions. The inset shows the log-log plot of average linear size of clusters ${\ell} = \langle {s} \rangle^{1/2}$ versus the shifted time $t - t^*$. The dashed lines correspond to linear function ${\ell} \sim (t - t^*)$. ](size_vs_t.pdf){width="0.98\columnwidth"}

The interface-controlled growth mechanism here indicates that new Mott insulating sites tend to be induced at the perimeter of the already formed Mott puddles. Physically, this is because the weakened electron bonding force in the vicinity of the Mott region allows the neighboring atoms to overcome the energy barrier more easily and tunnel into the insulating state through thermal fluctuations. One can also see this mechanism from the local slave-boson Hamiltonian Eq. (\[eq:H\_SB\]). The kinetic energy, represented by the $\mathcal{M}_\alpha \, \Delta_{i, \alpha}$ term, is reduced in the neighborhood of a Mott cluster due to the reduced $\Delta_{i, \alpha}$ parameter; see Eq. (\[eq:Delta\]). The competition with the Hubbard $\mathcal{U}$ thus results in a slave boson state with vanishing $d_i$. This local competition underlies the microscopic cluster growth process in the Gutzwiller MD simulations.

Avalanche behavior {#sec:avalanche}
==================

![(Color online) \[fig:scaled\_phi\] The Mott volume fraction $\tilde \phi_M$ normalized to its maximum versus the scaled time $(t - t_0)/w$ for varying $\Delta T$. The black dashed line correspond to the Avrami Eq. (\[eq:avrami\]) with $n = 3$. The red dash-dotted line is the Avrami equation with $n = 3.6$ shifted in time to fit the large time part of the transformation. ](scaled_phi.pdf){width="0.96\columnwidth"}

It is worth noting that the phase transformation of the 2D kinetic Ising model is well described by the standard Avrami equation with $n = 3$ [@shneidman99; @shneidman99b]. Moreover, cell dynamics simulations based on the stochastic 2D TDGL Eq. (\[eq:TDGL\]) with a scalar order parameter gives results quantitatively consistent with the $n = 3$ KJMA model [@iwamatsu08]. In our case, however, the standard KJMA kinetics only lasts a relatively small fraction of the total phase transformation time. After this initial stage of constant nucleation, the transformation speeds up significantly compared with the prediction of the Avrami equation. To further quantify the Mott transition kinetics, Fig. \[fig:scaled\_phi\] shows the normalized Mott volume fraction versus the scaled time $(t - t_0)/w$. Interestingly, by properly choosing the origin $t_0$ and the duration $w$ of the transformation, the numerical data from different $\Delta T$, shown in Fig. \[fig:phi\_M\](a), roughly collapse on the same curve, indicating a new universal kinetic behavior for the Mott transition. With the origin of the phase transformation fixed at $t = t_0$, there is two parameters $K$ and $n$ left for fitting the numerical data with Eq. (\[eq:avrami\]). As shown in Fig. \[fig:scaled\_phi\], the best-fitting result with the exponent set at $n = 3$ only captures the initial part of the transformation curve up to $\tilde{\phi}_M \sim 0.15$. Beyond that, one can see a clear acceleration of the phase transformation compared with the prediction of the KJMA model. In fact, by shifting the origin to a later time, Eq. (\[eq:avrami\]) with a larger exponent $n = 3.6$ (the red dash-dotted curve) seems to well describe the later stage of the phase transformation.

In addition to the increasing rate, this new stage of the transformation is characterized by kinetics similar to the Barkhausen effect [@durin05; @sethna01] in the magnetization or demagnetization processes of ferromagnets. It should be noted that the relatively smooth curves shown in Fig. \[fig:scaled\_phi\] were obtained by averaging over at least ten independent runs. Individual transformation curves as that shown in Fig. \[fig:jumps\] (right axis), however, exhibit step-like features in the accelerated regime of the phase transformation. These steps in $\phi_M$ are also marked by the delta-like peaks in the time dependence of the instantaneous transformation rate $d\phi_M/dt$, see Fig. \[fig:jumps\] (left axis). Microscopically, the Barkhausen effect is caused by sudden changes in the size and orientation of ferromagnetic domains [@urbach95; @cizeau97; @narayan96]. Similarly, the jumps in $\phi_M$ here also correspond to the local Mott transition of large domains of the lattice sizes, which is confirmed from the real-space snapshots of the double-occupancy in the vicinity of these steps. Moreover, we find that the formation of these large Mott clusters are not random, but tend to occur in the vicinity of the already existing Mott puddles.

![(Color online) The transformation rate $d\phi_M/dt$ as a function of time for a particular temperature-quench simulation from $T_i = 0.09$ to $T_f = 0.18$. The smooth black line indicates the averaged transformation rate. Also shown in the right axis is the time dependence of the Mott volume fraction $\phi_M$. The four distinct stages of the phase transformation are also outlined above. \[fig:jumps\] ](jumps.pdf){width="0.99\columnwidth"}

Quenched disorder, which introduces random pinning centers for the propagating domain walls, plays a crucial role in the Barkhausen effect [@sethna01; @urbach95; @cizeau97; @narayan96]. The disorder induces multiple local minima in the configuration space of domain-walls, and the jumps in magnetization correspond to transitions between nearly local minima. Thermal fluctuations, on the other hand, are not necessarily an important factor in the domain-wall dynamics. Remarkably, despite the absence of quenched disorder in our simulations, the transformation to the Mott phase still exhibits the avalanche-like behavior as shown in Fig. \[fig:jumps\]. The intermittent bursts of local Mott clusters are likely due to a self-generated disorder in the intermediate mixed-phase states.

This scenario can be qualitatively understood from the viewpoint of the energy landscape of the Mott order parameters. At the early stage of the phase transformation, the system moves toward a well defined global minimum in a relatively smooth valley through the nucleation and growth of the Mott clusters. As the transformed volume fraction $\phi_M$ increases, more and more lattice sites are disconnected from the rest due to the vanishing effective hopping $t^{\rm eff}_{ij} = \mathcal{R}_i \mathcal{R}_j t_{ij} \to 0$. These diluted bonds in turn introduce off-diagonal disorder to the electron Hamiltonian Eq. (\[eq:H\_renorm\]) and modify the localization of the electron wavefunctions. After integrating out the electrons, an increasingly rugged energy landscape for the slave-bosons is created by the enhanced disorder. Similar to the depinning transition of driven domain-walls, the jumps in the Mott volume fraction $\phi_M$ thus corresponds to the transitions between local minima of the resultant energy landscape.

![(Color online) \[fig:size\_phi\] The average cluster size $\langle {s} \rangle$ versus the Mott volume fraction $\phi_M$ for varying depth of the temperature quench. The inset shows the log-log plot of $\langle {s} \rangle$ versus time. ](size_phi.pdf){width="0.97\columnwidth"}

It is worth noting that numerically the avalanche behavior in our MD simulations comes mostly from the Gutzwiller iteration even within a single time-step. The Mott transition of the whole cluster is essentially a collective electronic response. This process is also similar to the avalanche behaviors in, e.g. the random field Ising model (RFIM) [@sethna93; @sethna05] or the depinning transition of domain walls [@urbach95; @cizeau97; @narayan96]. The evolution of these spin systems is often governed by adiabatic or athermal dynamics. The avalanche is initiated by a single spin-flip which triggers a chain-reaction-like response of neighboring spins in one time-step.

The effective energy functional $\mathcal{E}[\Phi_i]$ of the Mott order parameter can be formally obtained from Eq. (\[eq:free-energy\]) by integrating out the electron and the lattice degrees of freedom. Importantly, this process introduces long-range electron-mediated “interactions" between the local slave-boson variables $\Phi_i$, which could potentially create a complex landscape. This nonlocal interaction can be demonstrated by the effective spin model discussed in Appendix \[sec:spin-model\]. Assuming particle-hole symmetry, which is a good approximation in our case as the local electron density is very close to half-filling; see Fig. \[fig:hist\_db\], one can introduce classical spins $\mathbf S_i$ such that $S_i^x = d_i p_i + p_i d_i$ and $S_i^z = p_i^2 - d_i^2$. Interestingly, the effective Hamiltonian for these spin variables is the transverse-field Ising model with a nearest-neighbor exchange interaction $\mathcal{J}_{ij} = 8 t_{ij} (\rho_{ij} + \rho_{ji})$, where $\rho_{ij}$ is the reduced density matrix of the quasi-particles and depends nonlocally on other spins.

The nonlocal nature of the effective interaction also manifests itself in the peculiar growth dynamics of the Mott clusters. We found that the increase of the average cluster size $\langle s \rangle$ is governed neither by a power-law \[see inset of Fig. \[fig:size\_phi\]\] or exponential function. Instead, it increases exponentially with the Mott volume fraction $$\begin{aligned}
    \label{eq:exp_phi}
    \langle s \rangle \sim \exp(\gamma \phi_M),\end{aligned}$$ where $\gamma \approx 6.8$ – 6.9 is a constant weakly dependent on $\Delta T$. Since the fraction $\phi_M$ describes how many sites have been transformed to the Mott phase in the whole lattice, this dependence indicates that the growth of individual clusters in the avalanche regime depends on the [*global*]{} property of the system.

Before concluding this section, we outline in Fig. \[fig:jumps\] the four distinct kinetic behaviors observed in the Mott transition: (i) the initial incubation period, (ii) the constant nucleation regime, (iii) the accelerated transformation with avalanche behavior, and (iv) the coagulation at the late stage. In the incubation period, there are no Mott clusters of appreciable sizes and nearly negligible transformation rate. Once the system enters the quasi-steady state with a constant nucleation rate, the Mott fraction steadily increases with a $\phi_M \sim t^3$ dependence, consistent with the interface-controlled cluster-growth mechanism. It is also worth noting that the transition rate $d\phi_M/dt$ exhibits both positive and negative fluctuations, indicating that some Mott clusters would turn back to become metallic throughout the transformation. As the Mott volume fraction reaches $\phi_M \approx 0.15$, the phase transformation kinetics is characterized by avalanche behavior and a non-local cluster-growth phenomenon. And finally, the transition is slowed down at the late stage as the growth of the Mott cluster is blocked by the coagulation effect.

Structure of the mixed-phase states {#sec:structure}
===================================

![(Color online) \[fig:Sq\] The logarithm of the structure factor $S(\mathbf q)$ at (a), (b) $\phi_M = 0.008$ and (c), (d) $\phi_M = 0.4$ of the temperature quench shown in Fig. \[fig:trace\]. They correspond to simulation time $t = 15$ $t = 52$, respectively. ](sq.pdf){width="0.99\columnwidth"}

We next characterize the structure of the intermediate mixed-phase states during the phase transformation. In order to compare with results from the kinetic Ising model, we first map the electronic state to an effective Ising configuration by introducing a binary variable $\sigma_i$ such that $\sigma_i = +1$ if the corresponding site-$i$ is metallic, and $\sigma_i = -1$ if it is a Mott insulator. Again, numerically we define a site to be a Mott insulator if its double-occupancy $D_i < 10^{-8}$. Our central interest here is the time-dependent two-point correlation function $$\begin{aligned}
    G(\mathbf r_{ij}, t) = \langle \delta\sigma_i(t) \delta\sigma_j(t) \rangle = \langle \sigma_i(t) \sigma_j(t) \rangle - \langle \sigma(t) \rangle^2,\end{aligned}$$ where $\delta\sigma_i \equiv \sigma_i - \langle \sigma \rangle$ denotes the fluctuations relative to the “ferromagnetic” order. In this mapping, the metallic (Mott) phase corresponds to the ferromagnetic state $\langle \sigma \rangle > 0$ ($\langle \sigma \rangle < 0$). The Fourier transform of the correlation function gives the structure factor $$\begin{aligned}
    S(\mathbf q, t) = \sum_{\mathbf r} G(\mathbf r, t) e^{i \mathbf q \cdot \mathbf r} = \Bigl \langle \Bigl | \sum_i \delta \sigma_i(t) \,e^{i \mathbf q\cdot \mathbf r_i} \Bigr|^2 \Bigr\rangle,\end{aligned}$$ During the phase transformation, the presence of a finite ferromagnetic order introduces a huge peak at $\mathbf q = 0$. In order to see details of the structure factor, we plot in Fig. \[fig:Sq\] the logarithm of $S(\mathbf q)$ computed at two different simulation times of the temperature quench shown in Fig. \[fig:trace\]. An emergent circular symmetry of the structure factor indicates an isotropic phase-separated states. The central peak becomes sharper and stronger as the phase transformation progresses.

![(Color online) \[fig:sq\_cmp\] Log-log plot of the circularly averaged structure factor $S(\mathbf q)$ versus the $|\mathbf q|$ at four different simulation times. The black dashed lines are guide for the line. The $1/q^3$ power-law dependence is shown by the red long dashed line. ](sq_cmp.pdf){width="0.9\columnwidth"}

Fig. \[fig:sq\_cmp\] shows the circularly averaged structure factor as a function of the modulus of $\mathbf q$ at four different Mott volume fractions. Although not much information can be obtained from the small $\mathbf q$ behaviors, the large wavevector dependence is consistent with the Porod’s law which states that the structure factor for a two-phase system with interfaces of negligible thickness should behave as $S(q) \sim 1/q^{d+1}$ [@porod82], where $d$ is the spatial dimension. The deviation from the Porod’s law at the Brillouin zone boundary results from thermal fluctuations and finite width of the metal-insulator interfaces. The power-law $q$-dependence also indicates a short-distance singularity of the correlation function $G(\mathbf r) \approx (1 - \langle \sigma \rangle^2)\,(1 - A | \mathbf r|)$ at small $|\mathbf r|$, where $A$ is a constant. This linear-$r$ dependence at short distances is also confirmed in our simulations; see Fig. \[fig:gr\](a) for $G(r)$ computed at various Mott volume fractions $\phi_M < 0.5$.

The correlation length $\xi$ estimated from these correlation curves is found to increase monotonically with time before reaching the mid-point $\phi_M = 0.5$ of the phase transformation. A more systematic determination of $\xi$ will be discussed below. Importantly, despite the emergence of a transient power-law size distribution of the Mott clusters discussed in Sec. \[sec:nucleation\], the correlation length remains finite throughout the phase transformation, as evidenced by the log-log plot of the same correlation curves shown in Fig. \[fig:gr\](b). This result thus means that the probability of two separate sites being on the same large cluster is exceedingly small. A similar example is the 3D Ising model, in which a power-law distribution of geometrical spin clusters occurs at a temperature which is different from the critical point.

On the other hand, the correlation curves at different times seem to exhibit a self-similarity. To investigate potential universal behaviors, we compute the normalized correlation functions $\tilde G(r, t) \equiv G(r, t)/G(0, t)$ and plot them against a rescaled distance $r \to r/\xi$ in Fig. \[fig:gr\_scaled\]. The rather nice data points collapsing in the plot strongly suggests a universal correlation function for the mixed-phase state of the Mott transition. The scaling factor $\xi$, determined from the condition $\tilde G(\xi) = 1/e$, also provides a measure of the finite correlation length. The inset of Fig. \[fig:gr\_scaled\] shows the numerical $\xi$ versus the Mott volume fraction. Interestingly, the correlation length $\xi$ is at most of the order of 10 lattice constants.

![(Color online) \[fig:gr\] (a) Correlation function $g(r)$ at varying Mott volume fraction $\phi_M$ from the temperature quench simulation shown in Fig. \[fig:trace\]. (b) The same curves shown in log-log plot. ](gr.pdf){width="0.99\columnwidth"}

Fig. \[fig:gr\_scaled\] also shows the comparison of the numerical universal function with the analytical formula obtained by Sekimoto based on the interface-controlled mechanism within the KJMA theory [@sekimoto84; @sekimoto86]. In his derivation, for simplicity, the clusters of the stable phase are assumed to be $d$-dimensional spheres. As discussed in Sec. \[sec:interface\], the characteristic linear size $\ell$ of the clusters increases linearly with time. In this case, the radius of a spherical cluster thus follows $R(t) = v t$, where $v$ is the growth velocity. The correlation function $G(r, t) = 0$ for $r > 2vt$, and is given by [@sekimoto86] $$\begin{aligned}
    \label{eq:gr}
    G(r, t) =4  [1- \phi(t)]^2  \!\left\{\exp\!\left[\mathcal{I}\, v^d t^{d+1} \Psi_d\!\left(\frac{r}{2vt}\right)\right] - 1\right\}, \,\,\end{aligned}$$ for $r < 2 vt$. Here $\phi(t)$ is the transformed volume fraction in Eq. (\[eq:avrami\]), and the function $\Psi_d(x)$ in 2D is $$\begin{aligned}
    \Psi_2(x) \!=\! \frac{2}{3} \!\left[ \cos^{-1}\! x - 2 x \sqrt{1-x^2} + x^3 \ln\! \left(\frac{1 + \sqrt{1-x^2}}{x} \right) \right]\!. \nonumber\end{aligned}$$ It is worth noting that this analytical expression agrees very well with the simulation results of the kinetic Ising model [@ramos99]. After absorbing the prefactor into the exponential, we have $\tilde G(r) = \exp[C \, \Psi_2(r/\zeta)] - 1$. The two independent fitting parameters, $C$ and $\zeta$, are determined from the conditions $\tilde G(0) = 1$ and that the gradient $d\tilde G/dr|_0$ equals the slope of the linear segment of the numerical correlation function at small $r$. As shown in Fig. \[fig:gr\](b), although the short-range correlation of the two-phase states can be satisfactorily described by the KJMA model, the correlation function is appreciably enhanced at larger distances.

![(Color online) \[fig:gr\_scaled\] Data collapsing from the normalized $g(r)$ versus rescaled distance $r/\xi$. Also shown is the analytical correlation function based on the KJMA model $\tilde{G}_{\rm KJMA}(r) = \exp[C \Psi_2(r/\zeta)] - 1$. The inset shows the length scaling factor $\xi$ (in arbitrary units) versus the Mott volume fraction. ](gr_scaled.pdf){width="0.99\columnwidth"}

Conclusion and outlook {#sec:outlook}
======================

To summarize, we have systematically studied the kinetics of the first-order Mott transition by applying a novel Gutzwiller molecular dynamics (MD) method to the Hubbard model. Contrary to phenomenological approaches, such as the kinetic Ising model or time-dependent Ginzburg-Landau equation, often used in the study of Mott metal-insulator transition, our Gutzwiller MD simulation offers, for the first time, a microscopic description of the phase transformation dynamics. In this microscopic approach, the Langevin dynamics is used to model the lattice fluctuations, which give rise to the relaxation dynamics and the stochastic force for the order-parameter field represented by the coarse-grained double-occupancy density. To include the crucial electron correlation effects, the electron degrees of freedom are integrated out [*on-the-fly*]{} by solving the Hubbard model with the efficient Gutzwiller method at every time-step.

Our MD simulations reveal that the metal-to-insulator transformation is initiated by the nucleation of Mott droplets in a metallic background. The initial fluctuations of the Mott clusters are well described by the classical droplet model. After an incubation period, whose duration decreases with increasing depth of the temperature quench, the system enters a nonequilibrium quasi-steady state with constant nucleation rate, consistent with the picture of the classical nucleation theory. The total number of the nuclei as well as their characteristic linear size increase linearly with time in this constant nucleation regime. This in turns points to a underlying interface-controlled cluster-growth mechanism, i.e. the growth rate of the nuclei is proportional to the surface area, which is also directly confirmed in our MD simulations.

The time evolution of the transformed Mott volume fraction $\phi_M$ follows a universal function of the scaled time, after taking into account the initial incubation period. Interestingly, the classical theory of phase transition kinetics developed by Kolmogorov, Johnson, Mehl, and Avrami (KJMA model) only describes the early constant-nucleation stage of this universal kinetics of the Mott transition. As the Mott volume fraction exceeds $\phi_M^\dagger \approx 0.15$, the transformation speed picks up compared with the KJMA model. This accelerated transformation stage is also characterized by intermittent bursts of large-size Mott clusters and a nonlocal nature of the cluster growth. We argue that this avalanche behavior, similar to the Barkhausen effect, could be attributed to self-generated disorder for the order-parameter field. The effects of the emergent disorder in the mixed-phase states on the electron wavefunction has largely been unexplored; such study will help understanding the mechanism of the unusual transition dynamics. One might draw insight from the extensive research on the interplay between Anderson localization and Mott transition [@dobrosavljevic97; @dobrosavljevic03; @byczuk05; @aguiar09; @tanaskovic03; @andrade09; @andrade10]. We hope future studies will shed light on the peculiar relation Eq. (\[eq:exp\_phi\]) that describes an exponential dependence of the average cluster size on the Mott volume fraction in this dynamical regime.

Finally, we have also characterize the structure of the transient mixed-phase states. The structure factor shows the telltale Porod’s power-law at large wavevectors, indicating sharp metal-insulator interfaces in the phase-separated states. A universal two-point correlation function is obtained after properly scaling the distance with respect to the correlation length. Although a power-law cluster-size distribution emerges at roughly the midpoint of the phase transition, the correlation length remains finite through out the transformation.

A detailed study of the the metal-insulator coexistence regime in VO$_2$ was recently carried out by applying the critical cluster techniques to the experimental near-field images [@liu16]. Interestingly, a robust power-law cluster-size distribution $\mathcal{D}(s) \sim 1/s^\tau$ was observed over a wide temperature range with an exponent $\tau = 1.92 \pm 0.11$ [@liu16]. Other geometrical characterizations such as the volume fractal dimension and hull fractal dimension were also computed. After comparing these exponents with that of potential universality classes, the authors of Ref. [@liu16] concluded that the MIT in VO$_2$ is consistent with the system passing near criticality of the random field Ising model (RFIM) [@sethna93; @sethna05] as temperature is varied. Although the criticality of the Mott transition is expected to be the Ising universality class in the absence of quenched disorder, as in our case, the results of our cluster analysis in Sec. \[sec:nucleation\] seem to be consistent with both scenarios. We plan to clarify this issue and to study the crossover phenomena between the thermal and disorder-dominated Mott transitions by carrying out similar Gutzwiller MD simulations on Hubbard models with quenched randomness in the future.

Another recent nano-imaging experiment on V$_2$O$_3$ offers a first dynamical picture of a photo-induced Mott insulator-to-metal transition in real materials [@ronchi18]. The volume fraction of the transformed metallic phase can be estimated from the reflectivity data. Interestingly, a two-stage transformation was observed from the time evolution of the metallic volume fraction. The initial stage is well described by the KJMA model with an exponent $n = 2$, suggesting a 2D domain growth with pre-existing nuclei which are quickly exhausted [@christian75]. Contrary to our case, the transformation after this initial stage is slowed down. The reduced Avrami exponent, ranging from 0 to $n \approx 3/4$ depending on the temperature, was suggested to result from a martenistic transition characterized by highly anisotropic domain growth [@cahn56]. To shed light on these intriguing results, we plan to conduct a systematic study on the kinetics of Mott transition in the reverse direction, i.e. the transformation from the Mott insulator to the correlated metal. Another intriguing result is the striped domains in the phase separated states of V$_2$O$_3$ [@mcleod16; @ronchi18]. Similar stripe patterns have also been observed in the mixed-phase states of VO$_2$ [@liu13]. It is likely that structural transition or dimerization, which are not included in our model, plays an important role in the stripe pattern and the martenistic transformation dynamics.

The lattice dynamics of the Hubbard-Peierls model in this work is highly simplified. In future study, we plan to apply the Gutzwiller MD to more realistic models. One approach is to incorporate the Hubbard-type interactions into the so-called tight-binding MD method (TBMD) [@khan89; @wang89; @goedecker94; @koskinen09]. In this semi-empirical formalism, the effective tight-binding Hamiltonian is determined from the instantaneous atomic configuration through a hopping parameter $t(r_{ij})$ and a classical two-body potential $\phi(r_{ij})$, both dependent on the interatomic distance. For example, an $s$-orbital TBMD was employed in our previous work [@chern17] to construct a liquid Hubbard model. However, stabilization of a crystal phase requires anisotropic hopping with higher-order orbitals. Such TBMD-Hubbard approach without the restriction to an underlying fixed lattice also make it possible to investigate the effects and dynamics of structural defects in a Mott transition.

Finally, we remark that the limitations of the Gutzwiller and the slave-boson methods are well known. Notably among them are the neglect of quantum fluctuations and its inability to describe the incoherent electronic excitations. Moreover, while the spatial fluctuations and inhomogeneity are included in the real-space Gutzwiller method, as used in our simulations, the inter-site correlations are ignored. This drawback is shared by the single-site DMFT. While multi-site or cluster Gutzwiller has been proposed [@ayral17; @lee19], a self-consistent treatment is required for application to inhomogeneous electronic states. On the other hand, the Gutzwiller approach as a mean-field theory of the Mott transition offers a natural framework for multi-scale modeling of the correlated electron system. Also importantly, the GA is probably one of the few many-body techniques that can be feasibly combined with large-scale MD simulations. To go beyond the GA, the real-space DMFT [@potthoff99; @freericks04; @helmes08] as well as variational Monte Carlo (VMC) can be potentially used as the many-body solver of the quantum MD simulations. In particular, VMC has been used to study the Mott transition in disordered Hubbard model [@pezzoli09; @pezzoli10]. There have also been remarkable developments in the VMC-MD method and its applications in recent years [@magro96; @delaney06; @attaccalite08; @mazzola12]. One interesting direction for the future study would be to apply the VMC-MD to the same model Hamiltonian considered in this work. The time-consuming many-body calculations required in these sophisticated methods naturally limit the lattice sizes that can be feasibly simulated. As demonstrated in our recent works [@ma19; @suwa19], combining the advanced many-body techniques, such as real-space DMFT or VMC, with the machine learning methods offers a promising solution for large-scale dynamical simulations of correlated electron systems.

The author thanks C. D. Batista, K. Barros, and G. Kotliar for fruitful discussions and collaborations on related projects. The author also thanks K. Barros for the help with the kernel polynomial method, and N. Lanata and V. Dobrosavljević for insightful discussions on the real-space Gutzwiller/slave-boson calculation. This work is partially supported by the Center for Materials Theory as a part of the Computational Materials Science (CMS) program, funded by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

Gutzwiller/Slave Boson method {#sec:GA}
=============================

Here we outline details of the Gutzwiller solver used in the MD simulations. In terms of the slave boson coherent state $\Phi_i = (e_i, p_{i\uparrow}, p_{i\downarrow}, d_i)$ and the quasi-particle density matrix $\rho^{\rm qp}$, the free energy in Eq. (\[eq:free-energy\]) becomes $$\begin{aligned}
    \mathcal{F}[\rho^{\rm qp}, \Phi] &=& \sum_{\langle ij \rangle}\sum_\alpha \left(t_{ij}  \mathcal{R}^{\;}_{i\alpha} \mathcal{R}^*_{j\alpha}\, \rho_{j\alpha,i\alpha} 
    + t_{ji} \mathcal{R}^*_{i\alpha} \mathcal{R}^{\;}_{j\alpha} \, \rho_{i\alpha,j\alpha} \right) + U \sum_i D_i  
     + \,\, \sum_{i, \alpha} \mu_{i\alpha} \left(n_{i,\alpha} - \rho_{i\alpha, i\alpha} \right) \nonumber \\
    & & + T \, {\rm Tr}(\rho^{\rm qp} \ln \rho^{\rm qp}) 
     + \,\, T \sum_i \left[ |e_i|^2 \ln\frac{|e_i|^2}{\Pi_{i,e}} + |p_{i\uparrow}|^2 \ln \frac{|p_{i \uparrow}|^2}{\Pi_{i,\uparrow}} + + |p_{i,\downarrow}|^2 \ln \frac{|p_{i \downarrow}|^2}{\Pi_{i\downarrow}}
    + |d_i|^2 \ln \frac{|d_i|^2}{\Pi_{i, d}} \right].\end{aligned}$$ Here we have defined the two-point correlation function, $\rho_{i\alpha, j\beta} \equiv \langle c^\dagger_{j\beta} c^{\;}_{i\alpha} \rangle = {\rm Tr}(\rho^{\rm qp}\, c^\dagger_{j\beta} c^{\;}_{i\alpha} )$, which is also called the single-particle reduced density matrix. The particle number $n_{i\alpha}$, the local double occupancy, And the renormalization factors are given by \[cf. Eq. (\[eq:R\])\] $$\begin{aligned}
    \label{eq:R_2}
    n_{i, \alpha} = |p_{i, \alpha}|^2 + |d_i|^2, \qquad D_i = |d_i|^2, \qquad 
    \mathcal{R}_{i,\alpha} = \frac{ p^*_{i, \alpha} e^{\;}_i + d_i^* p^{\;}_{i, \overline{\alpha}} }{\sqrt{n_{i\alpha} (1 - n_{i\alpha})}}.\end{aligned}$$ The symbol $\overline{\alpha}$ denotes the time-reversed spin direction of $\alpha$. The uncorrelated local distribution $\Pi^0_i$ can be expressed in terms of the local electron density $$\begin{aligned}
    \label{eq:P0}
    \Pi_{i, e} = (1-n_{i, \uparrow})(1-n_{i, \downarrow}), \quad \Pi_{i, \uparrow} = n_{i, \uparrow} (1- n_{i, \downarrow}), \quad 
    \Pi_{i, \downarrow} = n_{i, \downarrow} (1- n_{i, \uparrow}), \quad \Pi_{i, d} = n_{i,\uparrow} n_{i, \downarrow}.\end{aligned}$$ In particular, both $\mathcal{R}_{i\alpha}$ and $\Pi_i$ should be viewed as functions of the slave boson amplitudes through $n_{i\alpha}$. Minimization of the $\mathcal{F}$ with respect to the quasi-particle density matrix amounts to solving the eigenvalue problem $H^{(1)} \varphi^{(m)} = \varepsilon_m \varphi^{(m)}$, where $H^{(1)}$ is the “first-quantized" matrix of the renormalized quasi-particle Hamiltonian given in Eq. (\[eq:H\_renorm\]), i.e. $\mathcal{H}^{\rm qp} = \sum_{i\alpha, j\beta} c^\dagger_{i\alpha} H^{(1)}_{i\alpha, j\beta} c^{\;}_{j\beta}$. The reduced density matrix can be expressed as $\rho = f_{\rm FD}(H^{(1)} - \mu)$, or explicitly in terms of the real-space eigenstates: $$\begin{aligned}
    \rho_{i\alpha, j\beta} = \sum_{m} f_{\rm FD}(\varepsilon_m - \mu) \varphi^{(m)*}_{i\alpha}\varphi^{(m)}_{j\beta},  \end{aligned}$$ where $f_{\rm FD}(\varepsilon)$ is the Fermi-Dirac function. However, calculation of $\rho$ based on exact diagonalization is infeasible for large-scale simulations. Instead, here we adopt the linear-scaling kernel polynomial method (KPM) [@barros13; @wang18] to solve the quasi-particle density matrix.

The minimization with respect to the variational variables can be achieved by solving the coupled equations $\partial \mathcal{F} / \partial \Phi^\dagger_i = 0$ subject to the normalization condition $ |e_i|^2 + |p_{i\uparrow}|^2 + |p_{i \downarrow}|^2 + |d_i|^2 = 1$. Explicit calculations lead to the following expressions $$\begin{aligned}
    \label{eq:dFde}
    \frac{\partial \mathcal{F}}{\partial e^*_i} = \sum_{\beta = \uparrow, \downarrow}  \Delta^*_{i,\beta} \frac{\partial \mathcal{R}^*_{i,\beta}}{\partial e^*_i} 
    + T \left( \ln\frac{|e_i|^2}{\Pi_{i, e}} + 1\right) e_i = \sum_{\beta = \uparrow, \downarrow} \frac{\Delta^*_{i,\beta}\, p_{i,\beta}}{\sqrt{n_{i\beta} (1 - n_{i\beta})}} + T \left( \ln\frac{|e_i|^2}{\Pi_{i, e}} + 1\right) e_i, \end{aligned}$$ $$\begin{aligned}
    \label{eq:dFdp}
    \frac{\partial \mathcal{F}}{\partial p^*_{i\alpha}} &=&  \Delta_{i,\alpha} \frac{\partial \mathcal{R}_{i,\alpha}}{\partial p^*_{i\alpha}}
    + \Delta^*_{i, \overline{\alpha}} \frac{\partial \mathcal{R}^*_{i, \overline{\alpha}}}{\partial p^*_{i\alpha}}
    + \left( \Delta_{i,\alpha} \frac{\partial \mathcal{R}_{i,\alpha}}{\partial n_{i\alpha}} + \Delta^*_{i, \alpha} \frac{\partial \mathcal{R}^*_{i, \alpha}}{\partial n_{i\alpha}} \right) \frac{\partial n_{i\alpha}}{\partial p^*_{i \alpha}} 
    + \mu_{i, \alpha} p_{i\alpha}
    + T \left( \ln\frac{|p_{i \alpha}|^2}{\Pi_{i, \alpha}} + 1\right) p_{i\alpha} \nonumber \\
    & & - T \left( |e_i|^2 \frac{\partial \ln \Pi_{i, e}}{\partial n_{i,\alpha}} 
    +  |p_{i\uparrow}|^2 \frac{\partial \ln \Pi_{i, \uparrow}}{\partial n_{i,\alpha}} 
    +  |p_{i\downarrow}|^2 \frac{\partial \ln \Pi_{i, \downarrow}}{\partial n_{i,\alpha}} 
    +  |d_i|^2 \frac{\partial \ln \Pi_{i, d}}{\partial n_{i,\alpha}}  \right) \frac{\partial n_{i,\alpha}}{\partial p^*_{i \alpha}} \nonumber \\
    & = & \frac{\Delta_{i, \alpha} \, e_i}{\sqrt{n_{i\alpha} (1- n_{i\alpha})} }+ \frac{ \Delta^*_{i, \overline{\alpha}} \, d_i}{\sqrt{n_{i \overline{\alpha}} (1- n_{i \overline{\alpha}})}}
    + \left( \Lambda_{i, \alpha} + \mu_{i, \alpha} \right) \,p_{i \alpha} 
    + T \left( \ln\frac{|p_{i \alpha}|^2}{\Pi_{i, \alpha}} + 1\right) p_{i\alpha},\end{aligned}$$ $$\begin{aligned}
    \label{eq:dFdd}
    \frac{\partial \mathcal{F}}{\partial d^*_i} &=& \sum_{\beta = \uparrow, \downarrow} \left[ \Delta_{i, \beta} \frac{\partial \mathcal{R}_{i, \beta}}{\partial d^*_i} 
    +  \left(\Delta_{i, \beta} \frac{\partial \mathcal{R}_{i, \beta}}{\partial n_{i, \beta}} + \Delta^*_{i, \beta} \frac{\partial \mathcal{R}^*_{i, \beta}}{\partial n_{i, \beta}} \right) 
    \frac{\partial n_{i, \beta}}{\partial d^*_i}
    +   \mu_{i,\beta} \, d_i \right]
    + U d_i + T \left( \ln\frac{|d_i|^2}{\Pi_{i, d}} + 1\right) d_i \nonumber \\
    & & - T \sum_{\beta = \uparrow, \downarrow} \left( |e_i|^2 \frac{\partial \ln \Pi_{i, e}}{\partial n_{i,\beta}} 
    +  |p_{i\uparrow}|^2 \frac{\partial \ln \Pi_{i, \uparrow}}{\partial n_{i,\beta}} 
    +  |p_{i\downarrow}|^2 \frac{\partial \ln \Pi_{i, \downarrow}}{\partial n_{i,\beta}} 
    +  |d_i|^2 \frac{\partial \ln \Pi_{i, d}}{\partial n_{i,\beta}}  \right) \frac{\partial n_{i,\beta}}{\partial d^*_{i}} \nonumber \\
    &=& \sum_{\beta = \uparrow, \downarrow} \left[ \frac{\Delta_{i, \beta} \, p_{i, \overline{\beta}}}{\sqrt{n_{i\beta}(1-n_{i\beta})}} + \left(\Lambda_{i,\beta} + \mu_{i,\beta} \right) \, d_i \right]
    + U d_i + T \left( \ln\frac{|d_i|^2}{\Pi_{i, d}} + 1\right) d_i\end{aligned}$$ where $\Delta_{i,\alpha}$ is given in Eq. (\[eq:Delta\]), and $\Lambda_{i,\beta}$ is an effective chemical potential $$\begin{aligned}
    \Lambda_{i, \alpha} =  \frac{ \left(n_{i \alpha} - \frac{1}{2} \right) ( \Delta_{i, \alpha} \mathcal{R}_{i, \alpha} +  \Delta^*_{i, \alpha} \mathcal{R}^*_{i, \alpha}) }{n_{i\alpha} (1- n_{i\alpha})}
    + T \Bigl[  |e_i|^2 \, (1 - n_{i, \overline{\alpha}}) - |p_{i, \alpha}|^2 (1-n_{i,\alpha}) + |p_{i, \overline{\alpha}}|^2 n_{i, \overline{\alpha}} - |d_i|^2 n_{i, \overline{\alpha}} \Bigr]. \end{aligned}$$ The set of coupled equations (\[eq:dFde\])–(\[eq:dFdd\]) is highly nonlinear and has to be solved with hard constraints $|\Phi_i|^2 = 1$. A more convenient and robust approach, as outlined in Ref. [@lanata12], is to reformulate the minimization to a set of nonlinear eigenvalue problem $\mathcal{H}^{\rm SB}_i[\Phi_i] \, \Phi_i = \epsilon_i \Phi_i$, where the effective slave-boson Hamiltonian $\mathcal{H}_i^{\rm SB}$ is $$\begin{aligned}
    \mathcal{H}^{\rm SB}_i = \sum_\alpha \frac{ \left(\Delta_{i, \alpha} \mathcal{M}_{\alpha} + \Delta^*_{i, \alpha} \mathcal{M}^{\dagger}_{\alpha} \right) }{\sqrt{n_{i\alpha} (1- n_{i\alpha})}} + U \mathcal{U} 
    + \sum_{\alpha} \left(\Lambda_{i,\alpha} + \mu_{i, \alpha} \right) \, \mathcal{N}_{\alpha} + T \Theta_{i}.\end{aligned}$$ This is the same as Eq. (\[eq:H\_SB\]) in the main text. The various matrices in the above expression are $$\begin{aligned}
    \mathcal{M}_{\uparrow} = \left[ \begin{array}{cccc}
    0 & 0 & 0 & 0 \\
    1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 
    \end{array} \right],  \quad
    \mathcal{M}_{\downarrow} = \left[ \begin{array}{cccc}
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 
    \end{array} \right],  \quad
    \mathcal{U} = \left[ \begin{array}{cccc}
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 
    \end{array} \right], \nonumber\end{aligned}$$ $$\begin{aligned}
    \mathcal{N}_{\downarrow} = \left[ \begin{array}{cccc}
    0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 
    \end{array} \right],  \quad
    \mathcal{N}_{\downarrow} = \left[ \begin{array}{cccc}
    0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 1 
    \end{array} \right], \quad
    \Theta_i = \mathbb{I}_4 +  \left[ \begin{array}{cccc}
    \ln(|e_i|^2/\Pi_{i, e}) & 0 & 0 & 0 \\
    0 & \ln(|p_{i\uparrow}|^2/\Pi_{i, \uparrow}) & 0 & 0 \\
    0 & 0 & \ln(|p_{i\downarrow}|^2/\Pi_{i, \downarrow}) & 0 \\
    0 & 0 & 0 & \ln(|d_i|^2/\Pi_{i, d}) 
    \end{array} \right]. \qquad\end{aligned}$$ The eigen-solutions of the renormalized tight-binding Hamiltonian $\mathcal{H}^{\rm qp}$ and the embedding Hamiltonians $\mathcal{H}^{\rm SB}_i$ have to be solved self-consistently. This can be achieved through the iteration procedure. Explicitly, let the solution after $k$-iteration be $\rho^{(k)}_{i\alpha, j\beta}$, $\Phi^{(k)}_i$, and $\mu^{(k)}_{i\alpha}$. The $(k+1)$-th iteration consists of the following steps: (i) New renormalization factors $\mathcal{R}_{i, \alpha}^{(k)}$ are computed from the slave-boson amplitudes $\Phi^{(k)}_i$ using Eq. (\[eq:R\_2\]). A new tight-binding Hamiltonian $\mathcal{H}^{(k+1)}_{\rm qp} = \mathcal{H}^{\rm qp}[\mathcal{R}^{(k)}, \mu^{(k)}]$ in turn is built from these renormalization factors. (ii) KPM method is used to solve the quasi-particle density matrix $\rho^{(k+1)}_{i\alpha, j\beta}$ from the Hamiltonian $\mathcal{H}^{(k+1)}_{\rm qp}$. (iii) A new set of kinetic factors $\Delta_{i, \alpha}^{(k+1)}$ is calculated from $\rho^{(k+1)}_{i\alpha,j\beta}$ using Eq. (\[eq:Delta\]). The new embedding Hamiltonians are then constructed according to $\mathcal{H}^{(k+1)}_{i, \rm SB} = \mathcal{H}^{\rm SB}_i[\Phi^{(k)}, \Delta^{(k+1)}, \mu^{(k)}]$. (iv) Exact diagonalization of $\mathcal{H}^{(k+1)}_{i, \rm SB}$ then gives a new set of slave-boson amplitudes $\Phi^{(k+1)}_i$. Finally, new values of the Lagrangian multipliers $\mu^{(k+1)}_{i,\alpha}$ are determined in order to satisfy the Gutzwiller constraints $n_{i\alpha} = \rho_{i\alpha, i\alpha}$. These steps are repeated until convergence is reached.

Paramagnetic Gutzwiller solutions {#sec:paramagnetic}
=================================

For non-magnetic solutions, we assume $p_{i, \uparrow} = p_{i, \downarrow}$, and $\rho_{i\alpha, j\beta} = \rho_{ij} \delta_{\alpha\beta}$. Since all derived quantities are now spin-independent, we will drop the spin indices. Futhermore, one can introduce a 3-vector for the slave-boson wavefunction $\Phi_i = (\phi_{i, 0}, \phi_{i, 1}, \phi_{i, 2}) = (e_i, \sqrt{2} p_i, d_i)$, where $p_{i} \equiv p_{i, \uparrow} = p_{i, \downarrow}$. The new wavefunction thus satisfies the same normalization condition $|\Phi_i|^2 = 1$. Using the appropriate gauge, both the renormalization factors $\mathcal{R}_i$ and the kinetic $\Delta_i$ can be made real. The free-energy now becomes $$\begin{aligned}
    \mathcal{F}[\rho^{\rm qp}, \Phi] = 2 \sum_{\langle ij \rangle}\mathcal{R}_i \mathcal{R}_j \, t_{ij} (\rho_{ij} + \rho_{ji}) + U \sum_i D_i 
    + 2 \sum_i \mu_i (n_i - \rho^{\;}_{ii}) 
    + 2T {\rm Tr}(\rho^{\rm qp} \ln \rho^{\rm qp}) + T \sum_{i}\sum_{m=0}^2 \phi_{i, m}^2 \ln\frac{\phi_{i,m}^2}{\Pi_{i,m}}, \quad\end{aligned}$$ where we have defined $\Pi_{i, 0} = \Pi_{i, e}$, $\Pi_{i, 1} = 2 \Pi_{i, \alpha}$ ($\alpha = \uparrow, \downarrow$), and $\Pi_{i, 2} = \Pi_{i, d}$. In terms of the $\phi$ variables, other relevant quantities are $$\begin{aligned}
    n_i = \frac{1}{2} \phi_{i, 1}^2 + \phi_{i, 2}^2, \qquad D_i = \phi_{i, 2}^2, \qquad \mathcal{R}_i = \frac{\phi_{i, 0} \phi_{i, 1} + \phi_{i, 1} \phi_{i, 2}}{\sqrt{2 n_i (1- n_i)}}.\end{aligned}$$ By expressing the $\phi_{i,m}$ in terms of $n_i$ and $D_i$ using the first two equations, one arrives at the Gutzwiller formula for the renormalization factor. The optimization of the quasi-particle density again amounts to solving the renormalized spinless Hamiltonian $$\begin{aligned}
    \mathcal{H}^{\rm qp} = \sum_{\langle ij \rangle} \mathcal{R}_i \mathcal{R}_j \, t_{ij} \left( c^\dagger_i c^{\;}_j + c^\dagger_j c^{\;}_i \right) - \sum_i \mu_i c^\dagger_i c^{\;}_i.\end{aligned}$$ And similarly, the minimization with respect to the $\phi$ variables subject to condition $|\Phi_i|^2 = 1$ can be recast into the eigen-solution of the following $3\times 3$ Hamiltonian $$\begin{aligned}
    \mathcal{H}^{\rm SB}_i = \frac{\sqrt{2}\Delta_i \mathcal{M}}{\sqrt{n_i (1-n_i)}} + U \mathcal{U} + (\Lambda_i + \mu_i) \mathcal{N} + T\Theta_i,\end{aligned}$$ where the various matrices now become $$\begin{aligned}
    \mathcal{M} = \left[ \begin{array}{ccc}
    0 & 1 & 0 \\
    1 & 0 & 1 \\
    0 & 1 & 0 
    \end{array} \right], 
    \quad
    \mathcal{U} = \left[ \begin{array}{ccc}
    0 & 0 & 0 \\
    0 & 0 & 0 \\
    0 & 0 & 1 
    \end{array} \right], 
    \quad
    \mathcal{N} = \left[ \begin{array}{ccc}
    0 & 0 & 0 \\
    0 & 1 & 0 \\
    0 & 0 & 2 
    \end{array} \right], 
    \quad
    \Theta_i = \mathbb{I}_3 + \left[ \begin{array}{ccc}
    \ln(\phi_{i,0}^2/\Pi_{i, 0}) & 0 & 0 \\
    0 & \ln(\phi_{i,1}^2/\Pi_{i, 1}) & 0 \\
    0 & 0 & \ln(\phi_{i,2}^2/\Pi_{i, 2}) 
    \end{array} \right], \quad\end{aligned}$$

Effective spin model {#sec:spin-model}
====================

In this Appendix, we discuss the effective spin-1/2 Hamiltonian of the Gutzwiller method based on the assumption of particle-hole symmetry and negligible charge fluctuations. This effective spin model can also be obtained using the slave-spin method for the Hubbard model. Here we consider the paramagnetic solutions with $n_{i\uparrow} = n_{i\downarrow} = 1/2$. The slave-boson amplitudes can now be represented using a pseudo-spinor $\Phi_i = (\phi_{i,1}, \phi_{i,2})$, where $\phi_{i, 1} = \sqrt{2} p_{i, \uparrow} = \sqrt{2} p_{i, \downarrow}$ and $\phi_{i, 2} = \sqrt{2} e_i = \sqrt{2} d_i$. Next we introduce spin-1/2 operator $\hat{\mathbf S} = \frac{1}{2} \bm\tau$, where $\bm\tau = [\tau_x, \tau_y, \tau_z]$ is the vector of Pauli matrices. The renormalization factors and the double-occupancy can both be written as the expectation values $$\begin{aligned}
    \mathcal{R}_i = 2 \phi_{i,0} \phi_{i,1} = \langle \Phi_i | \tau_x | \Phi_i \rangle \equiv 2  S^x_i, \qquad
    D_i = d_i^2 = \frac{1}{2} \phi_{i,2}^2 = \frac{1}{4} \left[ 1 - \langle \Phi_i | \tau_z | \Phi_i \rangle \right] \equiv \frac{1}{4} ( 1 - 2 S^z_i ).\end{aligned}$$ Here we have introduced classical spin variables $\mathbf S_i \equiv \frac{1}{2} \langle \Phi_i | \bm\tau | \Phi_i \rangle$. The free-energy can now be expressed as $\mathcal{F} = \mathcal{H}_{\rm eff} - T \mathcal{S}$, where the effective spin Hamiltonian and the entropy term are $$\begin{aligned}
    \mathcal{H}_{\rm eff} = \sum_{\langle ij \rangle} \mathcal{J}_{ij} \,S^x_i S^x_j - H \sum_i S^z_i,\end{aligned}$$ $$\begin{aligned}
    \mathcal{S} = - \sum_i \left[ \left(\frac{1}{2}+S^z_i \right) \ln\left(\frac{1}{2}+S^z_i \right) + \left(\frac{1}{2}-S^z_i \right) \ln\left(\frac{1}{2}-S^z_i \right) \right].\end{aligned}$$ The zero-temperature part here is essentially the quantum Ising model with the “transverse" field $H = U/2$. The effective exchange constant is $$\begin{aligned}
    \mathcal{J}_{ij} = 8 t_{ij} (\rho_{ij} + \rho_{ji}).\end{aligned}$$ The electrons are described by a renormalized tight-binding Hamiltonian $\mathcal{H}^{\rm qp} = \sum_{\langle ij \rangle} t^{\rm eff}_{ij} (c^\dagger_i c^{\;}_j + \mbox{h.c.})$, where the effective hopping depends on the spin configuration, i.e. $t^{\rm eff}_{ij} = 4 S^x_i\, S^x_j \, t_{ij}$. In the spin picture, the order parameter of the Mott transition is the $z$-component of the average spin; the first-order Mott transition is marked by a jump in $\langle S^z \rangle$ when crossing the transition temperature from the correlated metal to the Mott phase. Importantly, the nearest-neighbor exchange constant $\mathcal{J}_{ij}$ depends implicitly on spins at further neighbors through the nonlocal electronic degrees of freedom $\rho_{ij}$.

[99]{}

N. F. Mott, [*Metal-Insulator Transitions*]{} (Taylor & Francis, New York, 1997).

M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. [**70**]{}, 1039 (1998).

V. Dobrosavljevic, [*Introduction to Metal-Insulator Transitions*]{} (Oxford University Press, Oxford, 2012).

V. J. Emery and S. A. Kivelson, Superconductivity in Bad Metals, Phys. Rev. Lett. [**74**]{}, 3253 (1995).

P. Limelette, A. Georges, D. Jérome, P. Wzietek, P. Metcalf, J. M. Honig, Universality and Critical Behavior at the Mott Transition, Science [**302**]{}, 89 (2003).

L. Bartosch, M. de Souza, and M. Lang, Scaling Theory of the Mott Transition and Breakdown of the Grüneisen Scaling Near a Finite-Temperature Critical End Point, Phys. Rev. Lett. [**104**]{}, 245701 (2010).

T. Itou, E. Watanabe, S. Maegawa, A. Tajima, N. Tajima, K. Kubo, R. Kato, and K. Kanoda, Slow dynamics of electrons at a metal-Mott insulator boundary in an organic system with disorder, Sci. Adv. [**3**]{}, e1601594 (2017).

A. Pustogow, M. Bories, A. Löhle, R. Rösslhuber, E. Zhukova, B. Gorshunov, S. Tomić, J. A. Schlueter, R. Hübner, T. Hiramatsu, Y. Yoshida, G. Saito, R. Kato, T.-H. Lee, V. Dobrosavljević, S. Fratini, and M. Dressel, Quantum spin liquids unveil the genuine Mott state, Nature Mater. [**17**]{}, 773 (2018).

H. Terletska, J. Vucicević, D. Tanasković, and V. Dobrosavljević, Quantum Critical Transport near the Mott Transition, Phys. Rev. Lett. [**107**]{}, 026401 (2011).

G. Sordi, P. Sémon, K. Haule, and A.-M. S. Tremblay, Strong Coupling Superconductivity, Pseudogap, and Mott Transition, Phys. Rev. Lett. [**108**]{}, 216401(2012).

J. Vucicević, D. Tanasković, M. J. Rozenberg, and V. Dobrosavljević, Bad-Metal Behavior Reveals Mott Quantum Criticality in Doped Hubbard Models, Phys. Rev. Lett. [**114**]{}, 246402 (2015).

M. Abdel-Jawad, R. Kato, I. Watanabe, N. Tajima, Y. Ishii, Universality Class of the Mott Transition, Phys. Rev. Lett. [**114**]{}, 106401 (2015).

B. Lenz, S. R. Manmana, T. Pruschke, F. F. Assaad, and M. Raczkowski, Mott Quantum Criticality in the Anisotropic 2D Hubbard Model, Phys. Rev. Lett. [**116**]{}, 086403 (2016).

M. Imada, Universality classes of metal-insulator transitions in strongly correlated electron systems and mechanism of high-temperature superconductivity, Phys. Rev. B [**72**]{}, 075113 (2005).

J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. Kin-Lic Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Solutions of the Two-Dimensional Hubbard Model: Benchmarks and Results from a Wide Range of Numerical Algorithms, Phys. Rev. X [**5**]{}, 041041 (2015).

M. Motta, D. M. Ceperley, G. Kin-Lic Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Towards the Solution of the Many-Electron Problem in Real Materials: Equation of State of the Hydrogen Chain with State-of-the-Art Many-Body Methods, Phys. Rev. X [**7**]{}, 031059 (2017).

P. R. C. Kent and G. Kotliar, Toward a predictive theory of correlated materials, Special Section: Frontiers in Computation, Science [**361**]{}, 348 (2018).

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).

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).

E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. [**83**]{}, 349 (2011).

E. Dagotto, Complexity in strongly correlated electronic systems, Science [**309**]{}, 257 (2005).

D. A. Bonnell, D. N. Basov, M. Bode, U. Diebold, S. V. Kalinin, V. Madhavan, L. Novotny, M. Salmeron, U. D. Schwarz, and P. S. Weiss, Imaging physical phenomena with local probes: From electrons to photons, Rev. Mod. Phys. [**84**]{}, 1343

J. M. Atkin, S. Berweger, A. C. Jones, and M. B. Raschke, Nano-optical imaging and spectroscopy of order, phases, and domains in complex solids, Adv. Phys. [**61**]{}, 745 (2012).

M. Liu, A. J. Sternbach, and D. N. Basov, Nanoscale electrodynamics of strongly correlated quantum materials, Rep. Prog. Phys. [**80**]{}, 014501 (2017).

M. M. Qazilbash, M. Brehm, B.-G. Chae, P.-C. Ho, G. O. Andreev, B.-J. Kim, S. J. Yun, A. V. Balatsky, M. B. Maple, F. Keilmann, H.-T. Kim, and D. N. Basov, Mott transition in VO$_2$ Revealed by infrared spectroscopy and nano-imaging, Science [**318**]{}, 1750 (2007).

S. Lupi, L. Baldassarre, B. Mansart, A. Perucchi, A. Barinov, P. Dudin, E. Papalazarou, F. Rodolakis, J.-P. Rueff, J.-P. Itié, S. Ravy, D. Nicoletti, P. Postorino, P. Hansmann, N. Parragh, A. Toschi, T. Saha-Dasgupta, O. K. Andersen, G. Sangiovanni, K. Held and M. Marsi, A microscopic view on the Mott transition in chromium-doped V$_2$O$_3$, Nat. Commun. [**1**]{}, 105 (2010).

A. S. McLeod, E. van Heumen, J. G. Ramirez, S. Wang, T. Saerbeck, S. Guenon, M. Goldflam, L. Anderegg, P. Kelly, A. Mueller, M. K. Liu, Ivan K. Schuller, and D. N. Basov, Nanotextured phase coexistence in the correlated insulator V$_2$O$_3$, Nat. Phys. [**13**]{}, 80 (2016).

M. K. Liu, M. Wagner, E. Abreu, S. Kittiwatanakul, A. McLeod, Z. Fei, M. Goldflam, S. Dai, M. M. Fogler, J. Lu, S. A. Wolf, R. D. Averitt, and D. N. Basov, Anisotropic electronic state via spontaneous phase separation in strained vanadium dioxide films, Phys. Rev. Lett. [**111**]{}, 096602 (2013).

H. Madan, M. Jerry, A. Pogrebnyakov, T. Mayer, and S. Datta, Quantitative mapping of phase coexistence in Mott-Peierls insulator during electronic and thermally driven phase transition, ACS Nano [**9**]{}, 2009 (2015).

M. M. Qazilbash, A. Tripathi, A. A. Schafgans, Bong-Jun Kim, Hyun-Tak Kim, Zhonghou Cai, M. V. Holt, J. M. Maser, F. Keilmann, O. G. Shpyrko, and D. N. Basov, Nanoscale imaging of the electronic and structural transitions in vanadium dioxide, Phys. Rev. B [**83**]{}, 165108 (2011).

H. T. Stinson, A. Sternbach, O. Najera, R. Jing, A. S. Mcleod, T. V. Slusar, A. Mueller, L. Anderegg, H. T. Kim, M. Rozenberg, and D. N. Basov, Imaging the nanoscale phase separation in vanadium dioxide thin films at terahertz frequencies, Nat. Commun. [**9**]{}, 3604 (2018).

G. Mattoni, P. Zubko, F. Maccherozzi, A.J.H. van der Torren, D. B. Boltje, M. Hadjimichael, N. Manca, S. Catalano, M. Gibert, Y. Liu, J. Aarts, J.-M. Triscone, S. S. Dhesi, and A. D. Caviglia, Striped nanoscale phase separation at the metal-insulator transition of heteroepitaxial nickelates, Nat. Commun. [**7**]{}, 13141 (2016).

D. Preziosi, L. Lopez-Mir, X. Li, T. Cornelissen, J.-H. Lee, F. Trier, K. Bouzehouane, S. Valencia, A. Gloter, A. Barthélémy, and M. Bibes, Direct Mapping of Phase Separation across the Metal-Insulator Transition of NdNiO$_3$, Nano Lett. [**18**]{}, 2226 (2018).

A. Singer, J. G. Ramirez, I. Valmianski, D. Cela, N. Hua, R. Kukreja, J. Wingert, O. Kovalchuk, J. M. Glownia, M. Sikorski, M. Chollet, M. Holt, I. K. Schuller, and O. G. Shpyrko, Nonequilibrium Phase Precursors during a Photoexcited Insulator-to-Metal Transition in V$_2$O$_3$, Phys. Rev. Lett. [**120**]{}, 207601 (2018).

A. Ronchi, P. Homm, M. Menghini, P. Franceschini, F. Maccherozzi, F. Banfi, G. Ferrini, F. Cilento, F. Parmigiani, S. S. Dhesi, M. Fabrizio, J.-P. Locquet, and C. Giannetti, Nanoscale dynamics across the Mott transition in V$_2$O$_3$, arXiv:1807:03670 (2018).

M. Uehara, S. Mori, C. H. Chen, and S.-W. Cheong, Percolative phase separation underlies colossal magnetoresistance in mixed-valent manganites, Nature [**399**]{}, 560 (1999).

M. Fäth, S. Freisem, A. A. Menovsky, Y. Tomioka, J. Aarts, and J. A. Mydosh, Spatially inhomogeneous metal-insulator transition in doped manganites, Science [**285**]{}, 1540 (1999).

L. Zhang, C. Israel, A. Biswas, R. L. Greene, A. de Lozanne, Direct observation of percolation in a maganite thin film, Science [**298**]{}, 805 (2002).

M. C. Gutzwiller, Effect of Correlation on the Ferromagnetism of Transition Metals, Phys. Rev. Lett. [**10**]{}, 159 (1963).

M. C. Gutzwiller, Effect of Correlation on the Ferromagnetism of Transition Metals, Phys. Rev. [**134**]{}, A923 (1964).

M. C. Gutzwiller, Correlation of Electrons in a Narrow $s$ Band, Phys. Rev. [**137**]{}, A1726 (1965).

W. F. Brinkman and T. M. Rice, Application of Gutzwiller’s Variational Method to the Metal-Insulator Transition, Phys. Rev. B [**2**]{}, 4302 (1970).

J. S. Langer, Instabilities and pattern formation in crystal growth, Rev. Mod. Phys. [**52**]{}, 1 (1980).

J. D. Gunton, M. San Miguel, and P. S. Saint, [*The dynamics of first order phase transitions*]{}, in “Phase Transitions and Critical Phenomena”, edited by C. Domb and J. L. Lebowitz, vol. 8, pp.269-466 (Academic, New York, 1983).

S. Puri and V. Wadhawan, Ed. [*Kinetics of Phase Transitions*]{} (CRC Press, Taylor & Francis Group, London, 2009).

K. Binder, Theory of first-order phase transition, Rep. Prog. Phys. [**50**]{}, 783 (1987).

A. Onuki, [*Phase transition Dynamics*]{} (Cambridge, Cambridge Univ. Press 2002).

J. Marro, J. L. Lebowitz, and M. H. Kalos, Computer Simulation of the Time Evolution of a Quenched Model Alloy in the Nucleation Region, Phys. Rev. Lett. [**43**]{}, 282 (1979).

M. K. Phani, J. L. Lebowitz, M. H. Kalos, and O. Penrose, Kinetics of an Order-Disorder Transition, Phys. Rev. Lett. [**45**]{}, 366 (1980).

A. Sadiq and K. Binder, Kinetics of Domain Growth in Two Dimensions, Phys. Rev. Lett. [**51**]{}, 674 (1983).

D. Stauffer, A. Coniglio, and D. W. Heermann, Monte Carlo Experiment for Nucleation Rate in the Three-Dimensional Ising Model, Phys. Rev. Lett. [**49**]{}, 1299 (1982).

V. A. Shneidman, K. A. Jackson, and K. M. Beatty, Nucleation and growth of a stable phase in an Ising-type system, Phys. Rev. B [**59**]{}, 3579 (1999).

V. A. Shneidman, K. A. Jackson, and K. M. Beatty, On the applicability of the classical nucleation theory in an Ising model, J. Chem. Phys. [**111**]{}, 6932 (1999).

S. Wonczak, R. Strey, and D. Stauffer, Confirmation of classical nucleation theory by Monte Carlo simulations in the 3-dimensional Ising model at low temperature, J. Chem. Phys. [**113**]{}, 1976 (2000).

F. Soisson and G. Martin, Monte Carlo simulations of the decomposition of metastable solid solutions: Transient and steady-state nucleation kinetics, Phys. Rev. B [**62**]{}, 203 (2000).

K. Brendel, G. T. Barkema, and H. van Beijeren, Nucleation time in the two-dimensional Ising model, Phys. Rev. E [**71**]{}, 031601 (2005).

S. Ryu and W. Cai, Numerical tests of nucleation theories for the Ising models, Phys. Rev. E [**82**]{}, 011603 (2010).

K. Binder and P. Virnau, Overview: Understanding nucleation phenomena from simulations of lattice gas models, J. Chem. Phys. [**145**]{}, 211701 (2016).

S. Liu, B. Phillabaum, E. W. Carlson, K. A. Dahmen, N. S. Vidhyadhiraja, M. M. Qazilbash, and D. N. Basov, Random field driven spatial complexity at the Mott transition in VO$_2$, Phys. Rev. Lett. [**116**]{}, 036401 (2016).

T. Bar, S. K. Choudhary, M. A. Ashraf, K. S. Sujith, S. Puri, S. Raj, and B. Bansal, Kinetic spinodal instabilities in the Mott transition in V$_2$O$_3$: Evidence from hysteresis scaling and dissipative phase ordering, Phys. Rev. Lett. [**121**]{}, 045701 (2018).

S. Papanikolaou, R. M. Fernandes, E. Fradkin, P. W. Phillips, J. Schmalian, and R. Sknepnek, Universality of Liquid-Gas Mott Transitions at Finite Temperatures, Phys. Rev. Lett. [**100**]{}, 026408 (2008).

J. Lee and C.-H. Yee, Interfaces in coexisting metals and Mott insulators, Phys. Rev. B 95, 205126 (2017).

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).

N. Lanatá, H. U. R. Strand, Xi Dai, and Bo Hellsing, Efficient implementation of the Gutzwiller variational method, Phys. Rev. B [**85**]{}, 035133 (2012).

N. Lanatá, Y. Yao, C.-Z. Wang, K.-M. Ho, and G. Kotliar, Phase Diagram and Electronic Structure of Praseodymium and Plutonium, Phys. Rev. X [**5**]{}, 011008 (2015).

N. Lanatá, Y. Yao, X. Deng, V. Dobrosavljević, and G. Kotliar, Slave Boson theory of orbital differentiation with crystal field effects: Application to UO$_2$, Phys. Rev. Lett. [**118**]{}, 126401 (2017).

J. Bünemann, W. Weber, and F. Gebhard, “Multiband Gutzwiller Wave Functions for General On-Site Interactions", Phys. Rev. B [**57**]{}, 6896 (1998).

P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. [**49**]{}, 435 (1977).

K. Binder, Time-Dependent Ginzburg-Landau Theory of Nonequilibrium Relaxation, Phys. Rev. B [**8**]{}, 3423 (1973).

O. T. Valls and G. F. Mazenko, Nucleation in a time-dependent Ginzburg-Landau model: A numerical study, Phys. Rev. B [**42**]{}, 6614 (1990).

M. Iwamatsu, Direct numerical simulation of homogeneous nucleation and growth in a phase-field model using cell dynamics method, J. Chem. Phys. [**128**]{}, 084504 (2008).

R. Becker and W. Döring, Kinetische Behandlung der Keimbildung in übersättigten Dämpfen, Ann. Phys. [**24**]{}, 719 (1935).

Ya. B. Zeldovich, On the Theory of New Phase Formation: Cavitation, Acta Physicochem., USSR [**18**]{}, 1, (1943).

J. Frenkel, [*Kinetic Theory of Liquids*]{} (Oxford University Press, Oxford, 1946).

A. N. Kolmogorov, Statistical theory of nucleation processes, Bull. Acad. Sci. USSR, Sci. Mater. Nat. [**3**]{}, 355 (1937).

M. Avrami, Kinetics of Phase Change. I General Theory, J. Chem. Phys. [**7**]{}, 1103 (1939).

M. Avrami, Kinetics of Phase Change. II Transformation-Time Relations for Random Distribution of Nuclei, J. Chem. Phys. [**8**]{}, 212 (1940).

W. A. Johnson and R. F. Mehl, Reaction Kinetics in Processes of Nucleation and Growth, Trans. Metall. Soc. AIME, [**135**]{}, pp. 416-442 (1939).

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).

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).

D. Marx and J. Hutter, [*Ab initio molecular dynamics: basic theory and advanced methods*]{} (Cambridge University Press, Cambridge, 2009).

D. Plasienka, R. Martonak, and M. C. Newton, [*Ab initio*]{} molecular dynamics study of the structural and electronic transition in VO$_2$, Phys. Rev. B [**96**]{}, 054111 (2017).

G.-W. Chern, K. Barros, C. D. Batista, J. D. Kress, G. Kotliar, Mott transition in a metallic liquid - Gutzwiller molecular dynamics simulations, Phys. Rev. Lett. [**118**]{}, 226401 (2017)

M. Allen and D. Tildesley, [*Computer simulation of liquids*]{} (Oxford University Press, New York, 1989).

D. C. Rapaport, [*The art of molecular dynamics simulation,*]{} (Cambridge University Press, Cambridge, 2011).

M. Sandri, M. Capone, and M. Fabrizio, Finite-temperature Gutzwiller approximation and the phase diagram of a toy model for V$_2$O$_3$, Phys. Rev. B [**87**]{}, 205108 (2013).

W.-S. Wang, X.-M. He, D. Wang, Q.-H. Wang, Z. D. Wang, and F. C. Zhang, Finite-temperature Gutzwiller projection for strongly correlated electron systems, Phys. Rev. B [**82**]{}, 125105 (2010).

K. Barros and Y. Kato, Efficient Langevin simulation of coupled classical fields and fermions, Phys. Rev. B [**88**]{}, 235101 (2013).

Z. Wang, G.-W. Chern, C. D. Batista, and K. Barros, Gradient-based stochastic estimation of the density matrix, J. Chem. Phys. [**148**]{}, 094107 (2018).

G.-W. Chern, K. Barros, Z. Wang, H. Suwa, and C. D. Batista, Semiclassical dynamics of spin density waves, Phys. Rev. B [**97**]{}, 035120 (2018).

A. Mezio and R. H. McKenzie, Low quasiparticle coherence temperature in the one-band Hubbard model: A slave-boson approach, Phys. Rev. B [**96**]{}, 035121 (2017).

R. Bulla, T. A. Costi, D. Vollhardt, Finite-temperature numerical renormalization group study of the Mott transition, Phys. Rev. B [**64**]{}, 045103 (2001).

A. Liebsch, H. Ishida, and J. Merino, Mott transition in two-dimensional frustrated compounds, Phys. Rev. B [**79**]{}, 195108 (2009).

The slope of $dU_c/dT$ is shown to be positive for square lattice when inter-site correlations are included in the cluster DMFT calculations. However, as shown in Ref. [@liebsch09], cluster DMFT calculation on the frustrated triangular lattice finds a phase boundary with negative $dU_c/dT$.

D. Stauffer and A. Aharony, [*Introduction to Percolation Theory*]{}, 2nd ed. (CRC Press, 1994).

W. Janke and A. M. J. Schakel, Fractal structure of spin clusters and domain walls in the two-dimensional Ising model, Phys. Rev. E [**71**]{}, 036703 (2005).

J. W. Christian, [*The Theory of Transformations in Metals and Alloys*]{} (Pergamon Press, Oxford, 1975).

G. Durin and S. Zapperi, The Barkhausen effect, in [*The Science of Hysteresis*]{} Vol. II, edited by Giorgio Bertotti and I. Mayergoyz (Academic Press, 2005).

J. P. Sethna, K. A. Dahmen, and C. R. Myers, Crackling noise, Nature [**410**]{}, 242 (2001).

J. S. Urbach, R. C. Madison, and J. T. Markert, Interface depining, self-organized criticality, and the Barkhausen effect, Phys. Rev. Lett. [**75**]{}, 276 (1995).

P. Cizeau, S. Zapperi, G. Durin, and H. E. Stanley, Dynamics of a Ferromagnetic Domain Wall and the Barkhausen Effect, Phys. Rev. Lett. [**79**]{}, 4669 (1997).

O. Narayan, Self-similar Barkhausen noise in magnetic domain wall motion, Phys. Rev. Lett. [**77**]{}, 3855 (1996).

J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Hysteresis and hierarchies: Dynamics of disorder-driven first-order phase transformations, Phys. Rev. Lett. [**70**]{}, 3347 (1993).

J. P. Sethna, K. A. Dahmen, and O. Perković, Random-Field Ising Models of Hysteresis, in [*The Science of Hysteresis*]{} Vol. II, edited by Giorgio Bertotti and I. Mayergoyz (Academic Press, 2005).

G. Porod, Chap. 2 General theory in [*Small-Angle X-Ray Scattering*]{}, Eds. O. Glatter and O. Kratky (Academic Press, New York, 1982).

K. Sekimoto, Exact space correlation function corresponding to Avrami’s formula, J. Phys. Soc. Jpn. [**53**]{}, 2425 (1984).

K. Sekimoto, Evolution of the domain structure during the nucleation-and-growth process with non-conserved order parameter, Physica A [**135**]{}, 328 (1986).

R. A. Ramos, P. A. Rikvold, and M. A. Novotny, Test of the Kolmogorov-Johnson-Mehl-Avrami picture of metastable decay in a model with microscopic dynamics, Phys. Rev. B [**59**]{}, 9053 (1999).

V. Dobrosavljević and G. Kotliar, Mean Field Theory of the Mott-Anderson Transition, Phys. Rev. Lett. [**78**]{}, 3943 (1997).

V. Dobrosavljević, D. Tanasković, and A. A. Pastor, Glassy Behavior of Electrons Near Metal-Insulator Transitions, Phys. Rev. Lett. [**90**]{}, 016402 (2003).

K. Byczuk, W. Hofstetter, 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, G. Kotliar, Critical Behavior at the Mott-Anderson Transition: A Typical-Medium Theory Perspective, Phys. Rev. Lett. [**102**]{}, 156402 (2009).

D. Tanasković, V. Dobrosavljević, E. Abrahams, and G. Kotliar, Disorder Screening in Strongly Correlated Systems, Phys. Rev. Lett. [**91**]{}, 066603 (2003).

E. C. Andrade, E. Miranda, and V. Dobrosavljević, Electronic Griffiths Phase of the $d=2$ Mott transition, Phys. Rev. Lett. [**102**]{}, 206403 (2009).

E. C. Andrade, E. Miranda, and V. Dobrosavljević, Quantum Ripples in Strongly Correlated Metals, Phys. Rev. Lett. [**104**]{}, 236401 (2010).

J. W. Cahn, The kinetics of grain boundary nucleated reactions, Acta Metall. [**4**]{}, 449 (1956).

F. S. Khan and J. Q. Broughton, Simulation of silicon clusters and surfaces via tight-binding molecular dynamics, Phys. Rev. B [**39**]{}, 3688 (1989).

C. Z. Wang, C. T. Chan, and K. M. Ho, Empirical tight-binding force model for molecular-dynamics simulation of Si, Phys. Rev. B [**39**]{}, 8586 (1989).

S. Goedecker and L. Colombo, Efficient Linear Scaling Algorithm for Tight-Binding Molecular Dynamics, Phys. Rev. Lett. [**73**]{}, 122 (1994).

P. Koskinen and V. Mäkinen, Density-functional tight-binding for beginners, Comput. Mater. Sci. [**47**]{}, 237 (2009).

T. Ayral, T.-H. Lee, and G. Kotliar, Dynamical mean-field theory, density-matrix embedding theory, and rotationally invariant slave bosons: A unified perspective, Phys. Rev. B [**96**]{}, 235139 (2017).

T.-H. Lee, T. Ayral, Y.-X. Yao, N. Lanata, and G. Kotliar, Rotationally invariant slave-boson and density matrix embedding theory: Unified framework and comparative study on the one-dimensional and two-dimensional Hubbard model, Phys. Rev. B [**99**]{}, 115129 (2019).

M. Potthoff and W. Nolting, Surface metal-insulator transition in the Hubbard model, Phys. Rev. B [**59**]{}, 2549 (1999).

J. K. Freericks, Dynamical mean-field theory for strongly correlated inhomogeneous multilayered nanostructures, Phys. Rev. B [**70**]{}, 195342 (2004).

R. W. Helmes, T. A. Costi, and A. Rosch, Mott transition of Fermionic atoms in a three-dimensional optical trap, Phys. Rev. Lett. [**100**]{}, 056403 (2008).

M. E. Pezzoli, F. Becca, M. Fabrizio, and G. Santoro, Local moments and magnetic order in the two-dimensional Anderson-Mott transition, Phys. Rev. B [**79**]{}, 033111 (2009).

M. E. Pezzoli and F. Becca, Ground-state properties of the disordered Hubbard model in two dimensions, Phys. Rev. B [**81**]{}, 075106 (2010).

W. R. Magro, D. M. Ceperley, C. Pierleoni, and B. Bernu, [ Molecular Dissociation in Hot, Dense Hydrogen]{}, Phys. Rev. Lett. [**76**]{}, 1240 (1996).

K. T. Delaney, C. Pierleoni, and D. M. Ceperley, [ Quantum Monte Carlo Simulation of the High-Pressure Molecular-Atomic Crossover in Fluid Hydrogen]{}, Phys. Rev. Lett. [**97**]{}, 235702 (2006).

C. Attaccalite and S. Sorella, [ Stable Liquid Hydrogen at High Pressure by a Novel Ab Initio Molecular-Dynamics Calculation]{}, Phys. Rev. Lett. [**100**]{}, 114501 (2008).

G. Mazzola, A. Zen, and S. Sorella, Finite-temperature electronic simulations without the Born-Oppenheimer constraint, J. Chem. Phys. [**137**]{}, 134112 (2012).

J. Ma, P. Zhang, Y. Tan, A. Ghosh, and G.-W. Chern, Machine learning electron correlation in a disordered medium, Phys. Rev. B [**99**]{}, 085118 (2019).

H. Suwa, J. S. Smith, N. Lubbers, C. D. Batista, G.-W. Chern, and K. Barros, Machine learning for molecular dynamics with strongly correlated electrons, Phys. Rev. B [**99**]{}, 161107(R) (2019).
