Transient trapping into metastable states in systems with competing orders

TL;DR

This paper models the non-equilibrium dynamics of systems with competing orders after a rapid thermal quench, revealing how metastable states can form due to the dominance of faster-relaxing order parameters.

Contribution

It introduces a Ginzburg-Landau based relaxational dynamics framework to explain transient metastable states in systems with competing orders following a thermal quench.

Findings

Metastable states often form from the order parameter with faster relaxation.

Thermal fluctuations exponentially grow and influence the evolution.

The model explains laser-induced metastable states observed experimentally.

Abstract

The quench dynamics of a system involving two competing orders is investigated using a Ginzburg-Landau theory with relaxational dynamics. We consider the scenario where a pump rapidly heats the system to a high temperature, after which the system cools down to its equilibrium temperature. We study the evolution of the order parameter amplitude and fluctuations in the resulting time dependent free energy landscape. Exponentially growing thermal fluctuations dominate the dynamics. The system typically evolves into the phase associated with the faster-relaxing order parameter, even if it is not the global free energy minimum. This theory offers a natural explanation for the widespread experimental observation that metastable states may be induced by laser induced collapse of a dominant equilibrium order parameter.

Tables

Table 1. Table 1: Physical meaning of various time scales and Δ Δ \Delta .

Symbol	Physical meaning
$-t_{pump}$	When pump pulse arrives
$t_0$	When the temperature crosses $T_c$
$t_1,t_2$	$t_0$ for order I, II
$t_m$	Cooling (electron-phonon thermalization) time scale
$t_c$	The time that fluctuation becomes comparable to $1$
$t_{mu}$	The cooling time scale that makes $t_c=t_m$ (Eq. (40))
$t_{ms}$	See discussion around Eq. (47)
$t_{mc}$	The cooling time scale where mean field theory fails
$\Delta$	The critical ratio ${\gamma }_2{\alpha }_{L2}/\left({\gamma }_1{\alpha }_{L1}\right)$
$\delta$	$\left(1/\Delta -1\right)/2$

Equations

$$\frac{1}{{\gamma }_i}{\partial }_t{\psi }_i(\mathbf{r},t)=-\frac{1}{E_c}\frac{\delta F(t)}{\delta {\psi }_i(\mathbf{r},t)}+{\eta }_i(\mathbf{r},t).$$

$$F[{\psi }_1,{\psi }_2]=E_c\int d^D\mathbf{r}\left(f_1+f_2+f_c\right),$$

$$f_i=-{\alpha }_i{\psi }_i^2+({\xi }_{i0}\nabla {\psi }_i{)}^2+{\psi }_i^4$$

$$f_c=c{\psi }_1^2{\psi }_2^2$$

$$G_i(T)=\frac{T}{E_c{\xi }_{0i}^D}.$$

$${\rho }_k=\frac{1}{\sqrt{2\pi D_k(t)}}{e}^{-\frac{{\psi }_k^2}{2D_k(t)}}$$

$$⟨{\eta }_i(\mathbf{r},t){\eta }_i({\mathbf{r}}^{\prime },t^{\prime })⟩=\frac{2T(t)}{{\gamma }_i{E}_c}\delta (\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime })$$

$$\psi (\mathbf{r}{)}^{init}=\bar{\psi }+\sum _k{e}^{i\mathbf{k}\cdot \mathbf{r}}\delta {\psi }_k\equiv \bar{\psi }+V\int \frac{d^{D}k}{(2\pi {)}^D}{e}^{i\mathbf{k}\cdot \mathbf{r}}\delta {\psi }_k$$

$$\frac{1}{\gamma }{\partial }_t{\psi }_k=2{\alpha }_k(t){\psi }_k+{\eta }_k$$

$${\alpha }_k(t)=\alpha (t)-{\xi }_0^2{k}^2<0.$$

$${\psi }_k(t)={\psi }_k^{init}{e}^{S_k\left(t,-t_{pump}\right)}+\gamma {\int }_{-t_{pump}}^{t}d{t}^{\prime }{\eta }_k(t^{\prime })e^{S_k\left(t,t^{\prime }\right)}$$

$$S_k(t_a,t_b)=2\gamma {\int }_{t_b}^{t_a}dt{\alpha }_k(t)$$

$$D_k(t)=D_k^{init}{e}^{2S_k(t,-t_{pump})}+\frac{2\gamma }{E_{c}V}{\int }_{-t_{pump}}^{t}d{t}^{\prime }T(t^{\prime })e^{2S_k(t,t^{\prime })}$$

$$\bar{\psi }(t_0)=e^{-|{\alpha }_H|\gamma t_0}{e}^{-2|{\alpha }_H|\gamma t_{pump}}\bar{\psi }$$

$$D_k(t)=$$

$$\left(D_k(t_0)+\frac{2\gamma }{E_{c}V}{\int }_{t_0}^{t}d{t}^{\prime }{e}^{-2S_k(t^{\prime },t_0)}T(t^{\prime })\right)$$

$$⟨{\psi }_i(r=0{)}^2{⟩}_{t=t_c}\sim {\alpha }_i/2.$$

$${\xi }^2(t)={\xi }_0^2\left(8\gamma (t-t_0)+1/a\right)$$

$${⟨\psi (0)\psi (r)⟩}_t$$

$$=\frac{G/a}{{\left(16\pi \gamma \left(t-t_0\right)\right)}^{\frac{D}{2}}}{e}^{2S_0(t,t_0)}{e}^{-\frac{r^2}{2\xi (t{)}^2}}.$$

$$2S_0(t,t_0)=4{\alpha }_L\gamma t$$

$$4{\alpha }_L\gamma t_c=\ln \frac{1}{\zeta }+\frac{D}{2}\ln \left(4{\alpha }_L\gamma t_c\right)$$

$$\zeta =2(4\pi {)}^{-\frac{D}{2}}{\alpha }_L^{D/2-2}G$$

$$\frac{⟨{\psi }_2^2⟩}{⟨{\psi }_1^2⟩}=\frac{{\alpha }_{1L}}{{\alpha }_{2L}}{\left(\frac{{\gamma }_1}{{\gamma }_2}\right)}^{D/2}\frac{G_2}{G_1}{\left(\frac{1}{{\zeta }_1}{\left(\ln \frac{1}{{\zeta }_1}\right)}^{D/2}\right)}^{\frac{{\alpha }_{2L}{\gamma }_2}{{\alpha }_{1L}{\gamma }_1}-1}$$

$$\rho \left({\psi }_1(r),{\psi }_2(r)\right)=\frac{\mathrm{Exp}\left[-\frac{{\psi }_1^2}{2⟨{\psi }_1^2⟩}-\frac{{\psi }_2^2}{2⟨{\psi }_2^2⟩}\right]}{2\pi \sqrt{⟨{\psi }_1^2⟩⟨{\psi }_2^2⟩}}.$$

$${\gamma }_1^{-1}{\partial }_t{\psi }_1$$

$${\gamma }_2^{-1}{\partial }_t{\psi }_2$$

$${\psi }_1=\lambda {\psi }_2^{1/\Delta }$$

$$p_2$$

$$=\frac{1}{\pi }{2}^{\frac{1}{2\Delta }+\frac{1}{2}}\Gamma \left(\frac{1}{2}\left(1+\frac{1}{\Delta }\right)\right)\lambda ⟨{\psi }_1^2{⟩}^{\frac{1/\Delta -1}{2}}{\left(\frac{⟨{\psi }_2^2⟩}{⟨{\psi }_1^2⟩}\right)}^{\frac{1}{2\Delta }}$$

$$=\vartheta {\left(\ln \frac{1}{{\zeta }_1}\right)}^{-\frac{D}{2}\delta }{\zeta }_1^{\delta }$$

$$t_{life}\sim {\left(\frac{1}{p_2}\right)}^{\frac{1}{D}}\xi (t_c)/v\sim \frac{1}{\gamma }{\left(\ln \frac{1}{\zeta }\right)}^{\frac{1}{2}(1+\delta )}{\left(\frac{1}{\zeta }\right)}^{\frac{\delta }{D}}.$$

$$p_2\sim {\left(\ln \frac{1}{\zeta }\right)}^{-D\delta }{\zeta }^{2\delta }$$

Full text

I Introduction

Dynamical phase transitions Hohenberg and Halperin (1977); Polkovnikov et al. (2011); Bray (1994), in which systems are tuned through a phase transition by time variation of system parameters, are a fundamental topic of longstanding interest in many areas of science. For example, it is believed that cosmological expansion tuned the universe through the electroweak symmetry breaking transition Kibble (1976). Supercooled liquids are a widely studied terrestrial example. Spinodal decomposition Langer et al. (1975); Binder (1987); Carter and Johnson (1998); Alert et al. (2016) and Kibble-Zurek (KZ) Zurek (1996, 1985); Biroli et al. (2010); Chandran et al. (2012) theories have addressed important aspects of dynamical phase transition physics for systems characterized by an order parameter which is tuned through a first or second order transition respectively.

Systems with multiple competing or intertwined orders are of great current interest in condensed matter physics Grüner (1988); Fradkin et al. (2015). Examples include high Tc cuprates and transition metal dicalcogenides in which superconductivity and spin and/or charge density wave order compete and coexist as well as ‘colossal’ magnetoresistance manganites where ferromagnetic metal and charge ordered antiferromagnetic insulating states compete at low temperatures Tokura (2006); Zhang et al. (2016). Recent developments in “ultrafast” experimental technique Basov et al. (2011); Fausti et al. (2011); Nicoletti et al. (2014); Zhang et al. (2018a, b); Nicoletti et al. (2018); Cremin et al. (2019); Niwa et al. (2019); Zong et al. (2019); Suzuki et al. (2019); Rini et al. (2007) have made it possible to dynamically suppress one or more order parameters and study the subsequent evolution, raising the possibility of “steering” the order parameters into a desired metastable state.

The purpose of this paper is to provide theoretical insight into dynamical phase transitions in systems with multiple order parameters and in particular to draw attention to the crucial importance of the relative magnitudes of order parameter relaxation rates. We consider systems in which the relevant degrees of freedom are space-time dependent order parameter fields $\psi_{i}(\mathbf{r},t)$ defined from a fundamental theory by integrating out microscopic degrees of freedom such as electrons. For notational simplicity we deal here with a system with two real order parameters. Adding more real order parameters or making the order parameters complex does not alter our conclusions.

