# Dynamics of topological defects after a photo-induced melting of a charge-
density wave

Andrei E. Tarkhov<EMAIL_ADDRESS>Laboratory for the Physics of
Complex Quantum Systems, Moscow Institute of Physics and Technology,
Dolgoprudny 141700, Russia  Division of Genetics, Department of Medicine,
Brigham and Women’s Hospital and Harvard Medical School, Boston, MA 02115, USA
Skolkovo Institute of Science and Technology, Moscow 121205, Russia  A.V.
Rozhkov  Laboratory for the Physics of Complex Quantum Systems, Moscow
Institute of Physics and Technology, Dolgoprudny 141700, Russia  Institute for
Theoretical and Applied Electrodynamics, Russian Academy of Sciences, Moscow
125412, Russia  Alfred Zong  Massachusetts Institute of Technology, Department
of Physics, Cambridge, MA 02139, USA  University of California at Berkeley,
Department of Chemistry, Berkeley, CA 94720, USA  Anshul Kogar  Massachusetts
Institute of Technology, Department of Physics, Cambridge, MA 02139, USA  Nuh
Gedik  Massachusetts Institute of Technology, Department of Physics,
Cambridge, MA 02139, USA  Boris V. Fine<EMAIL_ADDRESS>Laboratory
for the Physics of Complex Quantum Systems, Moscow Institute of Physics and
Technology, Dolgoprudny 141700, Russia  Skolkovo Institute of Science and
Technology, Moscow 121205, Russia  Institute for Theoretical Physics,
University of Leipzig, Leipzig 04103, Germany

###### Abstract

Charge-density-wave order in a solid can be temporarily “melted” by a strong
laser pulse. Here we use the discrete Gross-Pitaevskii equation on a cubic
lattice to simulate the recovery of the CDW long-range phase coherence
following such a pulse. Our simulations indicate that the recovery process is
dramatically slowed down by the three-dimensional topological defects – CDW
dislocations – created as a result of strongly nonequilibrium heating and
cooling of the system. Overall, the simulated CDW recovery was found to be
remarkably reminiscent of a recent pump-probe experiment in LaTe3.