We consider two broad classes of behavior described by the equilibrium free energy landscapes sketched in Fig. 1 (a) strictly competing orders, where the free energy has two local minima, such that in each minimum only one of the two order parameters is nonzero, and Fig. 6 *intertwined order *parameters, where at the global minimum both order parameters are nonzero, but a metastable minimum exists in which only one of the order parameters is nonzero. When such systems are exposed to an experimentally relevant pump pulse they may be driven to a point in phase space near the origin ($o$) as indicated by the solid line trajectory in Fig. 1(a). We show that after the pump is turned off, the time evolution is dominated by the exponential amplification of very long wavelength fluctuations of the order parameters, and that even a modest difference in relaxation rates will drive the system to the minimum related to the faster evolving order, even if it is not the global free energy minimum. The probability of trapping into this metastable state is close to unity. The probability of going back to the equilibrium order scales as $p_{2}\sim\zeta^{\delta}$ where $\zeta\ll 1$ is the same Ginzburg parameter that controls the validity of static mean field theory and $\delta$ is a positive number defined latter. Nucleation dynamics operating on much longer time scales will lead eventually to relaxation to the global minimum Binder and Stauffer (1976); Hohenberg and Halperin (1977), but this physics is not explicitly considered here.

Transient dynamics in systems with a single order parameter have been extensively discussed; the literature is too large to review here but we note that a recent study of a quench to the superconducting state has derived from microscopics the model A dynamics used here Lemonik and Mitra (2017). Transient dynamics in systems with competing orders has been previously discussed in terms of deterministic dynamics of spatially uniform order parameters Kung et al. (2013); Ross Tagaras et al. (2019); Dolgirev et al. (2020). Our paper goes beyond the previous work by studying the formation and growth of spatial fluctuations and focusing on the difference in order parameter time constants.

Section II describes the physical picture and defines the formalism. Section C has the general solution to the dynamical problem. Section IV discusses the fast cooling limit which illustrates the essential physics. Section V analyzes the effect of slower cooling rate. Section VI discusses systems with intertwined orders. Section VII contains detailed comparison to recent experiments. Section VIII is a summary and conclusion, containing a discussion of the assumptions made and consequences of relaxing them. Appendices give detailed derivations of some of the formulas in the main text.

II Physical picture and formalism

We consider a system in which the important degrees of freedom are space-time dependent order parameter fields $\psi_{i}(\mathbf{r},t)$ obtained from a fundamental theory by integrating out quasiparticles. We assume that the order parameter fields evolve according to dissipative [relaxational Time-Dependent-Ginzburg-Landau (TDGL) or “Model A”] dynamics Gor’kov and Eliashberg (1968); Cyrot (1973); Hohenberg and Halperin (1977); Lemonik and Mitra (2017) defined by a free energy functional $F$ which is time dependent because of the applied pump field:

[TABLE]

Here $\gamma_{i}$ are the corresponding relaxation rates, $E_{c}$ is the condensation energy density and $\eta$ is a noise field determined by the microscopic degrees of freedom that were integrated out to obtain the order parameter theory. The free energy functionals are assumed to be of the general form

[TABLE]

as sketched in Fig. 1(a). Here the $\xi_{i0}$ are the bare coherence lengths, $D$ is the spatial dimension and

[TABLE]

describes the interaction between the two order parameters in the competing order case.

In our convention, $\psi_{i}$, $\alpha_{i}$, $c$ and $f_{i}$ are dimensionless, intensive and defined such that the quartic term in the free energy has coefficient $1$ and the $\alpha_{i}$ are of the order of unity at zero temperature. For cooperation ($c<0$) or weak competition ($0<c<2$), $F$ has a single minimum. For $c>2$, $F$ has two locally stable minima if $\alpha_{1}$ and $\alpha_{2}$ $>0$ and $\frac{2}{c}<\frac{\alpha_{1}}{\alpha_{2}}<\frac{c}{2}$ (see Appendix A). We study the $c>2$ (multiple minima) case in this paper. Following usual practice we assume that all parameters are temperature independent except the $\alpha_{i}=\kappa_{i}(T_{ci}-T)/T_{ci}$, which are positive at low temperatures, negative at high temperatures, vary smoothly with temperature and vanish at the respective critical temperatures $T=T_{ci}$ (we assume linear temperature dependence for simplicity). The nonequilibrium enters the formalism as a time dependence of the $\alpha_{i}$, determined by the time dependence of the effective temperature $T(t)$ of the microscopic degrees of freedom that were integrated out. We focus on the case $\alpha_{2}>\alpha_{1}$ but $\gamma_{2}\alpha_{2}<\gamma_{1}\alpha_{1}$ so minimum II is the equilibrium free energy minimum but the dynamics associated with minimum I is faster.

The condensation energy density $E_{c}$, in combination with $\xi_{i0}$, sets the relevant microscopic scales. An important dimensionless measure of the thermal fluctuations is

[TABLE]

The Ginzburg parameter defined in the conventional theory of critical phenomena is $G(T_{c})\alpha^{\frac{D-4}{2}}$ with mean field theory applying when the parameter is much less than unity. For example, in the weak coupling case of conventional superconductors, $G\sim\left(\mathrm{gap}/\mathrm{fermi\,\,energy}\right)^{D-1}$. The treatment that follows is formally valid in the $G\ll 1$ limit.

The stochastic Eq. (1) may be recast as a Fokker-Planck equation for a probability functional $\rho[\{\psi_{k}\}]$ that gives the distribution of fluctuations around the mean field value (see, e.g., Kramers (1940); Hohenberg and Halperin (1977); Risken (1996) and Appendix B). In the linearized approximation used below the probability functional is a direct product $\rho[\psi]=\prod_{k}\rho_{k}(\psi_{k})$ where

[TABLE]

is a Gaussian distribution for each Fourier mode $\psi_{k}$ of the field with time-dependent variance $D_{k}(t)=\langle\psi_{k}(t)\psi_{-k}(t)\rangle$ which we calculate below.

We are primarily interested in understanding experiments in which a system is highly excited by a pump pulse and the subsequent evolution is probed. We assume that the pump does not couple directly to the order parameters; rather, it excites microscopic degrees of freedom (e.g. electron quasiparticles or phonons) which thermalize very quickly (relative to the order parameter timescales) to a quasiequilibrium state described by an effective temperature $T(t)$ He and Millis (2016). The effective temperature is maintained by the pump at a high value $T_{H}$ for some time. After the pump is turned off, $T(t)$ evolves over a timescale $t_{m}$ to the true thermal equilibrium temperature $T_{L}$ determined by the bath. In most experiments, the bath is the lattice and $t_{m}$ is just the electron-phonon thermalization time scale which is typically of the order of picoseconds. This assumption is in essence the two temperature model of Rothwarf and Taylor Rothwarf and Taylor (1967) in which one set of degrees of freedom (e.g., electrons or a particular set of phonon modes) is excited to a high temperature and then relaxes back to the equilibrium temperature set by the rest of the system.

Within these assumptions, the instantaneous value of $T(t)$ determines the parameters $\alpha_{i}$ of the free energy and the noise. We take the noise correlators

[TABLE]

to be local in space and time and consistent with the fluctuation-dissipation theorem. Here the Boltzmann constant $k_{B}$ is set to unity and the average is over a probability distribution of the noise field. For the slow dynamics we consider here, the relevant frequency/momentum scale is well below those of the microscopic degrees of freedom that are integrated out, justifying the locality assumption of Eq. (6).

Representative time histories of $\alpha(t)$ are shown in Fig. 2. There are three time regimes: (a)pump on, covering the time interval $-t_{pump}<t<0$ in which the temperature $T=T_{H}$ and corresondingly the quadratic coefficient $\alpha=\alpha_{H}<0$; (b)relaxation, time $0<t<t_{m}$, during which $T$ evolves from $T_{H}$ through $T_{c}$ to $T_{L}$ while $\alpha$ evolves from $\alpha_{H}<0$ through $\alpha=0$ to $\alpha_{L}>0$; and (c) evolution, time $t>t_{m}$, $T=T_{L}$ and $\alpha=\alpha_{L}$. It is convenient also to introduce the time $t_{0}=\frac{T_{H}-T_{c}}{T_{H}-T_{L}}t_{m}$, at which $T=T_{c}$ and $\alpha(t)=0$ (the $t_{0}$ times for order parameters $\psi_{1}$ and $\psi_{2}$ are labeled as $t_{1}$ and $t_{2}$ in Fig. 2).

The physical picture is that for $t\in(-t_{pump},\,t_{0})$, the free energy landscape has its only minimum at $\psi=0$ (point $O$ in Fig. 1(a)). Thus the high temperature produced by the pump suppresses the mean field order parameter to nearly zero and the fluctuations remain bounded. After time $t_{0}$, point $O$ becomes unstable, the long wavelength fluctuations start to grow exponentially with time and the correlation length grows as $\xi\sim\xi_{0}\sqrt{8\gamma t}$. This process leads to the creation of large domains where most domains are in phase I as shown in Fig. 1(b).

III Dynamics

In this section we present a general solution of Eq. (1) with time dependent $\alpha$ and noise correlators as described above. The initial condition is

[TABLE]

where $\bar{\psi}$ is the initial order parameter, $\delta\psi_{k}$ represents thermal fluctuations about the initial ordered state and $V$ is the system volume. For simplicity of notation, here and henceforth we suppress the index $i$ labelling the different order parameters wherever possible. In the initial state, fluctuations are assumed small: $\bar{\psi}^{2}\gg\left<\delta\psi(r=0)^{2}\right>\equiv\sum_{k}\left<\delta\psi_{k}\delta\psi_{-k}\right>\sim G$. The pump acts to decrease $\bar{\psi}^{2}$ and increase the fluctuations. If $\bar{\psi}^{2}$ remains large compared to the mean square fluctuation amplitude, the state of the system is determined by a straightforward deterministic dynamics. This case is discussed briefly below, but our main interest is in situations in which the pump drives the initial order parameter to a value smaller than the root mean square fluctuation amplitude and the physics is determined by the evolution of the fluctuations.

Because Eq. (1) is first order in time it has no “memory”, so the evolution over one time regime fixes initial conditions for the next one. We will first consider the evolution over the decaying order parameter regime $t<t_{0}$; the resulting state of the system at $t_{0}$ is then the initial condition for the subsequent evolution.

In the pump-on regime the system is hot (temperature $T=T_{H}$) so the free energy is dominated by a large quadratic term which justifies the use of linearized dynamics even when the order parameter is not small. This means that in the pump-on regime we may study

[TABLE]

where

[TABLE]

We assume that the dynamics in the pump-on regime drives the order parameter to a small enough value that we may continue to use the linearized approximation throughout the $t<t_{0}$ regime and for some time into the growing fluctuations ($t>t_{0}$) regime. Conditions for the validity of this approximation will be presented below. The solution of Eq. (8) may be written as

[TABLE]

where the first term gives the propagation forward in time of the initial order parameter $\psi_{k}^{init}=\bar{\psi}\delta_{k,0}+\delta\psi_{k}$ with mean field part $\bar{\psi}$ and small fluctuations $\delta\psi$. The second term represents the propagation forward of fluctuations created by the noise after $-t_{pump}$. The accumulated phase $S$ is defined as

[TABLE]

and $S_{0}(t_{a},t_{b})$ has the interpretation as the signed area enclose by the solid lines and the time axis in Fig. 2. With Eq. (6), the square of Eq. (10) yields the fluctuation amplitude

[TABLE]

defined in Eq. (5) where $D_{k}^{init}$ is the initial fluctuation amplitude.

In the linear cooling profile approximation the initial mean field order parameter amplitude evolves as

[TABLE]

where the factor $e^{-2|\alpha_{H}|\gamma t_{pump}}$ gives the exponential suppression during the pump-on stage and the $e^{-|\alpha_{H}|\gamma t_{0}}$ factor gives the additional suppression during $t\in(0,\,t_{0})$. We assume that $|\alpha_{H}|\gamma(2t_{pump}+t_{0})$ is large enough that the mean field order parameter is reduced to a very small value at $t=t_{0}$, less than the mean square fluctuations. As shown in detail in section V, this assumption plus a small Ginzburg parameter implies that at time $t_{0}$, the system is well prepared in a disordered state with negligible mean field order parameter and small fluctuations.

We now consider the evolution of the distribution at times after $t_{0}$, where the system has cooled below the transition temperature ($\alpha(t)>0$) so long wavelength modes with $k<\sqrt{\alpha}\xi_{0}^{-1}$ grow exponentially with time. As long as the fluctuation does not become too large the linearized equation may be used for the dynamics so

[TABLE]

where as before the first term represents the propagation forward in time of the fluctuations existing at $t_{0}$ while the second term represents the additional contributions generated by the noise thereafter. A detailed analysis given in Appendix C shows that the second term in Eq. (14) is of the same order as the first in the situations of interest here.

The key observation is that long wavelength modes with $4\xi_{0}^{2}k^{2}\gamma(t-t_{0})<2S_{0}(t,t_{0})$ are exponentially amplified and we are interested in long times $\gamma(t-t_{0})\gg 1$ for which the growth is substantial ($e^{2S_{0}(t,t_{0})}\sim G^{-1}$). The exponential growth continues until the local mean square fluctuation amplitude of one of the order parameters becomes large enough that the nonlinearity becomes important to the dynamics, i.e., until $t$ reaches the crossover time $t_{c}$ defined by

[TABLE]

To compute $\langle\psi_{i}(0)^{2}\rangle$ we observe that at long times the important momentum dependence is controlled by the $e^{-4k^{2}\xi_{0}^{2}\gamma\left(t-t_{0}\right)}$ factor. Defining the time-dependent correlation length by

[TABLE]

we perform the momentum integral to obtain (up to an unimportant overall factor) the correlation function in real space

[TABLE]

The details of the pump and initial cooling enter Eq. (16) and (17) via the $a$ which varies from $\sim\alpha_{L}$ in the fast cooling limit to $\alpha_{KZ}$ in the slow cooling case (to be discussed in Sec. V). This variation leads to corrections that are subleading relative to the terms we consider and we will not explicity notate this dependence henceforth.

Eqs. (17) is the first key result. The exponential factors show that the fluctuations grow exponentially in time, and the spatial correlation is governed by the universal correlation length growth law $\xi\sim\xi_{0}\sqrt{8\gamma t}$ which does not depend on the equilibrium value of $\alpha$ or temperature-time profile.

IV Fast cooling limit

In this section we consider the fast cooling limit $t_{m}\rightarrow 0$ which illustrates the essential physics with minimal complexity. In the fast cooling limit the phase for the exponential amplification term in Eq. (17) is simply

[TABLE]

and Eq. (15) for the crossover time $t_{c}$ is

[TABLE]

where we have set $a$ in Eq. (17) to $\alpha_{L}$ without altering the leading behavior, and have defined

[TABLE]

which is in effect the usual Ginzburg parameter of the theory of critical phenomena. Similar results were previously presented by Lemonik and Mitra Lemonik and Mitra (2018) who noted the importance of the Ginzburg parameter in setting post-quench timescales. We see that $t_{c}$ is logarithmically large if $\zeta$ is small, i.e. if mean field theory works well for equilibrium. At time $t_{c}$ of the fast order I, the fluctuation of the slow order II is smaller by the ratio

[TABLE]

which is much less than unity if $\gamma_{2}\alpha_{2L}<\gamma_{1}\alpha_{1L}$ and $\zeta\ll 1$. At this stage of the evolution the two order parameters are independent and the joint distribution of local amplitudes at a position $r$ is the product of Gaussians:

[TABLE]

For $\gamma_{2}\alpha_{L2}<\gamma_{1}\alpha_{L1}$ the mean square values are very different, leading to the highly anisotropic joint distribution function shown in Fig. 3(a)(c). The probability distribution describing the space dependence of the fluctuations is derived in Appendix H and is plotted in Fig. 3(b)(d) for order I. We see that the fluctuations of order parameter I are highly correlated over scales out to $\xi(t)\gg\xi_{0}$ (the fluctuations of the slower order parameter II are correlated over slightly shorter distances).

Thus the physical picture at $t=t_{c}$ is of order parameter domains of typical size $\xi(t_{c})\approx\sqrt{\ln\frac{1}{\zeta}}\sqrt{\frac{2}{\alpha_{L}}}\xi_{0}\gg\xi_{0}$ within which the order parameters are nearly uniform and normally distributed and with the typical value of $\psi_{1}$ much larger than $\psi_{2}$, as illustrated in Fig. 1(b). Moreover, the typical value of local $\psi$ is now much larger than the new fluctuation scale induced by the noise, so that in the subsequent evolution the mean field dynamics dominates over the effect of the noise. Therefore, to study the subsequent evolution it suffices to consider the evolution within a domain, which is described by the uniform, deterministic TDGL equations, written here for the competing orders case:

[TABLE]

with initial conditions chosen from the joint probability distribution $\rho\left(\psi_{1}(0),\psi_{2}(0)\right)$. The issues associated with matching the solutions at the domain walls are a coarsening problem discussed briefly below.

The flow defined by Eq. (23) has a simple phase space structure with stable fixed points defined by the minima of $F$, as shown in Fig. 4. Each initial condition defines a trajectory that flows into one of the minima. For the physically relevant case $\alpha_{1L},\,\alpha_{2L}>0$ with $\frac{2}{c}<\frac{\alpha_{1L}}{\alpha_{2L}}<\frac{c}{2}$ there are four fixed points, with basins of attraction separated by a four-branched separatrix curve. We may estimate the position of the separatrix by matching the small $\psi$ regime, where the exponential growth requires

[TABLE]

to the requirement that the separatrix goes through the saddle point $(\psi_{1}^{2},\psi_{2}^{2})=(c\alpha_{2L}-2\alpha_{1L},\,c\alpha_{1L}-2\alpha_{2L})/(c^{2}-4)$. Here $\Delta=\gamma_{2}\alpha_{2L}/\gamma_{1}\alpha_{1L}<1$ and the coefficient is fixed by the matching condition.

By finding the relative weights of the probability distribution in the different basins of attraction we can estimate the relative volume fractions of the different order parameter domains. The volume fraction $p_{2}$ of order II domains is just the probability of $(\psi_{1},\,\psi_{2})$ lying in the basin of attraction of minimum II in Fig. 4, which at time $t=t_{c}$ can be estimated as

[TABLE]

where $\delta=(1/\Delta-1)/2>0$ and $\vartheta\sim 1$ can be found in Appendix I. Thus the proportion of $\psi_{2}$ domains is suppressed by a power law of the Ginzburg parameter $\zeta\ll 1$ and is negligibly small even if the time scales are just slightly different. Eq. (25) is the second key result of this paper.

*Life time of the metastable state—*Each domain then evolves to the appropriate minimum; the evolution takes a time of the order of $\frac{1}{\alpha_{L}\gamma}\ln\frac{1}{\zeta}$, after which the physical picture is of a set of domains, most of which have $\psi_{1}=\pm\sqrt{\alpha_{1L}/2}$ and $\psi_{2}\approx 0$ (i.e. are in phase I) while a small volume fraction of the sample are phase II domains where $\psi_{2}=\pm\sqrt{\alpha_{2L}/2}$ and $\psi_{1}\approx 0$, as illustrated in Fig. 1(b). However, long range order is not established. The subsequent evolution is determined by spontaneous nucleation of phase II regions in the dominant phase I domains, and by growth of the existing $\psi_{2}$ domains. The timescale for ultimate equilibration thus depends both on nucleation rates and on domain wall dynamics, both of which are beyond the scope of this paper. We do, however, provide a likely lower limit on the equilibration time by considering the free growth of $\psi_{2}$ domains, assuming no domain wall pinning and nucleation which is an exponentially slow process. The speed of domain wall motion is at the order of $v\sim\gamma\xi_{0}$ as long as the free energy difference $\delta f$ between the two minima is order one. Assuming the phase II domains are evenly distributed among the phase I domains we can estimate the equilibration time as

[TABLE]