Topological defects in the form of vortices play a prominent role in
nonequilibrium and equilibrium statistical physics of ordered systems [1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. In particular, vortices are known to
emerge through the Kibble-Zurek mechanism [14, 15, 16, 17], when the system is
cooled across a phase transition. In the Kibble-Zurek setting, the system is
initially at equilibrium, and the cooling process is supposed to be adiabatic
everywhere except for a close vicinity of the phase transition temperature. It
is still an outstanding question: What happens after a quench across a phase
transition, where the initial state is far from equilibrium, and the energy
flow from the system is so fast that the system does not reach quasi-
equilibrium at any time during the quench?

In the present work, we address the above question in the context of the
recent pump-probe experiment of Ref. [18] on the light-induced melting of a
charge-density wave (CDW) in $\text{LaTe}_{3}$. The experiment monitored the
recovery of the CDW order after its destruction by a laser pulse. The
experimental system was initially in the low-temperature ordered phase; then
the part of the system responsible for the CDW order was rapidly heated up and
cooled down, so that it did not have time to reach the high-temperature
equilibrium. It was observed [18] that the amplitude and the phase of the CDW
responded qualitatively differently: the amplitude recovered quickly, whereas,
for sufficiently high pump pulse intensity, the recovery of phase coherence
took much longer [19]. It was conjectured in Ref. [18] on the basis of
indirect experimental evidence that the observed slowdown of the CDW-order
recovery is due to the slow dynamics of the photoinduced topological defects –
CDW dislocations. Topological defects within the CDW phase were also invoked
in other contexts in Refs. [20, 21, 22, 23, 24, 25, 26, 27].

a

bCDW

DGPE

Vorticity

Figure 1:  Sketch of the mapping of a CDW dislocation
$u\cos\left[\mathbf{q}_{\rm CDW}\mathbf{r}_{m}+\phi(\mathbf{r}_{m})\right]$
onto a vortex of the DGPE $U(1)$ lattice variables
$u\exp\left[i\phi(\mathbf{r}_{m})\right]$. (a) 2D slice of the CDW dislocation
(left) vs. 2D slice of the DGPE vortex (right). (b) Three-dimensional
representations of a vortex loop: CDW (left) vs. DGPE lattice (middle) vs.
vorticity of the DGPE lattice defined in the text (right).

The translational symmetry breaking associated with the onset of a CDW can be
equivalently viewed as the breaking of the $U(1)$ symmetry with respect to the
phase of the CDW order parameter. In such a picture, the CDW dislocations are
mapped onto vortices in the underlying $U(1)$-ordered state (see Fig. 1).

In the present work, we simulate the principal features of the CDW phase
dynamics with the help of the discrete Gross-Pitaevskii equation (DGPE), which
also exhibits a $U(1)$ phase transition. We implement the heating-cooling
quench, perform the numerical imaging of vortices, trace their dynamics and,
thereby demonstrate that they are responsible for the slowdown of the order
parameter recovery. The basic character of our simulations suggest that the
observed non-equilibrium phenomenology is applicable to a broad class of
physical systems.

General considerations. — CDW is a modulation of electronic charge density
accompanied by a periodic modulation of ionic positions

$\mathbf{r}_{m}\to\mathbf{r}_{m}+u\,\mathbf{e}\cos\left(\mathbf{q}_{\rm
CDW}\mathbf{r}_{m}+\phi\right),$ (1)

where $\mathbf{r}_{m}$ is the position of the $m$’th ion, $\mathbf{q}_{\rm
CDW}$ the CDW wave vector, $u$ and $\phi$ the amplitude and the phase of the
CDW, and $\mathbf{e}$ the unit vector along the CDW lattice displacement
direction [19]. One can choose the complex number $u\exp(i\phi)$ as the CDW
order parameter. By allowing this order parameter to vary slowly in space and
time, $u(\mathbf{r},t)\exp\left(i\phi(\mathbf{r},t)\right)$, one can then
describe the low-frequency dynamics of a CDW.

In Fig. 2(a), we reproduce one of the experimental results of Ref. [18],
namely, the time dependence of the 2D-integrated intensity of the CDW
diffraction peak for different fluences of the laser pump pulse. This
intensity is proportional to the squared amplitude of the CDW order parameter.
For smaller pulse fluences, the order parameter initially becomes smaller but
never vanishes entirely and then quickly recovers to almost the initial value.
For higher fluences, the CDW first disappears, i.e. becomes totally melted,
but then it recovers over much longer times, in fact, not reaching the
equilibrium value on the timescale accessible in the experiment. The
subsequent theoretical analysis of Ref. [19] concluded that the above slow
recovery cannot be explained by the amplitude dynamics of the CDW order (see
also [28, 29]). The amplitude takes about $2$ ps to relax [18, 19] It is the
phase relaxation happening afterwards that we aim at describing in this paper.

a

b

Figure 2: (a) Time evolution of experimentally measured 2D-integrated
intensity of the CDW electron diffraction peak after photo-excitation, for
several laser photo-excitation densities (from Ref. [18]). (b) Numerically-
calculated 2D-integrated intensity
$I_{2}\equiv\sum_{k_{x},k_{y}}|\psi(k_{x},k_{y},0)|^{2}$ tracking the ordering
in the system [see Eq. (6)] as a function of time for different quench
intensities. The color of the curves encodes different quench intensities as
indicated on the colorbar. Each numerical data point represents the average of
$40$ quench realizations, and the error bars represent the statistical
uncertainty associated with this average. Inset: Variation of the energy
density during quenches. Horizontal line marks critical energy density
$\varepsilon_{c}$, see Fig. 3(a), above which the order is destroyed.

The model. — Our numerical model, DGPE, is a space-discretized version of the
continuous Gross-Pitaevskii equation (GPE). Both can be considered as possible
dynamical extensions of the static XY-model. The GPE describes the low-
temperature dynamics of superfluid Bose gases, while the DGPE is often used to
describe bosons on optical lattices. We note that both the spatially
inhomogeneous CDW order $u(\mathbf{r})\exp\left(i\phi(\mathbf{r})\right)$ and
the DGPE share the $U(1)$ structure of the dynamical variable, and both also
exhibit a spontaneous symmetry breaking transition into a $U(1)$-ordered
state. The actual dynamics of the inhomogeneous CDW is different and, in many
respects, more complex than that of DGPE. Real CDWs are dissipatively coupled
with the electrons and with the lattice degrees of freedom, while the DGPE
dynamics is time-reversible. Yet, due to the large number of degrees of
freedom involved in the simulations, DGPE delivers an effectively dissipative
environment for the smaller number of variables associated with the formation
of topological defects over larger spatial scales. It can be reasonably
expected that the detailed character of the thermal bath is not crucial here,
and, therefore, the thermal bath of the real CDW system can be efficiently
replaced by the bath associated with the microcanonical thermalizing dynamics
of the DGPE lattice.

To connect the CDW dynamics and the DGPE, we divide the crystal lattice into
cells of sizes $l_{x}\times l_{y}\times l_{z}$, where $l_{x}$, $l_{y}$, and
$l_{z}$ are of the order of the CDW coherence lengths for the respective
directions. The mapping is achieved by associating the DGPE lattice variable
$\psi_{j}$ with the average of the CDW order parameter over the CDW cell
around position $\mathbf{r}_{j}$:

$\psi_{j}\sim\left<u(\mathbf{r})e^{i\phi(\mathbf{r})}\right>_{j\text{th
cell}}.$ (2)

Given the crystal anisotropy of LaTe3, the coherence length for the direction
perpendicular to the Te layers, $l_{z}$, should be significantly smaller than
$l_{x}$ and $l_{y}$. We expect that $l_{x}$, $l_{y}$ and $l_{z}$ can be chosen
such that the renormalized phase stiffness between adjacent anisotropic unit
cells can be modelled with the help of the same constant $J$ [30, 31, 32] for
all three directions. Each cell is thus represented by a single site of the
cubic DGPE lattice. The modelling is classical, since, for large enough cells,
the quantum fluctuations can be neglected.

The time evolution of $\psi_{j}$ is to be modelled by the DGPE on a cubic
lattice:

$i\frac{d\psi_{j}}{dt}=-J\sum_{k\in\text{NN}(j)}\psi_{k}+U\left|\psi_{j}\right|^{2}\psi_{j},$
(3)

where the summation over $k$ runs through all nearest-neighbors $\text{NN}(j)$
of site $j$, coefficient $J$ is the “hopping” parameter and $U$ is the
interaction parameter. The lattice has $V$ sites. In the simulations, the
typical lattice dimensions were $50\times 50\times 50$. DGPE conserves the
total energy of the lattice
$E\equiv\sum_{j}\left(-J\sum_{k\in\text{NN}(j)}\psi_{k}\psi^{*}_{j}+\frac{U}{2}\left|\psi_{j}\right|^{4}\right)$
as well as the “total norm” $N\equiv\sum_{j}|\psi_{j}|^{2}$. We also define
the energy density $\varepsilon\equiv E/V-U/2$ and the “norm density” $n\equiv
N/V$. The character of the DGPE solutions is determined by the dimensionless
parameter $g\equiv Un/J$, which, for the reasons explained later, was chosen
to be equal to 10. Below, except for the period of quench, we simulate an
isolated DGPE lattice, which, when perturbed, is observed to exhibit dynamic
thermalization to a microcanonical temperature $T$ that was determined
numerically as in Ref. [33, 34]. The computed temperature is a monotonic
function of $\varepsilon$. For $g=10$, $J=1$, $n=1$, the DGPE exhibits a
spontaneous symmetry breaking transition at the critical temperature $T_{\rm
c}\approx 4.25$ into low-temperature ordered state characterised by the $U(1)$
order parameter
$\Psi(t)=|\Psi(t)|e^{i\phi(t)}\equiv\frac{1}{\sqrt{n}V}\sum_{j}\psi_{j}(t)$.
In equilibrium, $|\Psi|=|\Psi_{\rm eq}|(T)$, where $|\Psi_{\rm eq}|(T)$ is a
decreasing function of $T$ for $T<T_{\rm c}$, vanishing above $T_{\rm c}$.
Both $\varepsilon(T)$ and $|\Psi_{\rm eq}|^{2}(T)$ are plotted in Fig. 3(a).
The heat capacity $c_{v}(T)\equiv d\varepsilon/dT$ is shown in Fig. 3(b); it
exhibits a very recognizable cusp at $T=T_{\rm c}$ corresponding to
$\varepsilon=\varepsilon_{c}$.

Figure 3: Temperature dependence of various parameters for the 3D DGPE
lattice: (a) the equilibrium order parameter $|\Psi_{\rm{eq}}|^{2}$ and the
energy density $\varepsilon$, (b) specific heat $c_{v}$.

In principle, Eq. (3) can be generalized to include disorder representing the
CDW pinning centers. However, in the target experiment [18], the slowdown of
the recovery of the CDW coherence was observed on the spatial scale where no
evidence of the presence of the pinning centers was detected.

Parameters. — To finalize the mapping, we provide the order-of-magnitude
estimates for the three parameters of the cubic DGPE lattice, namely, $J$,
$n$, and $U$.

Since the CDW coherent domains are supposed to be larger than the crystal unit
cell, the CDW phase dynamics should occur at typical frequencies smaller than
the Debye frequency $\Omega_{\rm D}$. Hence, the paramater $J$, which
determines the characteristic frequency of the DGPE dynamics, can be roughly
estimated as

$J=\Omega_{\rm D}\frac{a}{l},$ (4)

where $a\equiv(a_{x}a_{y}a_{z})^{1/3}$ and $l\equiv(l_{x}l_{y}l_{z})^{1/3}$
are the geometrical averages of the crystal lattice periods
$a_{x},a_{y},a_{z}$ and the CDW coherence lengths $l_{x},l_{y},l_{z}$,
respectively. We use the mean-field theory [35], which is supposed to be
rather crude for LaTe3, to make the order-of-magnitude estimate $l=\hbar
v_{\rm F}/(1.76\pi k_{\rm B}T_{\rm c})$, where $v_{\rm F}=2\varepsilon_{\rm
F}/p_{\rm F}$ is the Fermi velocity, and $T_{\rm c}$ is the CDW phase
transition temperature. Thus $\varepsilon_{\rm F}\approx 1.5$ eV [36], $p_{\rm
F}\approx\frac{3\pi}{8}\frac{\hbar}{a}$, $T_{c}\sim 700$ K. As a result, we
get $l\sim 6a$. Inserting this together with $\Omega_{\rm D}\sim 4$ THz into
Eq.(4) gives $J\sim 4\text{ ps}^{-1}$.

The typical value for $n$ can be derived from the fact that, at the critical
temperature, the hopping energy is of order of the thermal energy $Jn\sim
k_{\rm B}T_{\rm c}$. Thus:

$n\sim\frac{k_{\rm B}T_{\rm c}}{J}\sim\hbar\frac{T_{\rm c}}{T_{\rm
D}}\frac{l}{a}.$ (5)

At the last step, we replaced $\hbar\Omega_{\rm D}$ in Eq.(4) with
$k_{B}T_{\rm D}$, where $T_{\rm D}\sim 200\rm{K}$ is the Debye temperature for
LaTe3 [37]. Thus $n/\hbar\sim 20$, which is consistent with the classical
character of our modelling.

Parameter $U$ can be determined by matching the speed of phasons, $c_{\rm ph}$
with the speed of low-frequency excitations of the DGPE lattice: $c_{\rm
ph}=l\sqrt{2JUn}$ [32]. In turn, $c_{\rm ph}$ can be crudely estimated as the
speed of sound [38]. It is further shown in [32] that, under these
assumptions, $U\sim 10\frac{J}{n}$. Therefore, our simulations were done with
$g=Un/J=10$.

a

b

Figure 4: (a,b) Snapshots of vorticity for two individual quenches at time
$t=20J^{-1}$: (a) weak quench, (b) strong quench, marked in panels (c,d) by
gray and red triangles, respectively [also shown in Fig. 2(b)]. Vorticity
$\mathbf{w}_{j}$ is depicted by (colored) arrows whose length is proportional
to $|\mathbf{w}_{j}|$ (see the text). For the visualization purposes, we apply
a smoothing filter to $\mathbf{w}_{j}$: at site $j$ we plot $\mathbf{w}_{j}$
averaged over $9$ sites comprising a cube with side $3$ encompassing site $j$.
The arrow color represents the position along one of the lattice axes, as
shown by the color bar. In (a) these arrows are visible as mere dots due to
smallness of $|\mathbf{w}_{j}|$. In (b) longer and hence brighter arrows merge
into distinctly visible worm-like bundles revealing vortex cores. Vortex lines
penetrating the entire system are clearly seen for the stronger quench.
(Further snapshots of vorticity for the later times of the stronger quench can
be found in the Supplement [32].) (c,d) Correlation between vorticity and the
slowing down of the order recovery for the simulated 3D DGPE lattice. (c) the
2D-integrated intensity $I_{2}(t)$, same as Fig. 2(b) but over a longer time
interval; (d) the average vorticity $Q_{av}(t)$. The color-coding and the
statistical ensemble behind the sampling are used the same as in Fig. 2(b).
(d, inset) Dependence of the order parameter recovery time $\tau_{r}$ on the
linear lattice size $V^{1/3}$.

Simulations of non-equilibrium quenches. — According to Ref. [19], the energy
deposited by the laser pulse into LaTe3 is first absorbed by the electrons,
then transferred partially to the “hot phonons”, drastically increasing their
temperature and thereby melting the CDW order. After that, the energy of the
hot phonons is transferred to the rest of the phonons, which have
significantly larger specific heat. This brings the overall temperature of the
system to nearly the same temperature as the one before the pulse. The phase
degrees of freedom we aim at describing with the help of the DGPE are supposed
to belong to hot phonons. We model the above process by first quickly pumping
the energy into the DGPE lattice (heating of hot phonons by electrons) and
then quickly removing the energy from the lattice (cooling of hot phonons by
the energy transfer to the main phonon bath). Technically the above scheme is
implemented by temporarily adding the energy-nonconserving terms to the right-
hand-side of Eq. (3) as explained in Refs. [34, 39, 32].

We simulated a set of quenches starting from the equilibrium state at
$T\approx 0.4T_{c}$ and then using varying amounts of energy deposited to and
then removed from the system, and averaged over $40$ random initial
conditions, thermalized before the quench for $t=600J^{-1}$ [32]. The time
dependencies of energy during the quenches is shown in the inset of Fig. 2(b).
These quenches imitated the experimental laser pulses of different fluences
used to produce the experimental plots in Fig. 2(a).

To compare our simulation with the diffraction experiment [18], we define the
two-dimensional integrated spectral weight $I_{2}(t)$ by the following
equality

$\displaystyle I_{2}=\sum_{k_{x},k_{y}}|\psi(k_{x},k_{y},0)|^{2},$ (6)

where $\psi({\bf k})=V^{-1}\sum_{j}\psi_{j}e^{i{\bf k}\cdot{\bf r}_{j}}$ is
the Fourier transform of $\psi_{j}$. Summation in Eq. (6) runs over quantized
momentum projections $k_{x,y}$ lying in the Brillouin zone, while $k_{z}=0$.
Definition (6) mimics the 2D-integrated diffraction intensity measured
experimentally in Ref. [18], see Fig. 2 (a). The square of the order parameter
amplitude $|\Psi|^{2}$ is the largest contributing term in the sum (6).
However, $I_{2}>|\Psi|^{2}$ since shorter-range order fluctuations enhance
$I_{2}$ relative to $|\Psi|^{2}$. Thus, one can think of $I_{2}$ as a quantity
tracking both long-range and short-range correlations in the DGPE system. We
have further verified [32] that $I_{2}(t)$ is proportional to the correlation
length of the CDW order, $l_{c}(t)$, as anticipated in Refs. [18, 19].

The main panel of Fig. 2(b), presents $I_{2}(t)$ before, during and after the
quench. The comparison between the simulation and the experiment reveals a
remarkable agreement: At small intensities of the quench, the spectral
intensity quickly recovers its pre-quench equilibrium value on the time scale
comparable to the duration of the quench itself. As the quench intensity
grows, the character of the dynamics undergoes a qualitative change: $I_{2}$
drops to zero and stays low for finite time interval. This happens when,
during the quench, the energy crosses the critical value $\varepsilon_{c}$
associated with the phase transition and indicated by the dashed line in the
inset of Fig. 2(b). Most importantly, the post-quench recovery is extremely
slow: the corresponding curves visually split away from the curves
representing weaker quenches.

Topological defects. — Here we investigate the possible connection between the
slowdown of the order recovery and the dynamics of topological defects. For 3D
DGPE systems with periodic boundary conditions, the topological defects form
closed-loop vortices [40]. In this respect, the 3D DGPE model is similar to
the three-dimensional $O(2)$, $U(1)$ and XY-models, where the percolation of
vortex loops is responsible for a phase transition [41, 42, 43].

To monitor the local topological charge (vorticity), we define the dual
lattice by translating the original lattice by vector $(1/2,1/2,1/2)$. Each of
the sites of the dual lattice is thus surrounded by $6$ square plaquettes of
the original lattice.

As an indicator of vorticity we use the finite-difference counterpart
$\mathbf{w}_{j}$ of the continuous curl of the current
$(\nabla\times\mathbf{v}_{s})_{j}$, where
$\mathbf{v}_{s}=-i(\psi^{*}\nabla\psi-\psi\nabla\psi^{*})$ (see [32]). In Fig.
4(a,b), we provide a space-resolved snapshots of field $\mathbf{w}_{j}$
generated in the course of the quenches shown in Fig. 2(b). The snapshots
correspond to the two states marked in Fig. 2(b) by red and gray triangles. In
principle, each site of the dual lattice has a non-zero vorticity
$\mathbf{w}_{j}$, which is the consequence of the discreteness of the lattice
field. However, in the absence of a real vortex, such a vorticity is very
small and does not form any line pattern, while the presence of a vortex is
indicated by the large values of $\mathbf{w}_{j}$ showing the vortex core as a
colored “tube” extended through the lattice.

Fig. 4(a) illustrates a weak quench, for which the order (quantified by
$I_{2}$) recovers quickly. One can see that the system does not exhibit any
colored tube indicative of a vortex. For more intensive quenches, on the
contrary, one can observe vortex tubes percolating through large volumes of
the system and forming clearly identifiable large vortex loops, see Fig. 4(b)
(see also Fig. S.7(a-c) and the video in the supplemental material [32]). As
one can see, for example, in the video, these loops initially form an
entangled network. The crossings of the loop lines then lead to smaller loops
separating from this network and collapsing, while one large loop spreading
over most of the lattice is eventually being formed. The time required for
that last loop to shrink and collapse is what determines the time scale of the
order recovery in Fig. 2(b).

We also introduce an auxiliary binary vorticity variable $Q_{j}$ for each dual
lattice site [32]: it is defined such that $Q_{j}=1$, if at least one of the 6
adjacent original-lattice plaquettes has non-zero vorticity; otherwise
$Q_{j}=0$. Finally, we define the average vorticity for the entire system as
$Q_{av}=\frac{1}{V}\sum_{j}^{V}Q_{j}$.

Further evidence that the relaxation of vorticity $Q_{av}(t)$ is directly
correlated with the order recovery after a quench is presented in Figs.
4(c,d): there the strong quenches exhibiting the slowdown of the $I_{2}$
recovery are accompanied by pronounced slowly relaxing tails of $Q_{av}(t)$.

In Fig. 4(c, inset), we show that the order parameter recovery $\tau_{r}$,
defined as the time to reach $75\%$ of the initial equilibrium value of
$|\Psi|^{2}$ for the strongest quench, grows with the lattice volume (number
of lattice sites) $V$. The simulations are consistent with the scaling
$\tau_{r}\propto V^{2/3}$. We have not attempted to derive this power law.
Yet, the growth of $\tau_{r}$ with the lattice size as such indicates that the
recovery of the order is contingent on the disappearance of the extensive
vortex loops that penetrate the entire system and need to drift by the
distances of the order of the lattice size before they annihilate.

Conclusions — In conclusion, our direct simulations of the DGPE lattice with a
CDW-like form of the order parameter are consistent with the conjecture of
Ref. [18] that the observed slowdown of the recovery of the CDW order after
laser-induced melting is due to the emergence and then slow disappearance of
topological defects in the order parameter texture. We were able to
numerically monitor the vorticity in the system and thereby establish the
correlation between the relaxation of vorticity and the recovery of the order.
In a broader context, our investigation illustrates the viability of using
DGPE for simulating the CDW dynamics even though DGPE as such was historically
introduced to describe superfluid systems.

Code availability — The source code is published in a GitHub repository [44].

Acknowledgments. — This work was supported by a grant of the Russian Science
Foundation (Project No. 17-12-01587). The work at MIT was mainly funded by the
U.S. Department of Energy, BES DMSE (data taking and analysis) and the Gordon
and Betty Moore Foundation’s EPiQS Initiative grant GBMF9459
(instrumentation).

## References

  * Berezinsky [1970] V. Berezinsky, Zh. Eksp. Teor. Fiz. 32, 493 (1970).
  * Berezinsky [1972] V. Berezinsky, Zh. Eksp. Teor. Fiz. 61, 610 (1972).
  * Kosterlitz and Thouless [1972] J. M. Kosterlitz and D. Thouless, J. Phys. C: Solid State Phys. 5, L124 (1972).
  * Svistunov [1991] B. Svistunov, J. Moscow Phys. Soc 1, 373 (1991).
  * Kagan  _et al._ [1992] Y. M. Kagan, B. Svistunov, and G. Shlyapnikov, JETP 75, 387 (1992).
  * Kagan and Svistunov [1994] Y. Kagan and B. Svistunov, Zh. Eksp. Teor. Fiz 105, 353 (1994).
  * Berloff and Svistunov [2002] N. G. Berloff and B. V. Svistunov, Phys. Rev. A 66, 013603 (2002).
  * Yoshino  _et al._ [2002] H. Yoshino, K. Hukushima, and H. Takayama, Phys. Rev. B 66, 064431 (2002).
  * Proctor  _et al._ [2014] T. C. Proctor, D. A. Garanin, and E. M. Chudnovsky, Phys. Rev. Lett. 112, 097201 (2014).
  * Garanin and Chudnovsky [2015] D. A. Garanin and E. M. Chudnovsky, The European Physical Journal B 88, 81 (2015).
  * Kobayashi and Cugliandolo [2016] M. Kobayashi and L. F. Cugliandolo, Phys. Rev. E 94, 062146 (2016).
  * Kumar  _et al._ [2017] M. Kumar, S. Chatterjee, R. Paul, and S. Puri, Phys. Rev. E 96, 042127 (2017).
  * Vasin and Vinokur [2019] M. Vasin and V. Vinokur, Physica A 525, 1161 (2019).
  * Kibble [1976] T. W. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Kibble [1980] T. W. Kibble, Phys. Rep. 67, 183 (1980).
  * Zurek [1985] W. H. Zurek, Nature 317, 505 (1985).
  * Zurek [1996] W. H. Zurek, Phys. Rep. 276, 177 (1996).
  * Zong  _et al._ [2019a] A. Zong, A. Kogar, Y.-Q. Bie, T. Rohwer, C. Lee, E. Baldini, E. Ergeçen, M. B. Yilmaz, B. Freelon, E. J. Sie, H. Zhou, J. Straquadine, P. Walmsley, P. E. Dolgirev, A. V. Rozhkov, I. R. Fisher, P. Jarillo-Herrero, B. V. Fine, and N. Gedik, Nat. Phys. 15, 27 (2019a).
  * Dolgirev  _et al._ [2020] P. E. Dolgirev, A. V. Rozhkov, A. Zong, A. Kogar, N. Gedik, and B. V. Fine, Phys. Rev. B 101, 054203 (2020).
  * Trigo  _et al._ [2021] M. Trigo, P. Giraldo-Gallo, J. N. Clark, M. E. Kozina, T. Henighan, M. P. Jiang, M. Chollet, I. R. Fisher, J. M. Glownia, T. Katayama, P. S. Kirchmann, D. Leuenberger, H. Liu, D. A. Reis, Z. X. Shen, and D. Zhu, Phys. Rev. B 103, 054109 (2021).
  * Wandel  _et al._ [2020] S. Wandel, F. Boschini, E. H. da Silva Neto, L. Shen, M. X. Na, S. Zohar, Y. Wang, G. B. Welch, M. H. Seaberg, J. D. Koralek, G. L. Dakovski, W. Hettel, M.-F. Lin, S. P. Moeller, W. F. Schlotter, A. H. Reid, M. P. Minitti, T. Boyle, F. He, R. Sutarto, R. Liang, D. Bonn, W. Hardy, R. A. Kaindl, D. G. Hawthorn, J. S. Lee, A. F. Kemper, A. Damascelli, C. Giannetti, J. J. Turner, and G. Coslovich, “Light-enhanced charge density wave coherence in a high-temperature superconductor,”  (2020), arXiv:arxiv:2003.04224 [cond-mat.supr-con] .
  * Kogar  _et al._ [2020] A. Kogar, A. Zong, P. E. Dolgirev, X. Shen, J. Straquadine, Y.-Q. Bie, X. Wang, T. Rohwer, I.-C. Tung, Y. Yang, R. Li, J. Yang, S. Weathersby, S. Park, M. E. Kozina, E. J. Sie, H. Wen, P. Jarillo-Herrero, I. R. Fisher, X. Wang, and N. Gedik, Nat. Phys. 16, 159 (2020).
  * Mihailovic [2019] D. Mihailovic, Appll. Sci. 9 (2019), 10.3390/app9050890.
  * Vogelgesang  _et al._ [2018] S. Vogelgesang, G. Storeck, J. Horstmann, T. Diekmann, M. Sivis, S. Schramm, K. Rossnagel, S. Schäfer, and C. Ropers, Nature Physics 14, 184 (2018).
  * Duan  _et al._ [2021] S. Duan, Y. Cheng, W. Xia, Y. Yang, C. Xu, F. Qi, C. Huang, T. Tang, Y. Guo, W. Luo, _et al._ , Nature 595, 239 (2021).
  * Storeck  _et al._ [2020] G. Storeck, J. G. Horstmann, T. Diekmann, S. Vogelgesang, G. von Witte, S. V. Yalunin, K. Rossnagel, and C. Ropers, Structural Dynamics 7, 034304 (2020), https://doi.org/10.1063/4.0000018 .
  * Zong  _et al._ [2021] A. Zong, A. Kogar, and N. Gedik, MRS Bull. 46, 720 (2021).
  * Maklar  _et al._ [2021] J. Maklar, Y. W. Windsor, C. W. Nicholson, M. Puppin, P. Walmsley, V. Esposito, M. Porer, J. Rittmann, D. Leuenberger, M. Kubli, M. Savoini, E. Abreu, S. L. Johnson, P. Beaud, G. Ingold, U. Staub, I. R. Fisher, R. Ernstorfer, M. Wolf, and L. Rettig, Nat. Commun. 12, 2499 (2021).
  * Zong  _et al._ [2019b] A. Zong, P. E. Dolgirev, A. Kogar, E. Ergeçen, M. B. Yilmaz, Y.-Q. Bie, T. Rohwer, I.-C. Tung, J. Straquadine, X. Wang, Y. Yang, X. Shen, R. Li, J. Yang, S. Park, M. C. Hoffmann, B. K. Ofori-Okai, M. E. Kozina, H. Wen, X. Wang, I. R. Fisher, P. Jarillo-Herrero, and N. Gedik, Phys. Rev. Lett. 123, 097601 (2019b).
  * Cardy [1996] J. Cardy,  _Scaling and renormalization in statistical physics_ , Vol. 5 (Cambridge university press, 1996).
  * Carlson  _et al._ [1999] E. Carlson, S. Kivelson, V. Emery, and E. Manousakis, Physical review letters 83, 612 (1999).
  * [32] Supplemental Material to this paper.
  * de Wijn  _et al._ [2015] A. S. de Wijn, B. Hess, and B. V. Fine, Phys. Rev. E 92, 062929 (2015).
  * Tarkhov [2020] A. E. Tarkhov,  _Ergodization dynamics of the Gross-Pitaevskii equation on a lattice_ , Ph.D. thesis, Skolkovo Institute of Science and Technology (2020).
  * Gruner [2018] G. Gruner,  _Density waves in solids_ (CRC press, 2018).
  * Yao  _et al._ [2006] H. Yao, J. A. Robertson, E.-A. Kim, and S. A. Kivelson, Phys. Rev. B 74, 245126 (2006).
  * Ru and Fisher [2006] N. Ru and I. Fisher, Phys. Rev. B 73, 033101 (2006).
  * Manley  _et al._ [2018] M. E. Manley, P. J. Stonaha, D. L. Abernathy, S. Chi, R. Sahul, R. P. Hermann, and J. D. Budai, Nat. Commun. 9, 1 (2018).
  * Tarkhov  _et al._ [2021] A. E. Tarkhov, A. Rozhkov, and B. V. Fine, arXiv preprint arXiv:2109.00961 (2021).
  * Rojas  _et al._ [2001] F. Rojas, S. Puri, and A. Bray, J. Phys. A: Math. Gen. 34, 3985 (2001).
  * Hulsebos [1994] A. Hulsebos, Nuclear Physics B-Proceedings Supplements 34, 695 (1994).
  * Kajantie  _et al._ [2000] K. Kajantie, M. Laine, T. Neuhaus, A. Rajantie, and K. Rummukainen, Phys. Lett. B 482, 114 (2000).
  * Kohring  _et al._ [1986] G. Kohring, R. E. Shrock, and P. Wills, Phys. Rev. Lett. 57, 1358 (1986).
  * [44] The code is published in a GitHub repository at https://github.com/TarkhovAndrei/DGPE.
  * Rugh [1997] H. H. Rugh, Phys. Rev. Lett. 78, 772 (1997).
  * Rugh [1998] H. H. Rugh, J. Phys. A: Math. Gen. 31, 7761 (1998).
  * Kamaletdinov and Berloff [2021] A. Kamaletdinov and N. G. Berloff, “Quantized vortices mediated annealing,”  (2021), arXiv:2109.05867 [cond-mat.dis-nn] .

##  I Supplemental Material

###  I.1 S.I. Details of the experiment of Ref.[18]

Here, for the convenience of the readers, we provide additional details of the
experiment of Ref.[18].

Figure S.1: Thermalization dynamics behind Fig. 3, $V=50^{3}$: energy density,
order parameter and the right-hand-side of Eq. (S.2)

The material studied, LaTe3, has a layered crystal structure consisting of
alternating Te bilayers and LaTe layers. The critical temperature of the CDW
transition in LaTe3 is about $700$ K. The CDW is localized primarily in the Te
bilayers.

In the experiment of Ref.[18], femtosecond laser pulses were applied to LaTe3
and then the CDW response was probed by (i) optical reflectivity measurements,
(ii) time-resolved angle-resolved photoemission spectroscopy and (iii) the
ultrafast electron diffraction. Fig. 2(a) of the main text shows the results
for the latter technique, namely, the 2D-integrated intensity, $I_{2}$, of the
CDW electron diffraction peak is plotted as a function of delay time for
various laser pulse fluences. The peak intensity is integrated for the two in-
plane wave-vector projections. The intensity is interpreted as being
proportional to the square of the CDW order-parameter multiplied by the CDW
correlation length for the out-of-plane direction (the one perpendicular to
the Te layers).

By examining the time evolution of the width of the CDW peak, the authors
Ref.[18] provided indirect evidence that the slow recovery of the intensity of
the CDW diffraction peak seen in Fig. 2(a) is caused by the emergence of
topological defects of the CDW order, followed by their subsequent slow
annihilation. Even though such a scenario is rather natural, no direct
evidence in real space in its favor was presented.

###  I.2 S.II. Microcanonical thermodynamics of the DGPE lattice: Numerical
aspects

###  I.3 Equilibrium properties

For the parameters used in the main text, the solutions of DGPE with fixed $E$
and $N$ were observed to exhibit the behavior consistent with the ergodic
hypothesis. Namely, the right-hand-side of Eq. (S.2) computed for a given
phase space trajectory and used below to define the temperature was found to
converge to a value independent of the initial conditions.

To describe the statistical properties of the DGPE lattice, we use the
standard definition of the microcanonical temperature:

$\frac{1}{T}\equiv\frac{\partial S}{\partial E},$ (S.1)

where $S\equiv\text{log}\,W(E)$ is the microcanonical entropy, with $W(E)$
being the $(V-2)$-dimensional volume of the energy shell at constant $N$. We
set the Boltzmann constant $k_{\rm B}=1$.

Figure S.2: Thermalization of pre-quench states for $40$ on-shell initial
conditions generated by the procedure described in Section S.II with
$V=50^{3}$.

The direct determination of the volume $W(E)$ in the many-dimensional phase
space is a computationally prohibitive task. It is, however, possible to use
Monte-Carlo sampling to determine directly $\frac{1}{W(E)}\frac{dW(E)}{dE}$,
which is equal to the right-hand-side of Eq.(S.1). We do this by adapting the
numerical recipe developed for classical spin lattices [33] and consistent
with the analytical approach of Refs. [45, 46]. The idea is based on the
observation that the rate of growth of a many-dimensional energy shell is
determined by the average local curvature of that shell. The latter, in turn,
can be locally sampled in the vicinity of each point by introducing small
isotropic random perturbations to the phase point coordinate and then
computing the variance of energy fluctuations $\langle\Delta E^{2}\rangle$ and
the mean fluctuation $\langle\Delta E\rangle$. We note that $\langle\Delta
E\rangle\neq 0$ due to the fact that, in general, there are exponentially more
states above the energy shell than below it. We use the above approach with
the additional constraint that the small perturbations should leave $N$
constant. The corresponding average is denoted as $\langle\ldots\rangle_{N}$.
With such a notation, the numerically computed expression for temperature
becomes

$\frac{1}{T}=\frac{2\langle\Delta E\rangle_{N}}{\langle\Delta
E^{2}\rangle_{N}}.$ (S.2)

In the limit $V\to\infty$, computing the averages entering Eq.(S.2) in the
vicinity of one randomly chosen sample point on an energy shell should be
sufficient to determine the temperature. However, since our numerically
simulated lattice is finite, we supplement the above procedure by the
additional averaging of the right-hand-side of Eq.(S.2) over points along a
randomly chosen phase space trajectory and then further over phase space
trajectories with randomly chosen initial conditions.

Figure S.3: Quench realizations behind the averages plotted in Figs. 2(b) and
4(a) – $40$ initial conditions, $V=50^{3}$. For stronger quenches $\kappa=1.0$
and $\kappa=1.3$, atypical slowly relaxing solutions appear – they indicate
the presence of large topological defects in the lattice. We do not always
observe the annihilation of defects, in some cases the transient resurgence of
vorticity occurs (black arrows indicate such quench realizations).

Technically, the calculation of equilibrium properties plotted in Fig. 3
consists of the following steps:

  1. (i)

choosing the initial configurations $\\{\psi_{j}\\}$ as a random Hilbert-space
vector normalised by the condition $\sum_{j}|\psi_{j}|^{2}\equiv N=V$. This
state has randomly chosen phases and amplitudes, and corresponds to
$\varepsilon\sim 0$ and $T\sim T_{c}$;

  2. (ii)

choosing an energy value for which the temperature is to be calculated;

  3. (iii)

activating the non-conservative term in the DGPE Eq. (I.5) with function
$K\equiv\operatorname{const}$ for a time interval sufficient to reach a
desired energy value (see Section S.III);

  4. (iv)

deactivating the non-conservative term, then calculating the conservative
dynamics of the DGPE for a sufficiently long time $t=100J^{-1}$ (see Fig.
S.1);

  5. (v)

choosing the next slightly smaller target energy value and then, starting from
the thermalized state at the final time point from the preceding step,
repeating steps (ii) to (v).

The overall procedure guarantees that the DGPE lattice has enough time to
thermalize at each energy shell.

###  I.4 Preparation of the initial state preceding the quenches.

In Fig. S.2, we illustrate the preparation of the initial thermalized state
for the quenches presented in Fig. 2(b) of the main text of the paper. We
start from the high-temperature state of step (i) above and then quickly drive
the system to $\varepsilon=-4.1$ corresponding to $T\approx 0.4T_{c}$. Such a
driving initially produces a strongly nonequilibrium state (in contrast to the
slow preparation procedure used in the previous subsection). Then, we track
the thermalization dynamics of $\varepsilon$,
$\frac{|\Psi|^{2}}{|\Psi_{\rm{eq}}|^{2}}$ and $\frac{2\langle\Delta
E\rangle_{N}}{\langle\Delta E^{2}\rangle_{N}}$ from Eq. (S.2). We observe that
the energy is conserved, the order parameter reaches its equilibrium value by
$t\sim 150J^{-1}$, while the parameter $\frac{2\langle\Delta
E\rangle_{N}}{\langle\Delta E^{2}\rangle_{N}}$ starts exhibiting the typical
behavior faster, namely, after times of the order of $t\sim 10J^{-1}$. By time
$t=200J^{-1}$, all the parameters reach their equilibrium regimes
independently of the initial conditions. For the quenches, presented in Fig.
2(b) of the main text, the conservative dynamics were run for $t=600J^{-1}$,
which is significantly longer than required for equilibration and
thermalization. Hence, before the start of the quenches in Fig. 2(b) the
system was well-thermalized, and the initial conditions for the quenches are
typical for the energy shell.

###  I.5 S.III. Numerical implementation of the DGPE lattice quench

The quench is performed on the DGPE lattice prepared in equilibrium at
$\varepsilon=-4.1$, $T/T_{c}\approx 0.4$. During the quench, the system’s
energy experiences quick rise and fall within a limited time interval. To
allow for the time variation of $E$, we modify [34, 39, 47] DGPE by adding a
non-energy-conserving term (the last one below):

$\displaystyle
i\frac{d\psi_{j}}{dt}=-\sum_{k\in\text{NN}(j)}\psi_{k}+\frac{Un}{J}\left|\psi_{j}\right|^{2}\psi_{j}-$
$\displaystyle-K(t)\psi_{j}D_{j},$ (S.3)

where

$\displaystyle
D_{j}^{\vphantom{*}}=D_{j}^{*}=i\sum_{k\in\text{NN}(j)}\left(\psi_{k}^{*}\psi_{j}-\psi_{k}\psi_{j}^{*}\right),$
(S.4)

and the function $K(t,\kappa)$ is to be specified below. Note that Eq. (I.5),
remains “gauge invariant” and thus conserves the norm $N$.

To verify the energy-nonconserving character of the added term one can use Eq.
(I.5) to obtain

$\dot{E}=-K(t)\sum_{j}D_{j}^{2}.$ (S.5)

The sign of $K(t)$ controls the increase or decrease of $E$: function $K$ must
be negative to heat the system, and positive to cool it. We additionally
require that the energies of the system before and after the quench are
identical.

All our simulated quenches use the function $K(t)$ having the following form

$\displaystyle K(t)$ $\displaystyle=$ $\displaystyle-\theta(t)\left[\ \kappa\
\frac{\gamma_{2}e^{-\gamma_{2}t}-\gamma_{1}e^{-\gamma_{1}t}}{\gamma_{2}-\gamma_{1}}\right.$
$\displaystyle-$
$\displaystyle\left.\theta(t-t^{*})\gamma_{3}\left(1-e^{-\gamma_{2}(t-t^{*})}\right)\frac{E(t)-E_{0}}{E^{*}-E_{0}}\right],$

where $\theta(t)$ is the Heaviside step-function, $\kappa$ is the parameter
controlling the strength of the quench (in our simulations, $\kappa$ spans the
range from 0.1 to 20), $\gamma_{1}=1$, $\gamma_{2}=0.3$ and $\gamma_{3}=0.01$
are numerical parameters,

$\displaystyle
t^{*}=\frac{1}{\gamma_{1}-\gamma_{2}}\log\left(\frac{\gamma_{1}}{\gamma_{2}}\right)$
(S.7)

is the time at which $K(t)$ changes sign, $E_{0}$ is the pre-quench energy of
the system, $E(t)$ is the energy at time $t$, and $E^{*}=E(t^{*})$.

The quench is principally determined by the first term in the square brackets
in Eq.(I.5) – the time-dependence of this term is shown in Fig. S.4. The
second term in the square brackets is added to guarantee that the energy after
the quench is equal to the initial energy $E_{0}$.

Before the quench, the system is prepared in an equilibrium low-energy state
as described in Section S.II. At the first stage of the quench the system is
subjected to the fast energy increase (heating). The duration that stage is
$t^{*}$. The heating is then followed by the cooling stage of duration
regulated by parameters $\gamma_{2}$ and $\gamma$. The resulting time
evolution of the energy density $\varepsilon(t)$ for various quench parameters
$\kappa$ is plotted in the inset of Fig. 2(b) of the main text.

Figure S.4: The principal factor in the first term in Eq. (I.5) as a function
of time. Figure S.5: Time dependence of the correlation length for the quench
illustrated in Fig. 2(b) and Fig. 4(c).

a

b

Figure S.6: System size scaling for the time evolution of (a) the
theoretically computed post-quench recovery, $I_{2}$, and (b) the correlation
length, $l_{c}$. The color of the curves represents different system sizes.

For further details one can refer to the source code published in a GitHub
repository [44].

###  I.6 S.IV. Estimate of parameter $U$ for the DGPE lattice

In order to estimate the nonlinearity parameter $U$, we match the velocity of
the phasons propagating within the real CDW system and the velocity of
acoustic excitations calculated theoretically for the DGPE model. To perform
the desired calculation, we treat phasons as weak perturbations distorting
spatially homogeneous zero-temperature solution $\Psi_{0}(t)$. Within such an
approach, the original DGPE variables $\psi_{j}$ should be expressed as the
sum $\psi_{j}(t)=\Psi_{0}(t)+\xi_{j}(t)$, where variables $\xi_{j}(t)$ are
small. The DGPE now reads

$\displaystyle
i\dot{\Psi}_{0}+i\dot{\xi}_{j}=-2dJ\Psi_{0}-J\sum_{m\in\text{NN}(j)}\xi_{m}+$
(S.8)
$\displaystyle+U\left|\Psi_{0}+\xi_{j}\right|^{2}\left(\Psi_{0}+\xi_{j}\right).$

Neglecting all terms containing $\xi$, we derive that the unperturbed zero-
temperature solution $\Psi_{0}(t)$ must satisfy

$i\dot{\Psi}_{0}=-2dJ\Psi_{0}+U\left|\Psi_{0}\right|^{2}\Psi_{0}.$ (S.9)

The solution to the latter equation can be trivially found

$\Psi_{0}(t)=\sqrt{n_{0}}\exp{(-it\mu)},$ (S.10)

where $n_{0}=|\Psi_{0}|^{2}$ and the “chemical potential” $\mu=-6J+Un_{0}$.

When the deviations from the homogeneous zero-temperature solution (S.10) are
small ($|\xi|\ll\sqrt{n_{0}}$), it is permissible to keep only the leading-
order terms in Eq. (S.8). This allows us to derive the following linear
differential equation

$i\dot{\zeta}_{j}=-J\sum_{m\in\text{NN}(j)}\zeta_{m}-\mu\zeta_{j}+U\left(2n_{0}\zeta_{j}+n_{0}\zeta_{j}^{*}\right).$
(S.11)

The new dynamical variables $\zeta_{j}$ introduced in this equation are
connected to $\xi_{j}$ by the relation $\xi_{j}=\zeta_{j}e^{-it\mu}$.

In Fourier space Eq. (S.11) becomes

$i\dot{\zeta}_{\mathbf{k}}=(\varepsilon_{\mathbf{k}}-\mu+2Un_{0})\zeta_{\mathbf{k}}+Un_{0}\zeta_{\mathbf{-k}}^{*},$
(S.12)

where $\varepsilon_{\mathbf{k}}$ is defined as

$\varepsilon_{\mathbf{k}}=-2J\left[\cos{(k_{x}l)}+\cos{(k_{y}l)}+\cos{(k_{z}l)}\right].$
(S.13)

Since Eq. (S.12) couples $\zeta_{\bf k}$ and $\zeta_{-\bf k}^{*}$, we need to
solve Eq. (S.12) together with the equation for $\zeta_{-\bf k}^{*}$, which
reads

$\displaystyle-i\dot{\zeta}_{\mathbf{-k}}^{*}=(\varepsilon_{\mathbf{k}}-\mu+2Un_{0})\zeta_{\mathbf{-k}}^{*}+Un_{0}\zeta_{\mathbf{k}}.$
(S.14)

Solving the system of equations (S.12) and (S.14), we find the spectrum of
excitations $\omega_{\bf k}$ for the DGPE lattice. Frequencies $\omega_{\bf
k}$ satisfy

${\rm
det}\,\begin{pmatrix}\omega_{\mathbf{k}}-\varepsilon_{\mathbf{k}}+\mu-2Un_{0}&-Un_{0}\\\
-Un_{0}&-\omega_{\mathbf{k}}-\varepsilon_{\mathbf{k}}+\mu-2Un_{0}\end{pmatrix}=0.$
(S.15)

The positive branch of the resulting dispersion relation is

$\omega_{\mathbf{k}}=\sqrt{(6J+\varepsilon_{\mathbf{k}}+Un_{0})^{2}-U^{2}n_{0}^{2}}.$
(S.16)

When $l|{\bf k}|\ll 1$, it can be further approximated as

$\omega_{\mathbf{k}}\approx\sqrt{\left(Jl^{2}\mathbf{k}^{2}\right)^{2}+2l^{2}JUn_{0}\mathbf{k}^{2}}\approx
c_{\rm ex}|{\bf k}|,$ (S.17)

where $c_{\rm ex}=l\sqrt{2JUn_{0}}$ is the speed of the acoustic excitations.

We map the DGPE onto the CDW dynamics by equating $c_{\rm ex}$ with the speed
of CDW phasons $c_{\rm ph}$, thereby arriving to the relation:

$c_{\rm ph}=l\sqrt{2JUn_{0}}$ (S.18)

used in the main text.

We then use Eq. (S.18) to estimate the DGPE nonlinearity term $U$ as follows:

$\displaystyle U=\frac{J}{2n}\left(\frac{c_{\rm
ph}}{Jl}\right)^{2}=\frac{J}{2n}\left(\frac{c_{\rm ph}}{\Omega_{\rm
D}a}\right)^{2}=$ (S.19)
$\displaystyle=(6\pi^{2})^{2/3}\frac{J}{2n}\left(\frac{c_{\rm
ph}}{c_{s}}\right)^{2},$

where $c_{s}\sim(6\pi^{2})^{-1/3}\Omega_{\rm D}a$ is the speed of sound.

As a further crude estimate, we now use $c_{\rm ph}\sim c_{s}$, thereby
obtaining $U\sim 10\frac{J}{n}$. (Ref. [38] showed that the CDW phason
velocity can be $2.8-4.3$ the sound velocity.)

###  I.7 S.V. Definition of the vorticity characteristics

The numerical calculation of the vorticity field $\mathbf{w}_{j}$ for each
dual lattice site is based on the fourth-order finite differentiation scheme
for each first-order derivative
$\left(df(x)/dx\right)_{j}=(f_{j+2}-8f_{j+1}+8f_{j-1}-f_{j-2})/12$.

In order to define the vorticity variable $Q_{j}$, we first introduce the
auxilliary variable $q_{p}$ for each plaquette surrounding the dual site $j$
as follows: (i) The DGPE lattice variables are parameterised as
$\psi_{j}=|\psi_{j}|\exp{(i\phi_{j})}$. (ii) For a given plaquette consisting
of vertices labelled as $1$, $2$, $3$ and $4$, we define

$\begin{split}q_{\text{p}}&=\frac{1}{2\pi}[\left.(\phi_{2}-\phi_{1})\right|_{[-\pi,\pi)}+\left.(\phi_{3}-\phi_{2})\right|_{[-\pi,\pi)}+\\\
&+\left.(\phi_{4}-\phi_{3})\right|_{[-\pi,\pi)}+\left.(\phi_{1}-\phi_{4})\right|_{[-\pi,\pi)}],\end{split}$
(S.20)

where the notation $(...)|_{[-\pi,\pi)}$ indicates that the phase differences
are “backfolded” into the interval $[-\pi,\pi)$. The variable $q_{\text{p}}$
can take only three possible values $-1$, $0$ and $1$. We then define the
variable $Q_{j}$ such that $Q_{j}=0$ if all adjacent plaquettes have
$q_{\text{p}}=0$, and $Q_{j}=1$ otherwise.

###  I.8 S.VI. Relation between $I_{2}$ and the correlation length of the CDW
order

The $\delta$-function peak due to the CDW order is known to be broadened by
the presence of topological defects, which make the peak to acquire a finite
width. The intensity of the $\delta$-peak is then spread over a range of small
wavevectors $q$, which in the real space then leads to a finite correlation
length of the CDW order. The same is supposed to be true in our simulations
around the ordering vector at $q=0$.

We have found empirically, that the low-$q$ intensity in the presence of
vortices scales as

$W(q)=A/q^{2}.$ (S.21)

Such a scaling is naturally expected given that the low-$q$ fluctuations
should be associated with the sound waves of a vortex liquid with dispersion
$\omega\propto q$. The equilibrium thermal amplitude of each mode $a_{q}$ is
supposed to obey the relation $a_{q}^{2}\omega^{2}=T$, from which Eq.(S.21)
follows given also that $W(q)\propto a_{q}^{2}$.

The singular character of the above shape precludes one from defining the
correlation length from the half-width of the peak at half maximum. Instead,
we define the correlation length as

$l_{c}=\frac{\pi}{q_{0}},$ (S.22)

where $q_{0}$ is the upper cutoff value for $W(q)$ given by Eq.(S.21). In
turn, the value of $q_{0}$ can be estimated in the isotropic thermodynamic
limit from the sum rule

$|\Psi_{0}|^{2}=\int_{0}^{q_{0}}W(q)4\pi q^{2}\ dq,$ (S.23)

which is just an expression of the fact that the intensity of the perfect
long-range order $|\Psi_{0}|^{2}$ becomes redistributed over small $q$ due to
the presence of the topological defects.

Equations (S.21), (S.22) and (S.23) lead to

$l_{c}=\frac{4\pi^{2}A}{|\Psi_{0}|^{2}}.$ (S.24)

In our simulations, we track $W(q)$ after a quench as a function of time, fit
the low-$q$ part using Eq.(S.21) and, thereby, obtain the value of the
constant $A$ as a function of time. The latter is then substituted to
Eq.(S.24) to obtain the time dependence of the correlation length.

The resulted time dependencies of $l_{c}$ are plotted for different quench
strengths in Fig. S.5 and for different lattice sizes in Fig. S.6(b). The
comparison of these plots with Figs. 4(c) and S.6(a) respectively illustrates,
that, as anticipated in Refs. [18, 19], the behavior of $I_{2}$ very closely
tracks the behavior of $l_{c}$.

###  I.9 S.VII. Additional snapshots of vorticity after a quench

In Fig. S.7, we present three supplemental snapshots of vorticity for a strong
quench from Fig. 4(b) at times $t=25,30$ and $35J^{-1}$. The subsequent
dynamics of vortices can be observed from these snapshots. The vortices tend
to decrease their size, and slowly collapse and disappear.

a

b

c

Figure S.7: Dynamics of vortices. Snapshots of vorticity for a strong quench
from Fig. 4(b) at times (a) $t=25J^{-1}$, (b) $t=30J^{-1}$ and (c)
$t=35J^{-1}$.

###  I.10 S.VIII. Simulations of anisotropic DGPE lattices

Many compounds demonstrating $U(1)$ ordering possess pronounced quasi-2D
anisotropy. For example, LaTe3 discussed in the main text may be modelled as a
collection of the Te planes where the CDW resides, with weak inter-plane
coupling. For such materials it is natural to inquire to what extent the
presented results are affected by the anisotropy. This question may look
particularly relevant given the incipient Berezinskii-Kosterlitz-Thouless
transition. Our simulations, as well as simulations performed in Ref. [31],
suggest that, unless the anisotropy is extremely high, its effects can be
accounted for by a simple rescaling of the model parameters.

To illustrate this point, we plot the equilibrium value of the order parameter
$|\Psi|$ as a function of temperature $T$, for various anisotropies, see Fig.
S.8 The four panels in the latter figure are plotted for four different values
of the anisotropy ratio $J_{z}/J_{x,y}$, where $J_{x,y,z}$ represent the
coupling constants along the corresponding axes, and $J_{z}\leq J_{x}=J_{y}$.
We notice that the general dependence $|\Psi|=|\Psi(T)|$ maintains its overall
mean-field-like appearance for all values of $J_{z}/J_{x,y}$, with small
corrections beginning to appear for $J_{z}/J_{x,y}=0.1$. More importantly, the
simple rescaling for the $U(1)$ ordering transition temperature in the
anisotropic model

$\displaystyle
T_{c}(J_{z}/J_{x,y})=\frac{1}{3}\left(2+\frac{J_{z}}{J_{x,y}}\right)T_{c}(1),$
(S.25)

works well for all anisotropies we had considered.

As for the Berezinskii-Kosterlitz-Thouless transition, the corresponding
critical point introduces significant new corrections only when the quasi-2D
anisotropy is extremely strong. We believe that our target experimental
material does not enter this regime.

a

b

c

d

Figure S.8: The scaling of critical temperature for the DGPE lattices with the
following anisotropy parameters: (a) $\frac{J_{z}}{J_{x,y}}=0.1$, (b)
$\frac{J_{z}}{J_{x,y}}=0.3$, (c) $\frac{J_{z}}{J_{x,y}}=0.5$, (d)
$\frac{J_{z}}{J_{x,y}}=0.7$. The critical temperature for each case was
estimated within the mean field approximation as
$\left(2+\frac{J_{z}}{J_{x,y}}\right)/3\cdot T_{c}$, where $T_{c}$ is the
critical temperature of the isotropic DGPE lattice.