If the order parameters are complex ones corresponding to $U(1)$ symmetry breaking, as in the case of superconductivity (SC) being order I and charge density wave (CDW) being order II, one just needs to double the degrees of freedom for each order before nonlinear dynamics is onset. As a result, Eq. (25) gives the offset of trapping probability as

[TABLE]

and Eq. (26) is modified to

[TABLE]

After $t_{c}$, the different SC regions are characterized by different order parameter phases which can continuously synchronize to leave behind topological vortices. The number density of the vortices scales as $n\sim 1/\xi(t_{c})^{2}\sim\left(\xi_{0}^{2}\frac{2}{\alpha_{L}}\ln\frac{1}{\zeta}\right)^{-1}$, different from the Kibble-Zurek scaling Zurek (1985) since the latter applies to the slow cooling limit.

V Finite cooling rate

We now ask how the physics is modified as the cooling time $t_{m}$ is increased from zero. The essential picture derived in the previous section of long length-scale domains of one or the other order parameter still applies, but because the time $t_{0}$ of transition from exponential decay to exponential growth is earlier for order II than that for order I, the $\psi_{2}$ fluctuations will have a longer period of growth than the $\psi_{1}$ fluctuations. The longer period of growth will compensate for the faster dynamics of $\psi_{1}$, meaning that the condition on the difference in relaxation rates required for the system to evolve to minimum I becomes more stringent. A second issue is that the cross over to nonlinear dynamics may occur at a time $t<t_{m}$ before the relaxation of the $\alpha$ to its equilibrium value is complete, meaning that the free energy landscape at the point of crossing to nonlinearity differs from the equilibrium one. For these reasons the behavior for given $t_{m}$ depends on the ratio of relaxation rates $\gamma_{2}/\gamma_{1}$ in a somewhat complicated manner. The various regimes are shown in Fig. 5.

V.1 The state at $t_{0}$

We first characterize the state at $t=t_{0}$. If the mean field order parameter is reduced to a small value, then the thermal fluctuations existing at time $t=-t_{pump}$ (the first term in Eq. (12)) will be reduced to a completely negligible level so the order parameter at $t=t_{0}$ and subsequent times is determined entirely by the random noise. We distinguish fast cooling ($|\alpha_{H}|\gamma t_{0}\leq 1$) and slow cooling ($|\alpha_{H}|\gamma t_{0}\geq 1$) regimes according to whether the cooling to $t_{0}$ after the pump is turned off has a significant effect on the order parameter. The linearized dynamics means that the corresponding distribution function is the product of Gaussians given in Eq. (5). Averaging the solution for $\psi_{k}(t_{0})$ over the noise using Eq. (6) shows that the fluctuation distribution half-width defined in Eq. (5) is

[TABLE]

The integral in Eq. (29) may easily be evaluated numerically, and in the linear quench approximation may be expressed exactly in terms of error functions (see Appendix D). Here we present results in important limits which explicate the basic physics. In the fast cooling limit the portion of the integral from $t=0\rightarrow t_{0}$ makes a negligible contribution and we find

[TABLE]

indicating that in the fast cooling limit the fluctuations at $t=t_{0}$ are those of the hot thermal state created by the pump, with distance $|\alpha_{H}|$ from criticality and correlation length $\xi_{H}=\xi_{0}/\sqrt{|\alpha_{H}|}$. In the slow cooling regime only times near $t_{0}$ are important and we find

[TABLE]

where the effective distance from criticality $\alpha_{KZ}$ and corresponding correlation length $\xi_{KZ}$ are given by

[TABLE]

which depend on the square root of the cooling rate, consistent with Kibble-Zurek scaling Zurek (1996, 1985) and the mean field exponents of the problem at hand.

The correlation function in real space is given by

[TABLE]

where the upper cutoff $k_{c}$ which we expect to be of the order of a few times $\xi_{0}^{-1}$ is required to make the integral finite as $r\rightarrow 0$. The momentum integral is dominated by large momenta for which $D_{k}\sim(k\xi_{0})^{-2}$, so the local fluctuation amplitude is

[TABLE]

where G (Eq. (4)) is to be evaluated at temperature $T=T_{H},T_{c}$ in the fast and slow quench limits respectively and $\xi$ takes the value appropriate for the relevant limit.

To summarize, if the initial pump and subsequent cooling are strong enough to drive the initial mean field order parameter to a level smaller than the local root mean square fluctuation then at time $t=t_{0}$ the order parameter fluctuations in both the rapid and slow quench cases are described by a Gaussian (mean field-like) probability characterized by a temperature ($T_{H}$ or $T_{c}$), a distance from criticality ($|\alpha_{H}|$ or $\alpha_{KZ}$) and the associated correlation length $\xi^{2}_{H,KZ}=\xi_{0}^{2}/\alpha_{H,KZ}$. The local mean square local fluctuation order parameters are of the order of $G$, justifying Eq. (17) even in the slow cooling case. Using Eqs. (13) and (36) we see the criterion for suppression of the mean field order parameter is roughly $|\alpha_{H}|\gamma\left(2t_{pump}+t_{0}\right)\gg\ln(1/G)$. The linearized analysis used here requires that $\left<\psi^{2}(r=0)\right>\ll\alpha$. A renormalization group-improved treatment will break down when $G\alpha(t)^{\frac{D-4}{2}}\sim 1$ (the Ginzburg criterion) which sets an upper limit on the cooling time $t_{m}\sim\zeta^{-1/(1-D/4)}/\gamma\equiv t_{mc}$ in the slow quench regime.

It is important to remember that the $T_{c}$ (and thus $t_{0}$) of the two different order parameters are different as are the bare correlation lengths and relaxation constants, so especially in the slow cooling limit the probability distribution functions of the two order parameters will differ, especially at long wavelengths, although Eq. (36) shows that the local fluctuation amplitudes, which are determined by short wavelength fluctuations, are not too different for the two order parameters.

V.2 The trapping condition

Now we analyze the effect of finite cooling rate on the condition for trapping into phase I. To simply the formulas, in the main text we focus on the exponential growth and neglect power-law prefactors, thus approximating $\langle\psi^{2}\rangle_{t}\sim Ge^{2S(t,t_{0})}$. Our main focus will be on establishing how small $\gamma_{2}$ must be relative to $\gamma_{1}$ for the system to evolve with high probability into the metastable minimum I. We will find dependence of the critical ratio $\Delta=\frac{\gamma_{2}\alpha_{L2}}{\gamma_{1}\alpha_{L1}}$ is a scaling function of the variable $t_{m}/t_{mu}$, where $t_{mu}$ is the cooling time at which the onset of nonlinearity $t_{c}$ coincides with the equilibration time $t_{m}$.

To begin the analysis we note that for $t_{m}>0$ and $t>t_{m}$ the accumulated phase becomes (after eliminating $t_{0}$ in favor of $\alpha_{L}$)

[TABLE]

and Eq. (19) for the crossover time becomes

[TABLE]

while Eq. (21) becomes

[TABLE]

The factor $\left(1/\zeta_{1}\right)^{\frac{\alpha_{2L}\gamma_{2}}{\alpha_{1L}\gamma_{1}}-1}$ is Eq. (21) with only the leading term in $t_{c}$ retained and the exponential factor expresses the additional growth of $\psi_{2}$ due $\alpha_{2}$ crossing zero earlier than $\alpha_{1}$. When

[TABLE]

we have $t_{c}=t_{m}$, i.e., the onset of nonlinearity occurs at $t=t_{m}$. Thus the onset of nonlinearity occurs before equilibration only for cooling rates very slow relative to the basic order parameter timescales by a factor of the order of the log of the Ginzburg parameter. Using

[TABLE]

we see that $\langle\psi_{2}^{2}\rangle/\langle\psi_{1}^{2}\rangle<1$ provided that $\Delta=\alpha_{2L}\gamma_{2}/(\alpha_{1L}\gamma_{1})$ is less than a critical value defined by

[TABLE]

as shown in Fig. 5. In the $t_{m}\rightarrow 0$ limit Eq. (42) reverts to the previous result $\alpha_{2}\gamma_{2}<\alpha_{1}\gamma_{1}$ (up to logarithmic corrections), but as $t_{m}$ increases the constraint on $\gamma_{2}$ becomes more stringent and when $t_{m}=t_{mu}$ Eq. (42) becomes

[TABLE]

For $t_{m}>t_{mu}$, $\psi_{1}$ reaches nonlinearity before $t_{m}$ which is the regime considered by Kibble-Zurek theory Zurek (1985). Note that relaxational dynamics predicts a logarithmic correction to Kibble-Zurek scaling, as described in Appendix F. In this time regime the accumulated phase may be written

[TABLE]

and after some algebra Eq. (15) for phase I can be written as

[TABLE]

which is obviously before $t_{m}$ if $t_{m}>t_{mu}$. The condition $S(t_{c},t_{2})<S(t_{c},t_{1})$ becomes

[TABLE]

which reduces to our previous Eq. (43) when $t_{m}=t_{mu}$ and drops as $1/t_{m}$ for large $t_{m}$ so that as the equilibration time becomes extremely long, the system would evolve to the equilibrium minimum unless the relxation rate $\gamma_{2}$ becomes exceptionally small.

Even if Eq. (46) is satisfied, the system will only evolve to the metastable minimum if $\alpha_{i}$ are such that the metastable minimum exists at the time order I crosses to nonlinearity, i.e., if $\alpha_{2}(t_{c})<\frac{c}{2}\alpha_{1}(t_{c})$ which yields

[TABLE]

If $t_{m}>t_{ms}$, as denoted by ‘?’ in Fig. 5, the following would happen in the cooling process: at $t_{c}$, the time for $\psi_{1}$ crossing over to nonlinear dynamics, the free energy landscape has not recovered enough such that order I is not yet a local minimum. Thus trapping into the order I state won’t necessarily happen even if $\langle\psi_{2}^{2}\rangle\ll\langle\psi_{1}^{2}\rangle$ at this time. We see that typically $t_{ms}$ cannot be too much larger than $t_{mu}$, unless either $T_{c2}-T_{c1}$ is very small or $\kappa_{1}\gg\kappa_{2}$ or $c\gg 2$. The various time scales are collected in Table 1.

Below the upper critical dimension $D=4$, our ‘time dependent fluctuation around mean field’ approach fails if at the predicted cross over time $t_{c}$, the system is still inside the critical regime $G\alpha(t)^{\frac{D-4}{2}}\sim 1$ (the Ginzburg criterion) where even renormalization group improved treatment breaks down. This imposes an ultimate upper limit

[TABLE]

for the cooling time, see Appendix F. This much larger time scale is deep inside the ‘?’ region and is also labeled in Fig. 5.

VI Intertwined order

To study the intertwined case we modify Eq. (2) to shift one of the minima away from one of the axes. One simple choice is to add a term to $f_{c}$ so

[TABLE]

and with, now, $0<c<2$ and $d_{1}>0$. We assume $T_{c2}>T_{c1}$ but that the difference in $T_{c}$ is not too large, and $d_{1}$ is not too small. In this case (see Appendix A for details), as temperature is lowered the system first enters a phase with only $\psi_{2}\neq 0$ and then at a lower temperature the phase with $\psi_{1}\neq 0,~{}\psi_{2}=0$ becomes locally stable although not the global minimum. At a still lower temperature the global minimum is gradually shifted to I+II where a nonzero $\psi_{1}$ component appears. If we identify $\psi_{2}$ with density wave order and $\psi_{1}$ with superconductivity, this scenario may describe stripe ordered cuprates (e.g., La2-xBaxCuO4 around $x=1/8$): the so called pair density wave (PDW) state Fradkin et al. (2015). The free energy analysis of the $\psi_{2}/\psi_{1}$ minimum has previously been discussed Fradkin et al. (2015); we have generalized the free energy so that it also includes a metastable phase with purely superconducting order and will argue that this generalization is needed to describe recent ultrafast experiments Cremin et al. (2019).

The considerations sketched in the previous sections carry over directly to the intertwined order case, as shown by the red line in Fig. 6. However, an additional interesting effect may occur if we relax the assumption that the pump heats up the bath appropriate to both order parameters. If the two orders couple to different microscopic degrees of freedom, then one may consider the case when only the free energy landscape for one order is changed. In particular, in the case of coupled superconducting and charge orders, one may imagine that the charge order couples to phonons much more strongly than do the electrons, so driving the phonons would affect the CDW much more strongly than the superconductivity. If the system starts in a minimum with both order parameters nonzero (as is the case for intertwined orders) and only $\alpha_{2}$ is driven to negative, leaving $\alpha_{1}$ positive, the transient free energy landscape will have only one minimum (I) and the mean-field dynamics will drive the system into it, as shown by the blue trajectory in Fig. 6. In this process, small fluctuations of the order parameter can be neglected and one can apply deterministic TDGL dynamics to the mean field order parameters. This mechanism does not require faster relaxation for $\psi_{1}$; all that is needed is that $\psi_{2}$ remains suppressed for long enough that the system evolves to the I minimum. This timescale is set by the time required for the order parameter to cross the basin boundary, $t_{s}\sim\frac{1}{\alpha\gamma}\ln\frac{\psi_{2m}}{\psi_{2b}}$ where $\psi_{2m}$ is the value of $\psi_{2}$ at the original point $\textit{{I}}+\textit{{II}}$ and $\psi_{2b}$ is its value at the intersection between the blue trajectory and the basin boundary. For shorter pump durations or for pumps that reduce $\alpha_{1}$ too much, the system would relax back to the global minimum as illustrated schematically by the green trajectory in Fig 6.

VII Experiments

Competing orders have been reported in many materials, and an increasing number of ultrafast experiments are appearing, including studies of competing charge density waves in tri-tellurides Kogar et al. (2020); Zong et al. (2019), ferromagnetic domain formation in charge ordered manganites Zhang et al. (2016), and charge and magnetic order in rare earth nickelates. Much attention has focused on reports of transient superconductivity appearing in materials that have low temperature nonsuperconducting density wave states but may reasonably be expected to have competing superconducting states Fausti et al. (2011); Nicoletti et al. (2014); Zhang et al. (2018a, b); Nicoletti et al. (2018); Cremin et al. (2019); Niwa et al. (2019); Suzuki et al. (2019). The appearance of long lived superconducting-like metastable states has been seen in cuprates Fausti et al. (2011); Nicoletti et al. (2014); Zhang et al. (2018a, b); Nicoletti et al. (2018); Cremin et al. (2019); Niwa et al. (2019) through optical pump-probe by several independent groups and in FeSe through time resolved APRES Suzuki et al. (2019). This widespread mystery has heretofore been theoretically addressed via explorations of models in which the nonequilibrium drive changes the microscopic Hamiltonian, creating new physics not existing in equilibrium Kennes et al. (2017); Babadi et al. (2017); Sentef et al. (2017); Chiriacò et al. (2018); Wang et al. (2018) and via TDGL analyses Kung et al. (2013); Ross Tagaras et al. (2019); Dolgirev et al. (2020) with deterministic uniform dynamics. We argue that in many cases, the physical picture developed here is more relevant.

VII.1 Cuprates

As an example, we consider the relevant parameters for the cuprate La1.675Eu0.2Sr0.125CuO4 (LESCO1/8), the first material found to exhibit such transient phenomenon. At $10\,\mathrm{K}$ where this compound is not superconducting due to competition of charge order, Fausti et al Fausti et al. (2011) pumped the system with a strong infrared pulse. A metastable state appeared within picoseconds and lived for at least nanoseconds. Most strikingly, it exhibits superconducting like terahertz optical response.

In discussing the application of our theory, the first issue is timescales. The timescales associated with gap recovery in cuprate superconductors are typically of the order of $\tau_{sc}\sim 1\,\mathrm{ps}$ Smallwood et al. (2012); similar timescales are reported in studies of transient enhancement of the photoresponse in LESCO1/8 Fausti et al. (2011) and Y-Bi2212 Giusti et al. (2019). Time resolved x-ray and electron diffraction experiments reported CDW relaxation timescales in the wide range from 4 to 1000 ps in Transition Metal Dichalcogenides (TMD) Eichberger et al. (2010); Erasmus et al. (2012) where the CDW order is coupled to the lattice. Timescales of only a few $ps$ were reported for the charge order in cuprates Hinton et al. (2013), but the time scale for the stripe order that strongly competes with the superconductivity is not know, and may be long because the stripe order couples strongly to the lattice Fausti et al. (2011). Assigning SC to order I and CDW to order II, we assume that $\tau_{sc}/\tau_{cdw}=1/3$ for LESCO1/8 and $\Delta=\gamma_{2}\alpha_{2L}/\gamma_{1}\alpha_{1L}\approx 1/3$ since the $\alpha_{iL}\sim 1$ for both CDW and SC.

The next issue is the Ginsburg parameter $G$. The coherence lengths are of the order of a few $nm$ and Gaussian fluctuations are observed for temperatures within $10\%$ of $T_{c}$, so $G$ is unlikely to be as small as it is in conventional materials. If there were no competition to CDW, the superconducting $T_{c}$ of LESCO1/8 is about $40\,\mathrm{K}$, indicating a zero temperature superconducting gap $\sim 14\,\mathrm{meV}$ which is perhaps a factor of $\sim 100$ less than the Fermi energy. We suggest that $\zeta\sim G\sim 10^{-2}$ (lower panel of Fig. 3) may be appropriate since the material is effectively two dimensional.

The experiment of Fausti et al Fausti et al. (2011) could then be interpreted as destruction of both orders followed by growth of fluctuations described in section IV. The pump along a-axis had a fluence of $1\,\mathrm{mJ/cm^{2}}$, the penetration depth was about $200\,\mathrm{nm}$ and almost all the photon energy was absorbed since the reflectivity at that frequency is nearly zero. Together with the electronic specific heat of $C_{el}/T\approx 3\,\mathrm{mJ\,k^{-2}mol^{-1}}$ Michon et al. (2019) and the lattice constant of $a=3.8\,\mathrm{\AA},\,c=13.2\,\mathrm{\AA}$ Radaelli et al. (1994), we estimate that the electronic system was transiently heated up to $T_{H}=2000\,\mathrm{K}$, much larger than the critical temperatures of both CDW ($80\,\mathrm{K}$) and SC. Thus it is reasonable to assume all orders are destroyed by the pump. The cooling time scale Smallwood et al. (2012); Fausti et al. (2011) should not be significantly larger than $\tau_{sc}$, thus the equations in the fast cooling limit should give reasonable estimations. Application of Eq. (27) in 2D yields $p_{2}\approx 5\times 10^{-6}$ as the volume fraction of CDW domains in the transient state, meaning most volume is transformed to the SC state. Eq. (28) predicts the lifetime of the metastable state to be about one nanosecond. These estimations qualitatively explain the phenomenon seen in LESCO1/8 but note that they depend sensitively on the values of $G$ and $\Delta$, see Fig. 7 and Eqs. (27) and (28).

The cuprate La1.885Ba0.115CuO4 has a density wave transition at $T_{co}=53\,\mathrm{K}$ followed by a second transition at $T_{c}=13\,\mathrm{K}$ to a state with both density wave order and weak superconductivity. Cremin et al Cremin et al. (2019) recently reported that upon moderate near infrared ($1.55\,\mathrm{eV}$) pump pulse (fluence of $0.1\,\mathrm{mJ/cm^{2}}$) along c-axis, the weak superconducting state may be converted to a long-lived strong superconducting like state within picoseconds. The key observation is that this state can be created only if the static system is in the weakly superconducting state below $T_{c}$ of bulk superconductivity. If the static temperature is even slightly above $T_{c}$, the strong superconducting like metastable state is not created. We interpret the result as suggesting that the equilibrium state is an intertwined state with both superconducting and charge order, and the pump couples more strongly to the charge order while it is not strong enough to kill both orders. The transient phenomenon is due to the second mechanism described in Section VI. The transition time is $t_{s}\sim\frac{1}{\alpha\gamma}\ln\frac{\psi_{2m}}{\psi_{2b}}\sim\,\mathrm{ps}$ if one uses $\gamma\sim 1\,\mathrm{ps^{-1}}$. As the static temperature gets close to $T_{c}$, the transition time diverges logarithmically since $\psi_{2b}$ approaches zero, explaining the key observation.

VII.2 Transition metal tri-tellurides

Competing phases occur in non-superconducting contexts. In equilibrium, LaTe3 and CeTe3 exhibit long ranged CDW order (denoted as c-CDW). However, when these systems are driven out of equilibrium by a sufficiently strong near infrared pump, a different CDW order (denoted as a-CDW), distinguished from the equilibrium one by the wavevector, appears Kogar et al. (2020); Zhou et al. (2019). This was interpreted as a-CDW states living in topological vortices created in the dominant cCDW state by the Kibble-Zurek mechanism. However, in the experiments by Kogar et al Kogar et al. (2020), stronger a-CDW signal was observed for a stronger pump while the Kibble-Zurek mechanism says the number of vortices depends only on the cooling time constant through $T_{c}$, a system parameter that does not quite depend on pump fluence. Our framework in section IV is an alternative explanation. Assume that the pump destroys the mean field order parameter of the c-CDW but not to zero, its recovery would suppress the growth of fluctuations of a-CDW. Therefore, a stronger pump would suppress c-CDW to a smaller value which gives more space for a-CDW fluctuations to grow. A quantitative application of our theory requires more information on the relaxation rates. Since both orders are CDWs in this case, their time scales should be comparable and we place LaTe3/CeTe3 close to $\Delta=1$ in Fig. 5.

VII.3 Manganites

Zhang and McLeod et al have reported that in charge ordered insulating films made of La2/3Ca1/3MnO3 Zhang et al. (2016); McLeod et al. (2019), exposure of pump radiation can create domains of ferromagnetic metallic order, which grow in size with successive pump pulses and at low temperature, do not revert to the ground state on measureable time scales. Analysis of these experiments requires extension of our theory to the case of first order free energy landscapes.

VIII Discussion

We presented a dynamical phase transition theory of a pumped system with competing orders, based on a Landau theory of non conserved order parameters with relaxational dynamics coupled to a quasi thermal bath. We focused on the case where an applied (“pump”) field changes the free energy landscape, thereby driving any mean field order parameters to vanishingly small values and studied in detail the growth of fluctuations after the pump is removed. We presented a general treatment valid for cooling rates that are fast or slow compared to the basic order parameter time scales, presented a scaling theory valid in the slow cooling limit, and connected our results to the Kibble-Zurek theory of systems quenched through a critical point.

We computed the probability distribution of order parameter fluctuations, identified the exponential growth of long wave length fluctuations characterized by a universal correlation length growth $\xi\sim\sqrt{8\gamma t}$, leading to a large domain structure. We showed that in physically reasonable cases, a modestly larger relaxation rate of the subdominant order can lead the system to a metastable state of domains, most of which are in this subdominant phase. We derived scaling functions for the volume fraction of different domains and the lifetime of the metastable state. Due to universality of the Landau theory, it applies to solid state systems with competing orders Fausti et al. (2011); Nicoletti et al. (2014); Zhang et al. (2018a, b); Nicoletti et al. (2018); Cremin et al. (2019); Niwa et al. (2019); Zong et al. (2019); Suzuki et al. (2019); Kogar et al. (2020); Giusti et al. (2019); Eichberger et al. (2010); Erasmus et al. (2012); Hinton et al. (2013), cold atoms Guardado-Sanchez et al. (2018) and even the early universe Kibble (1976).

The fluctuation theory developed here naturally explains the key features of the observations of metastable states Fausti et al. (2011); Suzuki et al. (2019): (a) the long life time ($>200\,\mathrm{ps}$) of the superconducting like state; (b) the finite frequency conductivity peaks Nicoletti et al. (2014); Hunt et al. (2015) maybe explained by the transiently created domain structure (Fig. 1(b)) which allows coupling of far field radiation to plasmonic modes. A further proof of consistency is that similar transient phenomenon are observed for both infrared and optical pumping, suggesting that the main effect of the pump is incoherent, related to heating of the microscopic degrees of freedom.

Next we discuss the assumptions underlying our approach. (a) We assume the existence of a well-defined, time-dependent temperature T(t) for the high energy electronic degrees of freedom throughout the full time evolution. Indeed, quasiparticle thermalization time ($\sim 10\,\mathrm{fs}$) is usually much shorter than the dynamics of the collective order parameter ($0.1\sim 1\,\mathrm{ps}$). See, e.g., Ref. He and Millis (2016). Therefore, for the slow dynamics of the order parameter we consider, it is legitimate to assume a well defined temperature for high energy degrees of freedom which acts as the bath for the order parameters. (b) Our theory relies on a small Ginzburg parameter $G$, which is the same parameter that controls the validity of mean field theory. Thus our theory applies wherever mean field does, e.g., for SC and CDW where $\mathrm{gap}/\mathrm{fermi\,\,energy}\ll 1$. Although the application of mean field arguments to strongly correlated systems such as cuprates maybe questionable, our estimation using $G\sim 10^{-2}$ still renders a trapping probability close to unity. (c) We worked in the context of Ginzburg-Landau theory with relaxational dynamics which holds close to $T_{c}$ Lemonik and Mitra (2017) or for superconductors rendered gapless by magnetic impurities Gor’kov and Eliashberg (1968). However, relaxational dynamics is not essential to our conclusions. If a second order time derivative is added to Eq. (1), it does not change the picture of exponential growth of the long wave length thermal fluctuations, as long as there is substantial damping to take away the energy. Extension of our approach to the under damped (Hamiltonian dynamics) case is of interest. (d) For the trapping into the metastable order I, we need it to relax faster than the equilibrium order II. This probably happens for competing SC and CDW orders since the latter often couples to the lattice which is heavy. Since the trapping probability crosses to unity exponentially as $\Delta$ crosses one, order I just needs to be moderately faster than order II, as shown by Fig. 7 where the crossing is quick already for a relatively large $G$. (e) We assumed the white noise in Eq. (6) to characterize thermal fluctuations. The underlying assumption is that there is a length scale separation between the long length scale physics of the order parameter dynamics and the presumably short length scale physics of the microscopic degrees of freedom that are integrated out to obtain the order parameter theory. If the microscopic degrees of freedom have a length scale comparable to order parameter length scales, then the partial differential equations we analyse should be replaced by integro-differential equations. Analysis of this situation is beyond the scope of our paper, but we believe that it would not change our main conclusions.

Our work defines directions for future research. On the theoretical side, detailed application of our theory to specific materials, extension to the case of conserved order parameters, to strongly interacting field theories ($G\sim 1$) and to the case of interfaces between domains with different orders Del Re et al. (2016) are all of interest. Similar conclusions are expected for quenching through a quantum critical point at zero temperature. Instead of thermal fluctuations, quantum fluctuations will be exponentially amplified. Also, generalization of our formalism to the case where the pump couples coherently to the order parameters, as would happen for Terahertz pumps is of interest. On the experimental side, the growing fluctuations and the induced domain structure (Fig. 1(b)) or the topological vortices can be ideally probed by time-space resolved techniques, e.g., ultra fast Terahertz near field microscopy or ultra fast scanning tunneling microscopy. Moreover, the growing SC fluctuation and thus superfluid density indicates increasing Drude weight in the non equilibrium optical conductivity, leading to novel effects on the collective modes Sun et al. (2016) and THz reflectivity.

Acknowledgements.

We acknowledge support from the Department of Energy under Grant DE-SC0018218. We thank M. M. Fogler, R. D. Averitt, D. N. Basov, M. Eckstein, W. Yang, D. Golez and Y. He for helpful discussions.

*Note added.—*A very recent paper Dolgirev et al. (2019) uses similar methods to study the post quench growth of order parameter phase fluctuations.

Appendix A Equilibrium Free Energy and Phase Diagram

A.1 Competing orders

The phase diagram is shown in Fig. 8. We write Eq. (2) for the spatially uniform case using Eq. (3) and writing $\psi_{1}^{2}=R^{2}cos^{2}\theta$, $\psi_{2}=R^{2}sin^{2}\theta$. We obtain

[TABLE]

Minimizing with respect to $cos2\theta$ gives

[TABLE]

so

[TABLE]

and minimizing over $R$ gives

[TABLE]

so

[TABLE]

and

[TABLE]

The alternative solution is to set one of the $\psi=0$, obtaining

[TABLE]

Thus we see that if $c>2$ then the mixed solution costs energy and lower energy solutions are $2\theta=0$ and $\pi$. Expanding around the $\theta=0$ solution we obtain

[TABLE]

or

[TABLE]

so we see that the minimum at $\psi_{2}=0$ is only stable if $\frac{c}{2}>\frac{\alpha_{1}}{\alpha_{2}}$; expanding around the other minimum changes the sign of the $\theta^{2}$ term and interchanges $\alpha_{2}$ and $\alpha_{1}$, justifying the inequalities presented in the main text.

A.2 Intertwined orders

We write Eq. (2) for the spatially uniform case using Eq. (3) and Eq. (49), now in their original form

[TABLE]

We suppose $0<c<2$ and $d_{1}>0$, assume $T_{c2}>T>T_{c1}$ and consider the physics as $T$ is decreased below $T_{c2}$. Initially we have a solution with $\psi_{2}^{2}=\frac{\alpha_{2}}{2}$ and $\psi_{1}=0$. As the temperature is decreased, $\alpha_{1}-c\alpha_{2}/2$ may become positive; if this occurs, a $\psi_{1}$ component is added to the solution with $\psi_{2}\neq 0$. This instability takes place at a temperature lower than $T_{c1}$ if $c>0$ and at a temperature higher than $T_{c1}$ if $c<0$. We interpret this mixed state as having intertwined order, since both order parameters are non-zero. As $T$ is decreased below $T_{c1}$ a second extremum (saddle point) appears at $\psi_{1}^{2}=\frac{\alpha_{1}}{2},~{}\psi_{2}=0$. After $\alpha_{1}$ becomes large enough such that $\frac{c\left|\alpha_{1}\right|}{2}+d\frac{\alpha_{1}^{2}}{4}>\alpha_{2}$, this saddle point becomes stable to variations in $\psi_{2}$ and thus a local minimum.

Appendix B The Fokker-Planck Equation and its Approximations

If one defines the probability functional $\rho[\psi]$, the stochastic TDGL equation (1) is equivalent to the Fokker-Plank equation Hohenberg and Halperin (1977) for the probability functional $\rho[\psi(\mathbf{r},t)]$:

[TABLE]

where $\partial_{\psi_{i}}$ should be understood as functional derivative: $\partial_{\psi_{i}}\equiv\frac{\partial}{\partial\psi_{i}(\mathbf{r},t)}$. The averages $\langle\psi^{2}\rangle$ taken throughout the paper is over the probability functional $\rho[\psi(\mathbf{r},t)]$. The order parameter can be written as a uniform field plus small fluctuations:

[TABLE]

Our mean field plus fluctuation theory can be viewed as an expansion in terms of the Ginzburg parameter $G$, or equivalently, the small noise term $\eta(\mathbf{r},t)$. The uniform background is the zeroth order term in the random noise $\eta(\mathbf{r},t)$. The TDGL equation thus leads to the coupled equations of the uniform background and the fluctuations

[TABLE]

where $()_{k}$ mean the Fourier component with momentum $k$, $\eta_{ik}(t)$ means the $k$ momentum component of the noise, $j\neq i$ represents the other order different from order $i$. Multiplying the second equation by $\psi_{-k}(t)$ and taking the average over the probability functional $\rho$, one obtains the equation of motion for the second moment:

[TABLE]

If one keeps only $O(\eta^{2})$ terms, the equation for the second moment $\langle\psi_{i}(k)\psi_{j}(-k)\rangle$ simplifies to

[TABLE]

where

[TABLE]

and repeated indices should be summed over. Note that $T_{v}=T/(E_{c}V)$ is the temperature normalized to the condensation energy of the whole volume. Therefore, the fluctuation just evolves in a time dependent quadratic potential determined by the curvature of the local free energy landscape taken at $\bar{\psi}(t)$. In principle, one could make a ‘mean field’ approximation by assuming the probability function $\rho$ is always a Gaussian function and easily take into account the fluctuation correction to Eq. (64). But this is out of the scope of this paper.

The linearized Fokker-Plank equation close to point $O$ reads

[TABLE]

where $\alpha$ and $T_{v}$ are time dependent and $\alpha_{k}(t)=\alpha(t)-\xi_{0}^{2}k^{2}$. This is a diffusion equation for the probability $\rho(\psi_{k})$ in the quadratic potential $-\alpha_{k}\psi(k)^{2}$ with diffusion constant $\gamma T_{v}$. The exact solution to Eq. (66) is the Gaussian function Eq. (5) with the variance satisfying (here and below we denote $\psi(k)$ as $\psi_{k}$ for simplicity)

[TABLE]

The solution to Eq. (67) is

[TABLE]

which could be obtained by squaring

[TABLE]

where $G_{k}(t,t^{\prime})=\Theta(t-t^{\prime})e^{S_{k}(t,t^{\prime})}$ is the Green’s function for Eq. (8) of the main text and $S_{k}(t,t^{\prime})=2\gamma\int_{t^{\prime}}^{t}dx\alpha_{k}(x)$ is the accumulated exponent.

Appendix C Dynamics

We evaluate Eq. (29) of the main text, which together with the initial condition term is

[TABLE]

where the accumulated exponent is

[TABLE]

and $T$ should be understood as the dimensionless quantity $T/(E_{c}V)$ here and in the following. The integral in Eq. (12) describes the contributions to the variance $D_{k}(t)$ from noise fluctuations that are created at time $t^{\prime}$ and then propagated forward by the equation of motion. The pump and cooling profile determines the time dependence of $D_{k}$ and the needed expressions may be straightforwardly evaluated for any pump and cooling profile. Within the linear cooling profile approximation of Figure 2, $S_{k}$ and $T$ are combinations of quadratic, linear and constant functions of time and analytic results can be written down in terms of error functions, Gaussians and exponentials. We present here further analytical work based on the linear cooling profile that brings insight.

At times $t<t_{0}$, $S_{k}<0$ for all $k$ so fluctuations created at a time $t<t_{0}$ decay as time increases to $t_{0}$. At times $t>t_{0}$ long wavelength fluctuations ($k^{2}\xi_{0}^{2}<\alpha(t)$) increase exponentially with time. Thus $t_{0}$ is a convenient reference point and we are interested in times $t$ greater than $t_{0}$. Separating the integral into times greater and less than $t_{0}$ and noting that $S(t,t^{\prime})=S(t,t_{0})+S(t_{0},t^{\prime})$ and that $S(t,t^{\prime})=-S(t^{\prime},t)$ we have

[TABLE]

with

[TABLE]

and

[TABLE]

The first term ($D^{(1)}$) describes the contribution to the variance of fluctuations created before $t_{0}$ and propagated forward to $t$ and the second term ($D^{(2)}$), which we have rearranged for later convenience, describes the additional contributions of fluctuations occurring after $t_{0}$. We are interested in the case in which the fluctuations at $t=t_{0}$ are very small, and we wish to focus on long times such that the growing modes have increased to an amplitude of the order of unity. In this circumstance a general asymptotic analysis is possible but for ease of writing we will focus on the linear cooling profile for which

[TABLE]

and

[TABLE]

We begin with $D^{(1)}$ which we rewrite as

[TABLE]

where $D_{H}$ describes the propagation forward in time of the fluctuations existing before the pump was turned on and created by the pump. Using the linear cooling profile formulas and the definition of $S$

[TABLE]

The requirement that the mean field order parameter be completely suppressed means that $e^{-2|\alpha_{H}|(2t_{pump}+t_{0})}\ll G/\alpha_{L}$ so

[TABLE]

which represents the hot thermal fluctuations created by the pump propagated to $t=t_{0}$.

We now turn to $D_{KZ}$ which represents the fluctuations created as the system cools from $t=0$ to $t=t_{0}$ after the pump is turned off:

[TABLE]

where in the second equality we have defined $u=(t_{0}-t)/t_{0}$. In the rapid cooling limit $|\alpha_{H}|\gamma t_{0}\ll 1$ $D_{KZ}$ is $\mathcal{O}(\gamma t_{0})$ and is much smaller than $D_{H}$. In the slow cooling limit the integral is dominated by small $u$, and we may extend the upper limit to infinity and rescale $u=v/(\sqrt{2|\alpha_{H}|\gamma t_{0}})$ obtaining

[TABLE]

where we have defined the important length and time scales

[TABLE]

We see that for $k^{2}\xi_{KZ}^{2}\ll 1$, $D_{KZ}\approx\frac{T_{c}}{2\alpha_{KZ}}\sqrt{\pi}$ and for $k^{2}\xi_{KZ}^{2}\gg 1$, $D_{KZ}\sim\frac{T_{c}}{2k^{2}\xi_{0}^{2}}$. This is the expected behavior of a critical theory with mean field exponents and an effective distance from criticality determined by the cooling rate, consistent with the general analysis of Kibble and Zurek Kibble (1976); Zurek (1996, 1985).

We now present a qualitative evaluation of $D^{(2)}(t)$. We have $S_{k}(t,t_{0})=S_{0}(t,t_{0})-2\xi_{0}^{2}k^{2}(t-t_{0})$ with

[TABLE]

We are interested in growing modes, for which $\frac{S_{0}(t,t_{0})}{t-t_{0}}>k^{2}\xi_{0}^{2}$; roughly these are those for which $\alpha_{L}>k^{2}\xi_{0}^{2}$. In the slow cooling case $\alpha_{L}\gamma(t_{m}-t_{0})=\frac{\alpha_{L}^{2}}{|\alpha_{H}|+\alpha_{L}}\gamma t_{m}>1$. In the fast cooling limit, we may set $t_{m}-t_{0}=0$ and write

[TABLE]

Here we have neglected $k$-dependence, which is on a scale that is not relevant at large enough $t$. In the ultra slow cooling case we have for $t<t_{m}$ and defining $u=\frac{t^{\prime}-t_{0}}{t_{m}-t_{0}}$

[TABLE]

with $u_{max}=\frac{t-t_{0}}{t_{m}-t_{0}}$. The argument of the exponential is maximized at $u=u^{\star}=\frac{\xi_{0}^{2}k^{2}}{\alpha_{L}}$ (for growing modes $u_{max}>2\xi_{0}^{2}k^{2}/\alpha_{L}$ so $u^{\star}$ is within integration range). Defining $u=u^{\star}+\frac{v}{2\sqrt{\gamma(t_{m}-t_{0})}}$ and noting that $t_{m}-t_{0}=t_{m}\frac{\alpha_{L}}{|\alpha_{H}|+\alpha_{L}}$ we have (after integrating over $v$ using saddle point approximation)

[TABLE]

Note that we have kept only half of the Gaussian integral. This is valid for the modes $\xi_{0}^{2}k^{2}\lesssim\frac{1}{4\gamma(t-t_{0})}\ll\frac{1}{\sqrt{8\gamma(t_{m}-t_{0})/\alpha_{L}}}$ which are the only relevant ones at long time $t-t_{0}$. The exponent is also small in this limit so to an adequate approximation we have

[TABLE]

which is the same as the $k\rightarrow 0$ limit of $D^{(1)}_{k}\approx D_{KZ}$. This is expected since $D^{(2)}_{k}$ can be interpreted as the fluctuations created after $t_{0}$ and back propagated to $t_{0}$. This is symmetric to $D^{(1)}_{k}$ for $k\rightarrow 0$ in the slow cooling limit.

Appendix D Exact solution in terms of error functions

The $D_{k}^{(1)}=D_{H}+D_{KZ}$ in Eq. (72) has the interpretation of the fluctuation prepared at time $t_{0}$. The exact form of $D_{H}$ is Eq. (86). In the linear cooling profile approximation used in this paper, the exact form of $D_{KZ}$ is

[TABLE]

where $\mathrm{Erf}[x_{1},x_{2}]=\frac{2}{\sqrt{\pi}}\int_{x_{1}}^{x_{2}}e^{-x^{2}}dx$ is the error function.

For the $D_{k}^{(2)}$ term, it is simpler to neglect the time dependence of the noise, i.e., take $T_{c}=T_{L}$, after which one obtains

[TABLE]

where we have defined

[TABLE]

One can take various limits of Eq. (97) and (98) to get the results in the previous section.

Appendix E The fast cooling limit: $t_{m}=0$

For $t_{m}=0$ and neglecting the initial condition at $t=t_{pump}$, Eq. (12) reduces to

[TABLE]

Summing up contribution from all the Fourier modes, one obtains the real space correlation function

[TABLE]

At long time $4\gamma\alpha_{l}t\gg 1$, $D_{k}$ is approximately a gaussian function in $k$ and the Fourier transform becomes

[TABLE]

where $\xi(t)=\xi_{0}\sqrt{8\gamma t}$ is the universal correlation growth law and

[TABLE]

is the Ginzburg parameter for critical phenomenon at equilibrium. Thus the fluctuation $\langle\psi_{i}(0)^{2}\rangle$ grows exponentially with time and $\psi_{1}$ grows faster due to a larger $\alpha_{1L}\gamma_{1}$. Setting $\langle\psi(0)^{2}\rangle=\alpha_{L}/2$ gives the crossover time

[TABLE]

At this time, keeping the first two terms in the $\ln$ expansion of $t_{c}$, the ratio between the fluctuations in the two directions is

[TABLE]

Appendix F The slow cooling case in competing order systems

The slow cooling case is characterized by a large $t_{m}$. After $t_{1}$, $\psi_{1}$ starts to grow exponentially while $\psi_{2}$ has been growing for a time of $t_{1}-t_{2}$. Assume both order parameters are in the exponential growing stage and nonlinearity is not yet on set, they obey the equation

[TABLE]

The crossover of $\psi_{i}$ to nonlinearity happens at $\left<\psi_{i}(0)^{2}\right>_{t}=\alpha_{i}(t)$ which yields

[TABLE]

where $x_{i}=2\gamma_{i}(\alpha_{iL}+|\alpha_{iH}|)(t-t_{i})^{2}/t_{m}$ and

[TABLE]

To leading order in the $\ln$ expansion we have for order I:

[TABLE]

At time $t_{c}$, the correlation length of the fluctuations is $\xi\approx 2\sqrt{2}\left(\frac{\gamma t_{m}}{2(\alpha_{1L}+|\alpha_{1H}|)}\ln\frac{1}{\zeta}\right)^{1/4}\xi_{0}\sim\left(\gamma t_{m}\ln\frac{1}{\zeta}\right)^{1/4}\xi_{0}$, which predicts a logarithmic correction to the $\xi\sim t_{m}^{1/4}$ Kibble-Zurek scaling Zurek (1985). Note that the logarithmic correction could be numerically large for very small $\zeta$. In 2D, the number density of topological vortices created is thus $n\sim 1/\xi(t_{c})^{2}\sim(\gamma t_{m}\ln\frac{1}{\zeta})^{-1/2}\xi_{0}^{-2}$. For the mean field plus fluctuation theory to be valid, we also require $G\alpha(t)^{D/2-2}\ll 1$ at the predicted $t_{c}$ for $D<4$, in other words:

[TABLE]

which yields

[TABLE]

Trapping into the metastable minimum requires that $\left\langle\psi_{2}(0)^{2}\right\rangle\ll\left\langle\psi_{1}(0)^{2}\right\rangle$ at $t_{m}$. Simply comparing the exponents yields

[TABLE]

where $\lambda_{d}=(t_{1}-t_{2})/t_{m}$. This imposes the criterion

[TABLE]

By assuming $\zeta=10^{-4}$, $\gamma_{1}=2\,\mathrm{ps}^{-1}$, $(\alpha_{1l},\alpha_{2l})=(1,1.1)$ and $(\alpha_{1h},\alpha_{2h})=(-1,-0.9)$, one obtains $t_{mu}\approx 4.6\,\mathrm{ps}$ and $t_{ms}\approx 460\,\mathrm{ps}$. Since the typical cooling time due to electron phonon thermalization ranges from $1\,\mathrm{ps}$ to $100\,\mathrm{ps}$, most ultrafast experiments are in the regime analyzed in this paper (Fig. 5).

Appendix G The pumping process

The pump brings $\alpha_{iL}$ to $\alpha_{iH}<0$ as shown in Fig. 2 which induces the order parameter dynamics from point II to $O$ in Fig. 1(a). At mean field level, this nonlinear dynamics is described by Eq. (62) and the uniform component obeys the exact solution

[TABLE]

The long time asymptotic form is $\bar{\psi}_{2}^{2}=\frac{-\alpha_{2H}/2}{(1-\alpha_{2H}/\alpha_{2L})}e^{4\alpha_{2H}\gamma_{2}(t+t_{pump})}$ which means $\bar{\psi}_{2}$ approaches zero exponentially but never reaches it during finite amount of time. At time zero, the pump is removed and $\bar{\psi}_{2}$ reaches a small value

[TABLE]

We first consider the fast cooling limit $t_{m}=0$. After time zero, $\bar{\psi}_{2}$ goes back towards minimum II following the dynamics

[TABLE]

Thus at time $t_{c}$ when $\psi_{1}$ fluctuation crossovers to nonlinearity, $\bar{\psi}_{2}$ has recovered by an exponential factor. The more accurate picture for the probability distribution is that of Fig. 3(a) but shifted in $\psi_{2}$ direction by the amount of $\bar{\psi}_{2}^{2}(t_{c})$. For the probability of trapping into phase I to be still close to one, we require that

[TABLE]

which yields

[TABLE]

and further leads to the criterion

[TABLE]

for the pump pulse in the leading order. Despite the logarithmic factor, this time scale can be made small with a larger $|\alpha_{2H}|$, i.e., a stronger pump will prepare the $\bar{\psi}_{20}$ with a smaller value at time zero.

If the cooling rate is finite, it is simpler to consider the case $t_{m}>t_{mu}$ such that the crossover happens at $t_{c}$ before $t_{m}$. The pumping time is effectively longer than $t_{pump}$ in this case since the suppression process of $\psi_{2}$ lasts until $t_{2}$, when $\alpha_{2}(t)$ crosses zero. Applying Eq. (117) to this case yields

[TABLE]

to leading order. For sufficiently large $|\alpha_{iH}|$, the right hand side of Eq. (120) becomes negative and thus the criterion is satisfied by any $t_{pump}>0$. This is because in the cooling process before $t_{2}$, the high temperature stage already suppresses $\psi_{2}$ well enough.

In realistic situation, the pump might not be strong enough and the proportion of phase I domains created, $p_{1}$, can be calculated as a function of pump fluence/duration. It should crossover sharply from [math] to $1$ at the boundary of Eq. (119) or Eq. (120) depending on which regime the cooling rate is in.

Appendix H Joint probability function

The probability that $\psi(0)=A$ and $\psi(r)=B$ is

[TABLE]

where

[TABLE]

and $\mathcal{N}$ is the normalization of the functional integral and $\delta$ is a functional delta function. Representing the delta functions by integrals gives

[TABLE]

or, Fourier transforming the real-space $\psi$

[TABLE]

We can now perform the integral over the $\psi_{k}$ and arrive at

[TABLE]

The sum over $k$ results in

[TABLE]

where

[TABLE]

at $t=t_{c}$. We finally perform the $\lambda$ integrals, getting

[TABLE]

Appendix I The coefficient $\lambda$ and $\vartheta$

The coefficients are:

[TABLE]

and

[TABLE]

Note that $\alpha_{i}$ should be interpreted as $\alpha_{iL}$ in the fast cooling limit.

Appendix J Non-equilibrium phase diagrams

‘Nonequilibrium’ phase diagrams Fig. 10 can be drawn for $\alpha_{iH}$. The system is originally at the blue dot $\alpha_{iL}$. If the pump brings the $\alpha_{iH}$ to any of the colored regions, metastable trapping into the SC (I) state could happen. However, different regions have different stories as described in the main text. For example, it $\alpha_{iH}$ is in region $1$ and $t_{pump}+t_{0}/2$ is much larger than $t_{d}=\frac{\Delta}{4|\alpha_{2H}|\gamma_{2}}\ln\frac{1}{\zeta_{1}}$, the system can be brought to disordered state by the pump. The subsequent dynamics of fluctuation will lead the system into the metastable SC state if the relaxation in the SC direction is substantially faster.

Appendix K Nano Granules

If the system is a nano granule whose size is smaller than a coherence length, one can neglect the spatial fluctuation and treat the order parameter as uniform, i.e., one could keep the $k=0$ mode only. The initial dynamics is an expansion of the Gaussian probability due to thermal noise, regardless of the flow direction due to the potential. After the time $4\alpha\gamma t\sim 1$, the dynamics starts to be dominated by the flow. At this time, $\psi_{i}^{2}\sim T_{v}/\alpha_{i}$ and is much smaller than $\alpha_{i}$ if the system volume $V$ is not too small. Thus nonlinearity is not yet onset. After passing the crossover point $\psi_{i}^{2}\sim T_{v}/\alpha_{i}$ to flow dynamics, the order parameter in each basin will be finally attracted to the corresponding minima, as shown in Fig. 4. One immediately observes that if $\psi_{i}$ distribution is still tiny at the crossover point, most of the $\psi$ lies inside basin I due to the nearly vertical shape of the basin boundary close to the origin.

To estimate of probability of trapping into phase I, one can draw a rectangle centered at $O$ with half lengths of $L_{i}=\sqrt{\langle\psi_{i}^{2}\rangle}$. The length of its edge embedded in basin II is $l_{2}=4\lambda L_{2}^{\gamma_{1}\alpha_{1}/(\gamma_{2}\alpha_{2})}$ while the total length is $l=4L_{1}+4L_{2}$. Therefore, the probability of trapping into phase I is roughly

[TABLE]

where $\delta=\frac{1}{2}\left(\gamma_{1}\alpha_{1}/(\gamma_{2}\alpha_{2})-1\right)>0$ and $\kappa$ is order one. Since $T_{v}$ is a very small number, this probability is almost unity. After trapped into it, the life time of this metastable state is is exponentially large: $T_{life}\sim\frac{1}{\gamma}e^{U/T_{v}}$ where $U$ is the dimensionless energy barrier between the global and metastable minima. The detailed calculation for the lifetime is described by Kramer’s theory Kramers (1940); Landauer and Swanson (1961).