# Supplemental information for “Fate of transient order parameter domain walls
in ultrafast experiments”

Lingxian Kong  International Center for Quantum Materials, School of Physics,
Peking University, Beijing 100871, China  Ryuichi Shindou  International
Center for Quantum Materials, School of Physics, Peking University, Beijing
100871, China  Zhiyuan Sun  Corresponding author<EMAIL_ADDRESS>State
Key Laboratory of Low-Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing 100084, China

###### Contents

  1. I Brief introduction about the formation of the order-parameter domain wall in ultrafast experiments.
  2. II The first-stage dynamics
    1. II.1 Eigenvalues of the $\hat{L}_{a}$ operators.
    2. II.2 Solving $C(\bm{R},t)=\langle\varphi(\bm{0},t)\varphi(\bm{R},t)\rangle$ based on Eq. (S4).
      1. II.2.1 An equation of motion for $\varphi(\bm{R},t)$.
      2. II.2.2 Derivation of Eqs. (4,6) in the main text.
    3. II.3 Determine the crossover time $t_{\text{c}}$.
  3. III The second-stage dynamics
    1. III.1 Spatial profile of the topological defect without the $U(1)$ symmetry ($\beta\neq 0$).
    2. III.2 Force balance of the defects
      1. III.2.1 2D systems
      2. III.2.2 3D systems
      3. III.2.3 Derivation of Eq. (S64) from Eq. (S63).
  4. IV Effects of the $\partial_{t}^{2}\psi$ term
    1. IV.1 Suppression of the $\partial_{t}^{2}\psi$ term when the temperature is near $T_{\rm c}$.
    2. IV.2 Quantitative modifications of the first and second stage dynamics by the $\partial_{t}^{2}\psi$ term.
    3. IV.3 Derivation of Eq. (S102, S103).
  5. V Coupled three-component field in 3D systems.
    1. V.1 Free energy and first-stage dynamics
    2. V.2 The second-stage dynamics
  6. VI Details of the numerical simulation

##  I Brief introduction about the formation of the order-parameter domain
wall in ultrafast experiments.

In this section, we introduce a mechanism of the domain wall formation [1, 2,
3]. It has been reported both theoretically and experimentally that an optical
pump can invert the charge-density-wave (CDW) order parameter in surface
regions of samples, generating domain walls of the order parameter. The
formation of the domain wall is driven by a thermal effect due to the optical
pump. When the optical pump acts on a sample surface (Fig. S1), layers near
the surface get hotter instantly, while layers inside a sample bulk remains
cold due to a finite penetration depth of the optical pump. Suppose that
during the optical pump, $t_{\text{pump}}<t<t_{\text{off}}$, the temperature
in the surface layers is elevated to a temperature with a spatial dependence
$T(z,t)=T_{0}+\Theta(t-t_{\rm pump})\Theta(t_{\rm
off}-t)(T_{\text{H}}-T_{0})e^{-(z_{\text{top}}-z)/z_{\text{p}}}$. Here
$\Theta(t)$ is the Heaviside step function, $z_{\text{top}}$ is a location of
the top surface, and $z_{\text{p}}$ is the penetration depth of the optical
pump. $T_{\text{H}}$ is the temperature on the top surface of the sample,
which is above an equilibrium critical temperature $T_{\text{c}}$. $T_{0}$ is
the temperature deep inside the sample, which is below $T_{\rm c}$.

Consider that the CDW free energy with the “Mexican hat” potential
$V(|\psi|)=E_{\text{c}}[|\psi|^{2}+\alpha(z,t)]^{2}$ depends on time $t$
through the coefficient $\alpha(z,t)$ and $\alpha(z,t)$ depends on a local
temperature $T(z,t)$, i.e. $\alpha(z,t)=\alpha(T(z,t))$. $\alpha(z,t)$ is
positive/negative if $T(z,t)$ is above/below a transition temperature
$T_{\text{c}}$. For $t<t_{\text{pump}}$, the order parameter in the whole
sample is in the equilibrium energy minima. Without loss of generality, we
assume that $\psi$ for $t<t_{\text{pump}}$ is real, e.g.
$\psi=\psi_{\text{eq}}=-\sqrt{|\alpha|}$. In experiments, the duration of the
optical pump is typically short, so that dynamics of the phase degree of
freedom can be ignored [1]. Therefore, during such a short duration, $\psi$
continues to be real, where the potential depends only on
$\text{Re}\psi\equiv\psi_{1}$ for $t_{\text{pump}}<t<t_{\text{off}}$. Fig. S2
shows schematically an evolution of the potential $V(\psi_{1})$ and how the
order parameter is inverted in the region of the surface layers. After the
pump, $\psi$ in the surface layers follows an oscillatory motion, while $\psi$
deep inside the bulk remains unchanged with
$\psi=\psi_{\text{eq}}=-\sqrt{|\alpha|}$. When the the pump is off at an
appropriate time, $\psi$ in the surface region find its minimum on the other
side of the double-well potential. Consequently, a domain wall shall be formed
at a certain time, $t=t_{0}>t_{\rm off}$. A $\psi$-configuration with single
domain wall at $z=0$ can given by
$\psi(\bm{r})=\sqrt{|\alpha|}\tanh\frac{z}{\xi}$ with
$\xi\equiv\xi_{0}/\sqrt{|\alpha|}$ being the coherence length.

Figure S1: Schematic pictures of the 3D sample under an optical pump, and
time($t$)-dependence of local temperature of a surface region. For
$t_{\text{pump}}<t<t_{\text{off}}$, the optical pump acts on the upper surface
of the sample, elevating the local temperature of the surface region to
$T_{\text{H}}$. For the bulk region away from the surface, the temperature
remains $T_{0}$. Figure S2:  Schematic picture of the order parameter
inversion in the surface region. The curves shows an evolution of a free-
energy potential for $\psi_{1}$ in the surface region,
$V(\psi_{1})=E_{\text{c}}[\psi_{1}^{2}+\alpha(t)]^{2}$. Before the pump
($t<t_{\rm pump}$), and after the pump ($t>t_{\rm off}$), $\alpha(t)$ is
negative, and the function has double wells. During the pump ($t_{\rm
pump}<t<t_{\rm off}$), $\alpha(t)>0$, and the potential is in the single-
minimum form. A ball in the pictures represents the real part $\psi_{1}$ of
the order parameter $\psi$. In the picture, $\psi_{1}$ is in the minimum with
the negative value for $t<t_{\rm pump}$. In the duration of the pump,
$\psi_{1}$ ‘slides down’ along the potential, and changes its sign. When the
pump is turned off, $V(\psi_{1})$ changes into the double-well form again, and
$\psi_{1}$ relaxes into the other well with the positive value.

In the following, we will explain the mechanism of the domain-wall formation
more formally [1, 2, 3]. To describe the oscillatory motion of $\psi$ in the
duration of the optical pump, we add an inertia term into the time-dependent
Ginzburg-Landau (TDGL) equation,

$\displaystyle\frac{1}{\gamma_{0}^{2}}\partial_{t}^{2}\psi+\frac{1}{\gamma}\partial_{t}\psi=\xi_{0}^{2}\nabla^{2}\psi-2\alpha(z,t)\psi-2|\psi|^{2}\psi.$
(S1)

Here we set $\alpha(z,t)=\alpha(T_{0})+\Theta(t-t_{\rm pump})\Theta(t_{\rm
off}-t)\alpha(T_{\rm H})e^{-(z_{\text{top}}-z)/z_{\text{p}}}$. Typically
$z_{\text{p}}\approx 20$nm [1, 2, 3]. By solving Eq. (S1) with the initial
condition $\psi=-\sqrt{|\alpha|}$ before the pump, we confirm that the domain
wall emerges after the duration of the optical pump, as shown in Fig. (S3).

Figure S3: Plot of $\psi$ in the time evolution of Eq. (S1) with the initial
condition $\psi=-1$. The pump is applied in $\gamma t\in[10,10.2]$. i.e.
$\alpha(z,t)=\alpha(T_{0})+\Theta(t-10)\Theta(10.2-t)\alpha(T_{\rm
H})e^{-(z_{\text{top}}-z)/z_{\text{p}}}$. Here $\alpha(T_{0})=-1$,
$\alpha(T_{\rm H})=7$, $z_{\text{p}}/\xi_{0}=20$. $\gamma_{0}=0.1\gamma$.

Several remarks on Eq. (S1) should be noted in order. To begin, the dynamics
of Eq. (S1) is in the mean-field level, which does not include the fluctuation
effect. In the absence of the fluctuation, $\psi$ remains to be real with the
real $\psi$ at the initial time, and the phase dynamics does not develop at
all. Such dynamics is approximately valid within the short duration $t_{\rm
pump}<t<t_{0}$; these durations are anyway short enough that the phase
fluctuation does not develop in practice [1, 3]. To study the further
evolution of the domain wall for $t>t_{0}$, however, it is crucial to include
the two-component nature of $\psi$ as well as the thermal fluctuation effect.
This is the major focus of this work.

Secondly, Eq. (S1) includes the inertia term, i.e. the second-order time-
derivative term $\partial_{t}^{2}\psi$, while the “model A” dynamics that we
mainly use in this work for $t>t_{0}$ does not include the inertia term. To
describe the formation of the domain wall during $t_{\rm pump}<t<t_{0}$, the
inertia term is crucial; without the inertia, $\psi$ in the surface region
falls into the $\psi\approx 0$ state and would stay there for the rest of the
duration. To study the further evolution of the domain wall for $t>t_{0}$, we
will mainly study the “model A” dynamics without the $\partial_{t}^{2}\psi$
term for simplicity. In Sec. IV of this supplemental material as well as the
end of the main text, we investigate the effects of the $\partial_{t}^{2}\psi$
term in the domain wall evolution for $t>t_{0}$. We argue that the inertia
term is suppressed when $T_{0}$ is close to $T_{\rm c}$, modifying only
quantitatively the domain wall evolution for $t>t_{0}$.

##  II The first-stage dynamics

In this section, we will discuss the first-stage dynamics for a general system
with and without the $U(1)$ symmetry. We consider the following free energy
and dynamics with the thermal noise;

$\displaystyle\frac{1}{\gamma}\partial_{t}\psi=-\frac{1}{E_{\text{c}}}\frac{\delta
F}{\delta\psi^{*}}+\eta$ (S2) $\displaystyle F=E_{\text{c}}\int d^{D}\bm{r}\ \
\xi_{0}^{2}|\nabla\psi|^{2}+(|\psi|^{2}+\alpha)^{2}+2\beta\psi_{2}^{2}.$ (S3)
$\displaystyle\langle\eta_{i}(\bm{r},t)\eta_{i^{\prime}}(\bm{r^{\prime}},t^{\prime})\rangle=\frac{2T}{\gamma
E_{\text{c}}}\delta_{ii^{\prime}}\delta(\bm{r}-\bm{r^{\prime}})\delta(t-t^{\prime})$

Here $\gamma^{-1}$ has dimension of time, $\alpha$, $\beta$ and $\psi$ are
dimensionless, and $E_{c}$ has a dimension of energy density. i.e.
[Energy][Length]-D. The system has the $U(1)$-symmetry for $\beta=0$, where
$F$ is invariant under a transformation $\psi\rightarrow\psi^{i\theta}$. The
symmetry reduces to $Z_{2}\otimes Z_{2}$ symmetry for $\beta\neq 0$, where $F$
is invariant under the transformations $\psi_{1(2)}\rightarrow\pm\psi_{1(2)}$.
At $t=t_{0}$, $\psi(\bm{r},t_{0})=\psi_{0}+\delta\psi$ with an initial
fluctuation $\delta\psi=\phi_{1}+i\phi_{2}$. The initial fluctuation comes
from the thermal noise accumulated during the time $t<t_{0}$. We will discuss
the nature of the initial fluctuation in the later section of this section. By
expanding small $\phi_{1}$ and $\phi_{2}$, we obtain Eq. (3) in the main text,

$\displaystyle\frac{1}{\gamma}\partial_{t}\phi_{a}$
$\displaystyle=-(\hat{L}_{a}-\xi_{0}^{2}\nabla_{\bm{R}}^{2})\phi_{a}+\eta_{a},$
(S4)

with $a=1,2$, ${\bm{r}}\equiv({\bm{R}},z)$,
$\nabla_{\bm{R}}^{2}=\nabla^{2}-\partial_{z}^{2}$ is the Laplacian in $D-1$
dimensions. and
$\hat{L}_{1}=-\xi_{0}^{2}\partial_{z}^{2}+2\alpha+6\psi_{0}^{2}$, and
$\hat{L}_{2}=-\xi_{0}^{2}\partial_{z}^{2}+2\alpha+2\beta+2\psi_{0}^{2}$.
According to Eq. (S4), eigenvectors of
$-\hat{L}_{a}+\xi^{2}_{0}\nabla^{2}_{\bm{R}}$ with positive eigenvalues
amplify the fluctuations, dominating the first-stage dynamics for $t>t_{0}$.
The eigenvectors are given by a product of eigenmodes of $\hat{L}_{a}$ and
plane waves $e^{i{\bm{K}}\cdot{\bm{R}}}$. Thus, the amplifying eigenvectors
must comprise of the eigenmodes of $\hat{L}_{a}$ with negative eigenvalues.

###  II.1 Eigenvalues of the $\hat{L}_{a}$ operators.

$\hat{L}_{1}$ and $\hat{L}_{2}$ take forms of one-dimensional scattering
problem with Poschl-Teller potential hole. The scattering problem is defined
by a following Hamiltonian with a parameter $\lambda$,

$\displaystyle\hat{H}(\overline{z};\lambda)=-\partial_{\overline{z}}^{2}-\frac{\lambda(\lambda-1)}{\cosh^{2}(\overline{z})}.$
(S5)

Setting $\overline{z}=z/\xi$, and $\xi\equiv\xi_{0}/\sqrt{|\alpha|}$, we
obtain $\hat{L}_{1}=|\alpha|\hat{H}(\overline{z};3)+4|\alpha|$,
$\hat{L}_{2}=|\alpha|\hat{H}(\overline{z};2)+2\beta$. The scattering problem
has bound-states solutions with negative energies as well as continuum states
with positive energies [4]. The negative eigenvalues are $-(\lambda-1-n)^{2}$
for non-negative integer $n$ smaller than $\lambda-1$: $0\leq n<\lambda-1$.
Thus, $\hat{L}_{1}$ has always positive eigenvalues. $\hat{L}_{2}$ can have
one negative eigenvalue ($-|\alpha|+2\beta$) with an eigenfunction
($\text{sech}(\overline{z})$) for $\beta<|\alpha|/2$.

###  II.2 Solving
$C(\bm{R},t)=\langle\varphi(\bm{0},t)\varphi(\bm{R},t)\rangle$ based on Eq.
(S4).

####  II.2.1 An equation of motion for $\varphi(\bm{R},t)$.

To see how the eigenmode of $\hat{L}_{2}$ with the negative eigenvalue
amplifies the fluctuation, let us set
$\phi_{2}(\bm{r},t)=\varphi(\bm{R},t)\text{sech}(z/\xi)$ with
${\bm{r}}=({\bm{R}},z)$. Eq. (S4) in the main text reads

$\displaystyle\frac{1}{\gamma}\partial_{t}\varphi(\bm{R},t)\text{sech}(z/\xi)=(|\alpha|-2\beta+\xi_{0}^{2}\nabla_{\bm{R}}^{2})\varphi(\bm{R},t)\text{sech}(z/\xi)+\eta_{2}(\bm{r},t)$
(S6)

By multiplying this by $\text{sech}(z/\xi)$ and taking an integral with
$\int\frac{dz}{2\xi}$, we obtain

$\displaystyle\frac{1}{\gamma}\partial_{t}\varphi(\bm{R},t)=(|\alpha|-2\beta+\xi_{0}^{2}\nabla_{\bm{R}}^{2})\varphi(\bm{R},t)+\eta_{\varphi}(\bm{R},t).$
(S7)

with
$\eta_{\varphi}(\bm{R},t)\equiv\int\frac{dz}{2\xi}\text{sech}(z/\xi)\eta_{2}(\bm{r},t)$
and

$\displaystyle\langle\eta_{\varphi}(\bm{R},t)\eta_{\varphi}(\bm{R^{\prime}},t^{\prime})\rangle=\int\frac{dz}{2\xi}\text{sech}(z/\xi)\int\frac{dz^{\prime}}{2\xi}\text{sech}(z^{\prime}/\xi)\langle\eta_{2}(\bm{r},t)\eta_{2}(\bm{r^{\prime}},t^{\prime})\rangle=\frac{T}{\gamma
E_{\text{c}}\xi}\delta(\bm{R}-\bm{R^{\prime}})\delta(t-t^{\prime}).$ (S8)

####  II.2.2 Derivation of Eqs. (4,6) in the main text.

To analyze the dynamics of $\varphi(\bm{R},t)$ for $t>t_{0}$ from Eq. (S7), we
follow the literature [5] and study the Fourier transform of Eq. (S7) with
$\tilde{\varphi}({\bm{K}},t)=\frac{1}{l^{D-1}}\int d^{D-1}{\bm{R}}\ \
\varphi({\bm{R}},t)e^{-i{\bm{K}}\cdot{\bm{R}}}$ and
$\tilde{\eta}_{\varphi}({\bm{K}},t)=\frac{1}{l^{D-1}}\int d^{D-1}{\bm{R}}\ \
\eta_{\varphi}({\bm{R}},t)e^{-i{\bm{K}}\cdot{\bm{R}}}$. 111In this work, we
use the following notation for the Fourier transform.  $\displaystyle f(x)$
$\displaystyle=\sum_{k}\tilde{f}(k)e^{ikx};\ \ \ \ \
\tilde{f}(k)=\frac{1}{l}\int f(x)e^{-ikx}dx.$

$\displaystyle\frac{1}{\gamma}\partial_{t}\tilde{\varphi}(\bm{K},t)=(|\alpha|-2\beta-\xi_{0}^{2}K^{2})\tilde{\varphi}(\bm{K},t)+\tilde{\eta}_{\varphi}(\bm{K},t).$
(S9)

Here $l$ is the system size. A general solution of Eq. (S9) is given by

$\displaystyle\tilde{\varphi}(\bm{K},t)=\tilde{\varphi}(\bm{K},t_{0})e^{S_{K}(t,t_{0})}+\gamma\int_{t_{0}}^{t}\tilde{\eta}_{\varphi}(\bm{K},t^{\prime})e^{S_{K}(t,t^{\prime})}dt^{\prime}$
(S10)

with
$S_{K}(t,t^{\prime})\equiv\gamma\int_{t^{\prime}}^{t}(|\alpha|-2\beta-\xi_{0}^{2}K^{2})ds$.
Note that
$\langle\tilde{\varphi}(-\bm{K},t)\tilde{\varphi}(\bm{K}^{\prime},t)\rangle=\delta_{{\bm{K}},{\bm{K}}^{\prime}}D_{K}(t)$.
An amplitude of $\tilde{\varphi}(\bm{K},t)$,
$D_{K}(t)\equiv\langle\tilde{\varphi}(-\bm{K},t)\tilde{\varphi}(\bm{K},t)\rangle$,
is calculated as

$\displaystyle D_{K}(t)$
$\displaystyle=\frac{1}{l^{D-1}}e^{2S_{K}(t,t_{0})}\Big{(}l^{D-1}D_{K}(t_{0})+\frac{T\gamma}{E_{\text{c}}\xi}\int_{t_{0}}^{t}dt^{\prime}e^{-2S_{K}(t^{\prime},t_{0})}\Big{)}$
(S11)
$\displaystyle=\frac{1}{l^{D-1}}e^{2S_{K}(t,t_{0})}\Big{(}l^{D-1}D_{K}(t_{0})+\frac{T}{2E_{\text{c}}\xi}\frac{1}{|\alpha|-2\beta-\xi_{0}^{2}K^{2}}(1-e^{-2S_{K}(t,t_{0})})\Big{)}.$

$C({\bm{R}},t)\equiv\langle\varphi(0,t)\varphi(\bm{R},t)\rangle$ is given by
an inverse Fourier transform of $D_{K}(t)$,

$\displaystyle C(\bm{R},t)$
$\displaystyle=\sum_{\bm{K}}D_{K}(t)e^{i\bm{K}\cdot\bm{R}}=l^{D-1}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}D_{K}(t)e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=C_{1}(\bm{R},t)+C_{2}(\bm{R},t),$ (S12)

where

$\displaystyle C_{1}(\bm{R},t)$
$\displaystyle\equiv\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K}(t_{0})e^{2(|\alpha|-2\beta)\gamma
t}e^{-2\xi_{0}^{2}\gamma(t-t_{0})K^{2}}e^{i\bm{K}\cdot\bm{R}}$ (S13)
$\displaystyle C_{2}(\bm{R},t)$
$\displaystyle\equiv\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}\frac{T}{2E_{\text{c}}\xi}\frac{e^{2S_{K}(t,t_{0})}-1}{(|\alpha|-2\beta)-\xi_{0}^{2}K^{2}}e^{i\bm{K}\cdot\bm{R}}$
(S14)

$C_{1}$ is dependent on the initial fluctuation at $t=t_{0}$, while $C_{2}$ is
a contribution of the thermal noises that are accumulated during
$t_{0}<t^{\prime}<t$. In the calculation below, we set $t_{0}=0$ for
convenience.

To calculate $C_{2}$ for larger $R$ and $t$, we follow an approximation in
Ref. 5, and obtain a simpler analytical expression. For $\gamma t\gg
1/(|\alpha|-2\beta)$, the first term in Eq. (S14) dominates over the second
term. Also, due to the $e^{-2\xi_{0}^{2}\gamma tK^{2}}$ factor in the
integrand, the integration over $\bm{K}$ converges well within the
$\xi_{0}^{2}K^{2}\ll|\alpha|-2\beta$ region in the $\gamma t\gg
1/(|\alpha|-2\beta)$ limit. Thus, the $K$-dependence in the denominator
becomes unimportant and could be neglected. Thus, we have

$\displaystyle C_{2}(\bm{R},t)$
$\displaystyle\approx\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}\frac{T}{2E_{\text{c}}\xi}\frac{1}{(|\alpha|-2\beta)-\xi_{0}^{2}K^{2}}e^{2(|\alpha|-2\beta)\gamma
t}e^{-2\xi_{0}^{2}\gamma tK^{2}}e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle\approx\frac{T}{2E_{\text{c}}\xi(|\alpha|-2\beta)}e^{2(|\alpha|-2\beta)\gamma
t}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}e^{-2\xi_{0}^{2}\gamma
tK^{2}}e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=\frac{T}{2E_{\text{c}}\xi\xi_{0}^{D-1}(|\alpha|-2\beta)}\frac{e^{2(|\alpha|-2\beta)\gamma
t}}{\left(\sqrt{8\pi\gamma t}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma
t}}$
$\displaystyle=\frac{G}{2\sqrt{|\alpha|}(1-2\beta/|\alpha|)}\frac{e^{2(|\alpha|-2\beta)\gamma
t}}{\left(\sqrt{8\pi\gamma t}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma
t}}$
$\displaystyle=\frac{\zeta|\alpha|}{2(1-2\beta/|\alpha|)}\frac{e^{2(|\alpha|-2\beta)\gamma
t}}{\left(\sqrt{8\pi|\alpha|\gamma
t}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma t}}.$ (S15)

Here we introduced two dimensionless quantities

$\displaystyle G$ $\displaystyle\equiv T/(E_{\text{c}}\xi_{0}^{D}),$ (S16)
$\displaystyle\zeta$ $\displaystyle\equiv G|\alpha|^{D/2-2},$ (S17)

that measures the strength of the thermal fluctuation. Since $t>t_{0}$ in the
calculation above, $T=T_{0}$ in $G$ and $\zeta$ appeared above. Especially,
$\zeta$ plays a similar role as the so-called Ginzburg parameter does in
theories of equilibrium critical phenomenon [6]. A similar $C_{2}(\bm{R},t)$
appears in the literature [5].

To calculate $C_{1}$, we need the expression of the initial fluctuation at
$t=t_{0}$, i.e. $D_{K}(t_{0})$. The fluctuation originates from the thermal
fluctuation at equilibrium for $t<t_{\rm pump}$. It also evolves during the
time domain $t_{\rm pump}<t<t_{0}$. To determine $D_{K}(t_{0})$, we calculate
the thermal fluctuation at equilibrium for $t<t_{\rm pump}$ first.

$\bullet$ Thermal fluctuation at equilibrium : Following the “model A”
dynamics, the thermal fluctuation at equilibrium can be also calculated in
terms of the linearized equation of motion for small imaginary part $\psi_{2}$
around the uniform ground state $\psi=-\sqrt{|\alpha|}$;

$\displaystyle\frac{1}{\gamma}\partial_{t}\psi_{2}$
$\displaystyle=\xi_{0}^{2}\nabla^{2}\psi_{2}-2\beta\psi_{2}+\eta_{2}.$ (S18)

The equation can be solved in the momentum space with
$\tilde{\psi}_{2}({\bm{k}},t)=\frac{1}{l^{D}}\int
d^{D}{\bm{r}}\psi_{2}({\bm{r}},t)e^{-i{\bm{k}}\cdot{\bm{r}}}$ and
$\tilde{\eta}_{2}({\bm{k}},t)=\frac{1}{l^{D}}\int
d^{D}{\bm{r}}\eta_{2}({\bm{r}},t)e^{-i{\bm{k}}\cdot{\bm{r}}}$,

$\displaystyle\frac{1}{\gamma}\partial_{t}\tilde{\psi}_{2}(\bm{k},t)$
$\displaystyle=-(2\beta+\xi_{0}^{2}k^{2})\tilde{\psi}_{2}(\bm{k},t)+\tilde{\eta}_{2}(\bm{k},t),$
$\displaystyle\tilde{\psi}_{2}(\bm{k},t)$
$\displaystyle=\tilde{\psi}(\bm{k},-\infty)e^{-\gamma(2\beta+\xi_{0}^{2}k^{2})\int_{-\infty}^{t}dt^{\prime}}+\gamma\int_{-\infty}^{t}\tilde{\eta}_{2}(\bm{k},t^{\prime})e^{-\gamma(2\beta+\xi_{0}^{2}k^{2})(t-t^{\prime})}dt^{\prime}.$
(S19)

Here we regard that the linear system size along $z$ direction is also given
by $l$. Taking the correlation at $t=-\infty$ to be zero, we obtain

$\displaystyle\langle\tilde{\psi}_{2}(\bm{k},t)\tilde{\psi}_{2}(\bm{k^{\prime}},t)\rangle=\delta_{\bm{k},-\bm{k^{\prime}}}\frac{T}{l^{D}E_{\text{c}}(2\beta+\xi_{0}^{2}k^{2})}$
(S20)

and

$\displaystyle\langle\psi_{2}(\bm{r},t)\psi_{2}(\bm{r^{\prime}},t)\rangle$
$\displaystyle=\sum_{\bm{k},\bm{k^{\prime}}}e^{i\bm{k}\cdot\bm{r}}e^{i\bm{k^{\prime}}\cdot\bm{r^{\prime}}}\langle\tilde{\psi}_{2}(\bm{k},t)\tilde{\psi}_{2}(\bm{k^{\prime}},t)\rangle=\frac{2T\gamma}{E_{\text{c}}}\int\frac{d^{D}\bm{k}}{(2\pi)^{D}}e^{i\bm{k}(\bm{r}-\bm{r^{\prime}})}\frac{1}{2\gamma(2\beta+\xi_{0}^{2}k^{2})}$
(S21)

From Eq. (S21), we obtain

$\displaystyle\langle\psi_{2}(\bm{r},t)\psi_{2}(\bm{r},t)\rangle$
$\displaystyle=\begin{cases}\frac{T}{E_{\text{c}}}\int^{1/a}_{0}\frac{k^{2}dk}{2\pi^{2}}\frac{1}{(2\beta+\xi_{0}^{2}k^{2})},&\text{in
3D,}\\\
\frac{T}{E_{\text{c}}}\int^{1/a}_{0}\frac{kdk}{2\pi}\frac{1}{(2\beta+\xi_{0}^{2}k^{2})},&\text{in
2D.}\end{cases}$ (S22)

where $1/a$ stands for a cutoff for the integral in the ultraviolet (UV)
region. With the UV cutoff, the integral is convergent in 3D. Meanwhile, the
integral in 2D has an infrared (IR) divergence for $\beta=0$. The divergence
indicates that the fluctuation is so strong in 2D that the U(1) symmetry
breaking does not occur at finite temperature [6, 7]. Thus, we discuss only
(i) $\beta\neq 0$ case in 2D, and (ii) $\beta=0$ and $\beta\neq 0$ cases in
3D.

We will next discuss the evolution of the fluctuation during $t_{\rm
pump}<t<t_{0}$.

$\bullet$ Evolution of the thermal fluctuation during $t_{\rm pump}<t<t_{0}$:
In the time slot, due to the pump, $\psi_{1}$ acquires an evolution from
$\psi_{1}(\bm{r},t_{\rm pump})=-1$ to
$\psi_{1}(\bm{r},t_{0})=\sqrt{|\alpha|}\tanh(z/\xi)$ as shown in Fig. S2,S3.
The “model A” dynamics gives the following linearized equation of motion of
$\psi_{2}$.

$\displaystyle\frac{1}{\gamma}\partial_{t}\psi_{2}$
$\displaystyle=\xi_{0}^{2}\nabla^{2}\psi_{2}-2[\psi_{1}^{2}(z,t)+\alpha(z,t)+\beta]\psi_{2}+\eta_{2},$
(S23)

with

$\displaystyle\langle\eta_{2}(\bm{r},t)\eta_{2}(\bm{r^{\prime}},t^{\prime})\rangle$
$\displaystyle=\frac{2T(z,t)}{\gamma
E_{\text{c}}}\delta(\bm{r}-\bm{r^{\prime}})\delta(t-t^{\prime}).$ (S24)

Here $\psi_{1}^{2}(z,t)$ follows the solution of Eq. (S1), $\alpha(z,t)$ and
$T(z,t)$ take the form appeared in Eq. (S1). Apply the Fourier transform to
the $x$-$y$ coordinates,
$\psi_{2}(\bm{r},t)\rightarrow\tilde{\psi}(\bm{K},z,t)$,
$\eta_{2}(\bm{r},t)\rightarrow\tilde{\eta}_{2}(\bm{K},z,t)$ with the momentum
${\bm{K}}$ in the $x$-$y$ plane, we obtain the following equation

$\displaystyle\frac{1}{\gamma}\tilde{\psi}_{2}(\bm{K},z,t)=s(\bm{K},z,t)\tilde{\psi}_{2}(\bm{K},z,t)+\tilde{\eta}_{2}(\bm{K},z,t)$
(S25)

with

$\displaystyle
s(\bm{K},z,t)=-\xi_{0}^{2}K^{2}+\tilde{\psi}^{-1}_{2}(\bm{K},z,t)\xi_{0}^{2}\partial_{z}^{2}\tilde{\psi}_{2}(\bm{K},z,t)-2[\psi_{1}^{2}(z,t)+\alpha(z,t)+\beta].$
(S26)

Therefore, we have the following solution for
$\tilde{\psi}_{2}(\bm{K},t_{0})$,

$\displaystyle\tilde{\psi}_{2}(\bm{K},z,t_{0})$
$\displaystyle=e^{\gamma\int_{t_{\rm
pump}}^{t_{0}}s(\bm{K},z,t^{\prime})dt^{\prime}}\tilde{\psi}(\bm{K},z,t_{\rm
pump})+\gamma\int_{t_{\rm
pump}}^{t_{0}}e^{\gamma\int_{t^{\prime}}^{t_{0}}s(\bm{K},z,t^{\prime\prime})dt^{\prime\prime}}\tilde{\eta}_{2}(\bm{K},z,t^{\prime})dt^{\prime}$
$\displaystyle\equiv\mathcal{T}(\bm{K},z,t_{0},t_{\rm
pump})\tilde{\psi}_{2}(\bm{K},z,t_{\rm pump})+\gamma\int_{t_{\rm
pump}}^{t_{0}}\mathcal{T}(\bm{K},z,t_{0},t^{\prime})\tilde{\eta}_{2}(\bm{K},z,t^{\prime})dt^{\prime}$
(S27)

with $\mathcal{T}(\bm{K},z,t_{1},t_{2})={\rm
exp}[\gamma\int_{t_{2}}^{t_{1}}s(\bm{K},z,t)dt]$. Since the formation of the
domain wall is a fast process, e.g. $\gamma(t_{0}-t_{\rm pump})\sim 1$ (see
Fig. S3), the $\mathcal{T}$-factor gives a result of order 1 for long-
wavelength modes with $\xi_{0}^{2}K^{2}\leq 1$. Therefore, for a practical
calculation, we neglect the $\mathcal{T}$-factor in later discussion.

$\bullet$ Calculation of $D_{K}(t_{0})$: To calculate
$D_{K}(t_{0})\equiv\langle\tilde{\varphi}(-\bm{K},t_{0})\tilde{\varphi}(\bm{K},t_{0})\rangle$,
we read out $\tilde{\varphi}(\bm{K},t_{0})$ from
$\tilde{\psi}_{2}({\bm{K}},z,t_{0})$ by $\tilde{\varphi}(\bm{K},t_{0})=\int
u(z)\tilde{\psi}_{2}(\bm{K},z,t_{0})dz$ with
$u(z)\equiv(2\xi)^{-1}\text{sech}(z/\xi)$ being the mode amplified with time,
which is localized around the domain wall.

$\displaystyle\tilde{\varphi}(\bm{K},t_{0})$ $\displaystyle\approx\int
u(z)\tilde{\psi}_{2}(\bm{K},z,t_{\rm pump})dz+\int u(z)\gamma\int_{t_{\rm
pump}}^{t_{0}}\tilde{\eta}_{2}(\bm{K},z,t^{\prime})dt^{\prime}dz\equiv\tilde{\varphi}_{a}(\bm{K},t_{0})+\tilde{\varphi}_{b}(\bm{K},t_{0}).$
$\displaystyle D_{K}(t_{0})$
$\displaystyle\approx\langle\tilde{\varphi}_{a}(-\bm{K},t_{0})\tilde{\varphi}_{a}(\bm{K},t_{0})\rangle+\langle\tilde{\varphi}_{b}(-\bm{K},t_{0})\tilde{\varphi}_{b}(\bm{K},t_{0})\rangle\equiv
D_{K,a}(t_{0})+D_{K,b}(t_{0}).$ (S28)

Here $\tilde{\varphi}_{a}$ and $\tilde{\varphi}_{b}$ come from the thermal
fluctuation in the equilibrium state for $t<t_{\rm pump}$ and the thermal
fluctuation accumulated in $t_{\rm pump}<t<t_{0}$, respectively. Now we
calculate $D_{K,a}(t_{0})$ and $D_{K,b}(t_{0})$.

$\bullet\star$ Calculation of $D_{K,a}(t_{0})$:

$\displaystyle l^{D-1}D_{K,a}(t_{0})$
$\displaystyle=\int\frac{dk_{z}}{2\pi}\frac{T}{E_{\text{c}}}\frac{1}{2\beta+\xi_{0}^{2}(k_{z}^{2}+K^{2})}\\!\
|l\tilde{u}(k_{z})|^{2}$
$\displaystyle=\int\frac{dk_{z}}{2\pi}\frac{T}{E_{\text{c}}}\frac{1}{2\beta+\xi_{0}^{2}(k_{z}^{2}+K^{2})}\\!\
\Bigg{|}\int\frac{dz}{2\xi}{\rm
sech}\big{(}\frac{z}{\xi}\big{)}e^{ik_{z}z}\\!\ \Bigg{|}^{2}$
$\displaystyle=\int\frac{dk_{z}}{2\pi}\frac{T}{E_{\text{c}}}\frac{\pi^{2}}{4}\frac{\text{sech}^{2}(k_{z}\xi\pi/2)}{2\beta+\xi_{0}^{2}(k_{z}^{2}+K^{2})}$
$\displaystyle=\frac{T}{E_{\text{c}}\xi|\alpha|}\frac{\pi^{2}}{4}\frac{1}{\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}}\int\frac{dq_{z}}{2\pi}\frac{\text{sech}^{2}(q_{z}\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}\pi/2)}{1+q_{z}^{2}}$
(S29)

with $k_{z}\xi\equiv q_{z}\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}$. Due to a
factor of $e^{-2\xi_{0}^{2}\gamma tK^{2}}$ in the integrand of Eq. (S13),
$C_{1}({\bm{R}},t)$ in the limit of $\gamma t\gg 1/|\alpha|$ is dominated by
smaller $K$, and only those $D_{K}(t_{0})$ in $\xi_{0}^{2}K^{2}\ll|\alpha|$
matters in the long-time behavior of $C_{1}({\bm{R}},t)$. Besides, since
$2\beta<|\alpha|$ and $\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}\leq 1$, the
$q_{z}$-integral in Eq. (S29) always gives a result on the order of $1$. See
the estimations below).

$\displaystyle\int\frac{\text{sech}^{2}(q_{z})}{1+q_{z}^{2}}dq_{z}\approx
1.41\ \ <\ \
\int\frac{\text{sech}^{2}(q_{z}\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}\pi/2)}{1+q_{z}^{2}}dq_{z}\
\ <\ \ \int\frac{1}{1+q_{z}^{2}}dq_{z}=\pi.$ (S30)

We approximate it to be $\pi$ for convenience and obtain

$\displaystyle
l^{D-1}D_{K,a}(t_{0})\approx\frac{T}{E_{\text{c}}\xi|\alpha|}\frac{\pi^{2}}{8}\frac{1}{\sqrt{(2\beta+\xi_{0}^{2}K^{2})/|\alpha|}}.$
(S31)

$\bullet\star$ Calculation of $D_{K,b}(t_{0})$:

$\displaystyle l^{D-1}D_{K,b}(t_{0})$ $\displaystyle=\gamma\int_{t_{\rm
pump}}^{t_{0}}dt\ \ \int dz\ \ \frac{2T(z,t)}{E_{\rm
c}}u^{2}(z)\stackrel{{\scriptstyle z_{\rm
p}\gg\xi}}{{\approx}}\gamma\int_{t_{\rm pump}}^{t_{0}}dt\ \ \int dz\ \
\frac{2T(0,t)}{E_{\rm c}}u^{2}(z)=\gamma\int_{t_{\rm pump}}^{t_{0}}dt\ \
\frac{T(0,t)}{E_{\rm c}\xi}\equiv D_{b}.$ (S32)

The calculation gives a constant result independent of $K$.

Substituting $D_{K}(t_{0})=D_{K,a}(t_{0})+D_{K,b}(t_{0})$ in Eq. (S13), we
obtain $C_{1}({\bm{R}},t)=C_{1,a}(\bm{R},t)+C_{1,b}(\bm{R},t)$ with

$\displaystyle C_{1,a}(\bm{R},t)$
$\displaystyle=\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K,a}(t_{0})e^{2(|\alpha|-2\beta)\gamma
t}e^{-2\xi_{0}^{2}\gamma tK^{2}}e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=\frac{T}{E_{\text{c}}\xi_{0}}e^{2(|\alpha|-2\beta)\gamma
t}\frac{\pi^{2}}{8}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}\frac{1}{\sqrt{(2\beta+\xi_{0}^{2}K^{2})}}e^{-2\xi_{0}^{2}\gamma
tK^{2}}e^{i\bm{K}\cdot\bm{R}}.$ (S33)

$\displaystyle C_{1,b}(\bm{R},t)$
$\displaystyle=\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K,b}(t_{0})e^{2(|\alpha|-2\beta)\gamma
t}e^{-2\xi_{0}^{2}\gamma tK^{2}}e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=D_{b}e^{2(|\alpha|-2\beta)\gamma
t}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}e^{-2\xi_{0}^{2}\gamma
tK^{2}}e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=\frac{D_{b}|\alpha|^{D/2-2}}{\sqrt{|\alpha|}\xi_{0}^{D-1}}\frac{|\alpha|}{\sqrt{8\pi|\alpha|\gamma
t}^{D-1}}e^{2(|\alpha|-2\beta)\gamma t}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma
t}}.$ (S34)

By the definition of $D_{b}$, the factor
$\frac{D_{b}|\alpha|^{D/2-2}}{\sqrt{|\alpha|}\xi_{0}^{D-1}}\sim\zeta$. The
integration of $C_{1,a}(\bm{R},t)$ could be evaluated in some relevant limits.
In the following, we discuss $C_{1,a}(\bm{R},t)$ in these limits.

$\bullet$ $\gamma t\gg 1/\beta$ limit: In this limit, the $K$ integral is
dominated by those $K$ in $\xi_{0}^{2}K^{2}\ll\beta$ due to a factor of
$e^{-2\xi_{0}^{2}\gamma tK^{2}}$ in the integrand of Eq. (S33). Thus, we can
estimate the integral as follows

$\displaystyle C_{1,a}(\bm{R},t)$
$\displaystyle\approx\frac{T}{E_{\text{c}}\xi_{0}}e^{2(|\alpha|-2\beta)\gamma
t}\frac{\pi^{2}}{8}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}\frac{1}{\sqrt{2\beta}}e^{-2\xi_{0}^{2}\gamma
tK^{2}}e^{i\bm{K}\cdot\bm{R}}$ $\displaystyle=Ge^{2(|\alpha|-2\beta)\gamma
t}\frac{\pi^{2}}{8}\frac{1}{\sqrt{2\beta}}\frac{1}{\sqrt{8\pi\gamma
t}^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma t}}$
$\displaystyle\approx\zeta|\alpha|\frac{\sqrt{|\alpha|/\beta}}{\sqrt{8\pi|\alpha|\gamma
t}^{D-1}}e^{2(|\alpha|-2\beta)\gamma t}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma
t}}.$ (S35)

Combining with Eqs. (S15, S34), we obtain Eq. (6) in the main text.

$\bullet$ U(1) symmetric case in the 3D system ($\beta=0$ and $D=3$): The 3D
${\bm{K}}$-integral in Eq. (S33) can be exactly evaluated in the polar
coordinate with the Bessel functions,

$\displaystyle C_{1,a}(\bm{R},t)$
$\displaystyle=\frac{T}{E_{\text{c}}\xi_{0}^{2}}e^{2|\alpha|\gamma
t}\frac{\pi^{2}}{8}\int_{0}^{+\infty}\frac{dK}{2\pi}e^{-2\xi_{0}^{2}\gamma
tK^{2}}J_{0}(KR)=\frac{T}{E_{\text{c}}\xi_{0}^{2}}e^{2|\alpha|\gamma
t}\frac{\pi^{2}}{16}\frac{1}{\sqrt{8\pi\xi_{0}^{2}\gamma
t}}I_{0}(\frac{R^{2}}{16\xi_{0}^{2}\gamma
t})e^{-\frac{R^{2}}{16\xi_{0}^{2}\gamma t}}$ $\displaystyle\approx
Ge^{2|\alpha|\gamma t}\frac{1}{\sqrt{8\pi\gamma
t}}I_{0}(\frac{R^{2}}{16\xi_{0}^{2}\gamma
t})e^{-\frac{R^{2}}{16\xi_{0}^{2}\gamma t}}=\zeta|\alpha|e^{2|\alpha|\gamma
t}\frac{1}{\sqrt{8\pi|\alpha|\gamma t}}I_{0}(\frac{R^{2}}{16\xi_{0}^{2}\gamma
t})e^{-\frac{R^{2}}{16\xi_{0}^{2}\gamma t}}.$ (S36)

Here $J_{0}(x)$ is the Bessel function of the first kind, and $I_{0}(x)$ is
the modified Bessel function of the first kind. A comparison between Eq. (S36)
and Eq. (S15) with $D=3$ and $\beta=0$ shows that $C_{1,a}({\bm{R}},t)$ always
dominates over $C_{1,b}({\bm{R}},t)$ and $C_{2}({\bm{R}},t)$ in the long-time
limit, leading to Eq. (4) in the main text. From Eq. (S36),
$\lim_{x\rightarrow 0}I_{0}(x)=1$, and
$\lim_{x\rightarrow+\infty}=e^{x}/\sqrt{2\pi x}$, we can see that the spatial
dependence of the correlation in the long-time limit ($\gamma t\gg 1$) shows a
crossover from an exponential decay to a power-law decay;

$\displaystyle
C_{1,a}({\bm{R}},t)\propto\begin{cases}e^{-R^{2}/16\xi_{0}^{2}\gamma
t},&\text{for}\\!\ \\!\ R\ll 4\xi_{0}\sqrt{\gamma t}\\\
\frac{4\xi_{0}\sqrt{\gamma t}}{R},&\text{for}\\!\ \\!\ R\gg
4\xi_{0}\sqrt{\gamma t}.\end{cases}$ (S37)

The crossover is consistent with the power-law spatial correlation
$\langle\psi_{2}(\bm{r})\psi_{2}(\bm{r^{\prime}})\rangle\sim
1/|\bm{r}-\bm{r^{\prime}}|$ at long distances in the initial equilibrium
fluctuation in the $3D$ systems.

###  II.3 Determine the crossover time $t_{\text{c}}$.

When $|\psi_{2}({\bm{r}},t)|^{2}=C(0,t){\rm sech}^{2}(z/\xi)$ becomes on the
same order as $|\psi_{1}({\bm{r}},t)|^{2}=|\alpha|\tanh^{2}(z/\xi)$, the
first-stage dynamics crossover into the second-stage dynamics. The crossover
time scale $t_{c}$ is defined as
$C({\bm{0}},t_{c})=C_{1,a}({\bm{0}},t)+C_{1,b}({\bm{0}},t)+C_{2}({\bm{0}},t)\sim
C_{1,a}({\bm{0}},t)+2C_{2}({\bm{0}},t)=|\alpha|$. Note that the integrals in
Eq. (S33) can be calculated exactly at $\bm{R}=0$ for both 2D and 3D cases,

$\displaystyle
C_{1,a}(0,t)=\begin{cases}\zeta|\alpha|e^{2|\alpha|\gamma(t-t_{0})}\frac{\pi\sqrt{\pi}}{32}\frac{1}{\sqrt{2|\alpha|\gamma(t-t_{0})}}(1-\text{erf}(\sqrt{4\beta\gamma(t-t_{0})})),&\text{in
3D systems,}\\\
\zeta|\alpha|e^{2(|\alpha|-\beta)\gamma(t-t_{0})}\frac{\pi}{16}K_{0}(2\beta\gamma(t-t_{0})),&\text{in
2D systems.}\end{cases}$ (S38)

Here $\text{erf}(x)$ and $K_{0}(x)$ are the error function and the modified
Bessel function of the second kind, respectively. For the long-time limit
($\gamma|t-t_{0}|\gg 1/(|\alpha|-2\beta)$), $C_{2}(0,t)$ is evaluated in Eq.
(S15).

In the general case, we can use Eq. (S15,S38) to calculate $C(0,t_{\text{c}})$
and determine $t_{\text{c}}$ by $C(0,t_{\text{c}})=|\alpha|$. The results of
$\gamma t_{\text{c}}$ mentioned in the main text are approximated results in
the limits discussed.

##  III The second-stage dynamics

When $t$ reaches the crossover time $t_{c}$, randomly distributed
positive-$\psi_{2}$ and negative-$\psi_{2}$ domains are located on the $z=0$
interface. Since $\psi$ is a continuous function of $\bm{r}$, the boundary
between the positive-$\psi_{2}$ and negative-$\psi_{2}$ domains is nothing but
the zero of $\psi_{2}(\bm{R},z=0)$. Furthermore, $\psi_{1}(\bm{R},z=0)$
remains to be $0$. As a result, $|\psi|=0$ on the boundary of the domains on
the $z=0$ interface, which suggests the existence of topological defects.
Writing $\psi=|\psi|e^{i\theta}$, the profile of $\theta$ around these defects
could be analyzed as follows. On the interface, $\theta(\bm{R},z=0)$ takes the
value of $\pi/2$ and $3\pi/2$ in the positive-$\psi_{2}$ and
negative-$\psi_{2}$ domains. Far away from the interface, the boundary
condition at $z=\pm\infty$ set $\psi_{1}=\pm\sqrt{|\alpha|}$ and $\psi_{2}=0$,
i.e. $\theta=0/\pi$ at $z=+\infty/-\infty$.

Based on the analysis above, on the 1D interface of 2D systems, a clockwise
closed path around a zero of $\psi_{2}$ obtains a phase winding
$2\pi$($-2\pi$) if the slope $\partial_{y}\psi_{2}(y,z=0)$ is
negative(positive). Therefore, a zero of $\psi_{2}$ with a negative(positive)
$\partial_{y}\psi_{2}$ is nothing but a vortex(antivortex). Since the
positive-$\psi_{2}$ and negative-$\psi_{2}$ domains appear alternatively along
the $y$ axis, the vortex and antivortex also appear alternatively. In other
words, a vortex(antivortex) has to lie in between two neighboring
antivortex(vortex). See Fig. (S4)(a) for a schematic plot of the
$\psi_{2}(y,z=0)$ function, the $\theta(y,z)$ profile and the
vortex/antivortex. On the 2D interface of 3D systems, the zeros of $\psi_{2}$
form 1D lines called vortex strings. The orientations of the vortex strings
are defined by a right-hand rule with respect to the direction of
$\nabla\theta$ around the strings, as shown in Fig. S4(b). In any intersecting
2D plane parallel to $z$, $\theta$ distributes in a similar way as in the 2D
systems, and vortex/antivortex appears alternatively on the $z=0$ interface in
the plane. Since the positive-$\psi_{2}$ and negative-$\psi_{2}$ domains are
closed, the line of zeros of $\psi_{2}$, which is nothing but the boundary
between the domains, are closed loops. Therefore, the vortex strings form
loops.

In the second stage dynamics of 2D systems, the vortex moves under an
attractive force from neighboring antivortices, a frictional force from the
damping term $\gamma^{-1}\partial_{t}\psi$ and Langevin force from the thermal
noise $\eta$ in the “model A” dynamics. Due to the attractive force, the
neighboring vortex and antivortex move toward each other and annihilate. The
annihilation processes lead to a coarsening phenomenon of the topological
defects in the interface [8]. In the second stage dynamics of 3D systems, the
vortex string contracts under a tension corresponding to the local curvature,
a frictional force from the damping term $\gamma^{-1}\partial_{t}\psi$ and
Langevin force from the thermal noise $\eta$ in the “model A” dynamics. It
contracts and annihilates, leading to a coarsening phenomenon. In this
section, we derive an equation of motion of the defects, and analyze
quantitatively the coarsening process in the second stage dynamics.

Figure S4:  Schematic plot of the profile of $\psi_{2}$ after the first-stage
dynamics on the $z=0$ interface. (a) On the 1D interface of 2D systems, the
nodes of $\psi_{2}$ generate a vortex-antivortex-pair array denoted by
$\otimes/\odot$. The inset is the corresponding profile of $\theta$. (b) On
the 2D interface of 3D systems, $\psi_{2}$ are positive/negative in the
black/white regions. The borders between the black and white regions, where
$\psi=0$, generate vortex strings. The vortex strings have specific
orientations marked by the blue arrows. The orientations are defined by a
right-hand rule concerning the direction of $\nabla\theta$ around the strings,
illustrated in the inset.

###  III.1 Spatial profile of the topological defect without the $U(1)$
symmetry ($\beta\neq 0$).

We first determine a spatial profile of the defect to derive the equation of
motion of a topological defect from the “model A” dynamics. For the
$U(1)$-symmetric case ($\beta=0$), the spatial profile of the vortex or vortex
string is well known [8]. In this section, we will determine the spatial
profile of the vortex without the $U(1)$ symmetry ($\beta\neq 0$). Especially,
we will clarify the spatial profile far from the vortex core, and the profile
near the core. In the following, we will focus on a 2D spatial profile of a
vortex pinned in the 1D interface of the 2D systems. A vortex loop in the 3D
systems can be regarded as a stack of the 2D spatial profile of the vortex
along the third direction. We also limit ourselves to the ‘nearly U(1)
symmetric’ case ($\beta\ll|\alpha|$) for analytical convenience.

The spatial profile of a single vortex forms a local minimum of the Ginzburg-
Landau free energy functional of
$\psi({\bm{r}})=|\psi(\bm{r})|e^{i\theta({\bm{r}})}$, i.e. Eq. (S3). In the
limit of $\beta\ll|\alpha|$, it is a good approximation to neglect the
amplitude dynamics and fix the amplitude $\psi=\sqrt{|\alpha|}$, and study the
free energy functional of $\theta({\bm{r}})$ with ${\bm{r}}\equiv(y,z)$;

$\displaystyle F$ $\displaystyle=E_{\text{c}}|\alpha|\int d^{2}\bm{r}\ \
\xi_{0}^{2}|\nabla\theta|^{2}+2\beta\sin^{2}\theta.$ (S39)

Or equivalently, $\delta F/\delta\theta=0$;

$\displaystyle\xi_{0}^{2}\nabla^{2}\theta=\beta\sin(2\theta).$ (S40)

With the boundary condition $\lim_{z\rightarrow+\infty}\theta_{0}=0$,
$\lim_{z\rightarrow-\infty}\theta_{0}=\pi$, the equation leads to the
following two solutions,

$\displaystyle\theta_{0}^{\pm}(z)=\text{arg}(\sinh(\sqrt{2\beta}z/\xi_{0})\pm
i)$ (S41)

These solutions are nothing but the Sine-Gordon solitons. The soliton has a
finite energy per unit length along the $y$-direction,

$\displaystyle\sigma_{0}=E_{\text{c}}|\alpha|\int dz\ \
\xi_{0}^{2}|\nabla\theta_{0}^{\pm}|^{2}+2\beta\sin^{2}\theta_{0}^{\pm}=4E_{\text{c}}\xi_{0}|\alpha|\sqrt{2\beta}.$
(S42)

In the second-stage dynamics, $\theta(y,z)=\theta^{+}_{0}(z)$ and
$\theta(y,z)=\theta^{-}_{0}(z)$ appears alternatively along the $y$-direction,
connected by vortices. Consider a vortex at $(y,z)=(0,0)$. The 2D spatial
profile of the vortex can be described by combining the two soliton solutions
at $y=0$; $\theta(y,z)=\theta_{0}^{+}(z)$ for $y\rightarrow+\infty$ and
$\theta(y,z)=\theta_{0}^{-}(z)$ for $y\rightarrow-\infty$ [9]. Such spatial
profile of the vortex is generally given by

$\displaystyle\theta(y,z)=\begin{cases}\theta^{+}_{0}(z)+s(y,z)&{\rm for}\ \
\\!\ y>0,\\\ \theta^{-}_{0}(z)+s(y,z)&{\rm for}\ \ \\!\ y<0.\\\ \end{cases}$
(S43)

Note that at $y$ approaches $\pm\infty$, the vortex profile approaches the
Sine-Gordon soliton profile due to small but finite $\beta$. At $|z|$
approaches $+\infty(-\infty)$, $\theta=0(\pi)$. Therefore, $s$ in Eq. (S43)
satisfies the boundary condition $\lim_{y\rightarrow\pm\infty}s(\bm{r})=0$.
$s$ and $\lim_{z\rightarrow\pm\infty}s(\bm{r})=0$. Thus, we can determine a
form of $s(y,z)$ for the larger $|y|$ by substituting
$\theta(\bm{r})=\theta_{0}(z)+s(\bm{r})$ into Eq. (S40) and linearizing the
equation in small $s$;

$\displaystyle\xi_{0}^{2}\nabla^{2}s=2\beta\cos(2\theta_{0})s.$ (S44)

The linearized equation is decomposed into an operator in $y$ and an operator
in $z$,

$\displaystyle\frac{\xi_{0}^{2}}{2\beta}\partial_{y}^{2}s(y,z)=[-\frac{\xi_{0}^{2}}{2\beta}\partial_{z}^{2}-\frac{2}{\cosh^{2}(\sqrt{2\beta}z/\xi_{0})}+1]s(y,z)\equiv[H(z^{\prime};2)+1]s(y,z),$
(S45)

with $z^{\prime}\equiv\sqrt{2\beta}z/\xi_{0}$. $H(z^{\prime};\lambda)$ is
nothing but the 1D Hamiltonian with the Poschl-Teller potential hole mentioned
in Eq. (S5). Thus, in terms of the eigenmodes of $H(z^{\prime};2)$, $s(y,z)$
can be expanded as

$\displaystyle s(y,z)=c_{0}h_{0}(z^{\prime})+\int dk^{\prime}\\!\
c(k^{\prime})\\!\
e^{-\sqrt{{k^{\prime}}^{2}+1}|y^{\prime}|}h(k^{\prime},z^{\prime}).$ (S46)

with $k^{\prime}\equiv\xi_{0}k/\sqrt{2\beta}$,
$y^{\prime}=y\sqrt{2\beta}/\xi_{0}$, $z^{\prime}=z\sqrt{2\beta}/\xi_{0}$, and
$c^{*}(k^{\prime})=c(-k^{\prime})$. Here $h_{0}(z^{\prime})\propto{\rm
sech}(z^{\prime})$ is the bound-state eigenmode of $H(z^{\prime},2)$ with
negative eigenenergy $-1$.
$h(k^{\prime},z^{\prime})\propto\big{(}1+i/(k^{\prime})\tanh(z^{\prime})\big{)}e^{ik^{\prime}z^{\prime}}$
is a continuum eigenstate of $H(z^{\prime},2)$ with positive eigenenergy
$(k^{\prime})^{2}$. From the boundary condition at $y=\pm\infty$, $c_{0}=0$.
$c(k^{\prime})$ is chosen to satisfy the boundary condition
$\lim_{z\rightarrow\pm\infty}s=0$.

The energy for the Sine-Gordon soliton with a single vortex is given by

$\displaystyle E(R_{\rm d})=\int^{R_{\rm d}/2}_{-R_{\rm d}/2}dy\\!\
\sigma(y),$

with

$\displaystyle\sigma(y)=E_{c}|\alpha|\int dz\ \
\xi_{0}^{2}|\nabla(\theta_{0}^{\pm}+s)|^{2}+2\beta\sin^{2}(\theta_{0}^{\pm}+s).$

Here $R_{\rm d}$ is a size of the vortex along the 1D interface, and
$\theta^{+}_{0}$ and $\theta^{-}_{0}$ in the right hand side are for an energy
density per length, $\sigma(y)$, in the $y>0$ and $y<0$ regions, respectively.
Since the size of the vortex is constrained by its neighboring antivortex,
whose 2D spatial profile takes the same form as Eq. (S43) with
$y\rightarrow-y$, $R_{\rm d}$ can be also regarded as a distance between the
two defects along the 1D interface.

For large $|y|$, the energy density per length, $\sigma(y)$, can be further
evaluated by an expansion in small $s(y,z)$. The zeroth order term in $s$ is
the energy density of the domain wall without the vortices, $\sigma_{0}$. The
first-order term in $s$ vanishes because the soliton solution
$\theta^{\pm}_{0}$ is a local energy minimum of the free energy. The second-
order term gives the energy density of a vortex per length for large $|y|$;

$\displaystyle\sigma(y)$
$\displaystyle\approx\sigma_{0}+E_{\text{c}}|\alpha|\int dz\ \
\xi_{0}^{2}(\partial_{y}s)^{2}+\xi_{0}^{2}(\partial_{z}s)^{2}+2s^{2}\beta\cos(2\theta_{0})$
$\displaystyle=\sigma_{0}+E_{\text{c}}|\alpha|\int dz\ \
\xi_{0}^{2}(\partial_{y}s)^{2}+\xi_{0}^{2}\partial_{z}(s\partial_{z}s)-\xi_{0}^{2}(s\partial_{z}^{2}s)+2s^{2}\beta\cos(2\theta_{0})$
$\displaystyle=\sigma_{0}+E_{\text{c}}|\alpha|\int dz\ \
\xi_{0}^{2}(\partial_{y}s)^{2}+\xi_{0}^{2}(s\partial_{y}^{2}s)+\xi_{0}^{2}\partial_{z}(s\partial_{z}s)$
$\displaystyle=\sigma_{0}+2E_{\rm
c}|\alpha|\xi_{0}\sqrt{2\beta}\int|c(k^{\prime})|^{2}(k^{\prime
2}+1)e^{-2\sqrt{k^{\prime 2}+2}|y^{\prime}|}\
dk^{\prime}+E_{\text{c}}|\alpha|\xi_{0}^{2}(s\partial_{z}s|^{z=+\infty}_{z=-\infty})$
$\displaystyle=\sigma_{0}+2E_{\rm
c}|\alpha|\xi_{0}\sqrt{2\beta}\int|c(k^{\prime})|^{2}(k^{\prime
2}+1)e^{-2\sqrt{k^{\prime 2}+2}|y|\sqrt{2\beta}/\xi_{0}}\ dk^{\prime}$ (S47)

From the second line to the third line, we use Eq. (S44). From the third line
to the fourth line, we use Eq. (S46) and an orthonormal condition of
$h(k^{\prime},z^{\prime})$; $\int
dz^{\prime}h(k^{\prime}_{1},z^{\prime})h^{*}(k^{\prime}_{2},z^{\prime})=\delta(k^{\prime}_{1}-k^{\prime}_{2})$.
From the fourth line to the fifth line, we use the boundary condition
$\lim_{z\rightarrow\pm\infty}s=0$.

The free energy can be decomposed into an energy of the vortex, $E_{d}(R_{\rm
d})$, and the energy of the soliton, $\sigma_{0}R_{\rm d}$,

$\displaystyle E(R_{\rm d})=\sigma_{0}R_{\rm d}+\int^{R_{\rm d}/2}_{-R_{\rm
d}/2}dy\ \ [\sigma(y)-\sigma_{0}]\equiv\sigma_{0}R_{\rm d}+E_{d}(R_{\rm d}).$
(S48)

By Eq. (S47), $\sigma(y)-\sigma_{0}$ decays exponentially for larger $|y|$.
Therefore, the energy of the vortex is convergent in the larger $R_{\rm d}$
limit; $\lim_{R_{\rm d}\rightarrow+\infty}E_{\text{d}}(R_{\rm d})=E_{0}\propto
E_{\text{c}}|\alpha|\xi_{0}^{2}$.

Since $R_{\rm d}$ can also be regarded as the distance between vortex and
antivortex along the 1D interface, $dE_{\rm d}/dR_{\rm d}$ gives an estimate
of an action-reaction force between the two defects. The force is attractive
and it is dependent on the distance $R_{\rm d}$. For larger $R_{\rm d}$ with
$R_{\rm d}\gg\xi_{0}/\sqrt{2\beta}$, the magnitude of the force decays
exponentially,

$\displaystyle F_{\rm d}=\frac{dE_{\text{d}}}{dR_{\rm d}}=\sigma(R_{\rm
d}/2)\sim E_{\rm c}|\alpha|\xi_{0}\sqrt{2\beta}e^{-R_{\rm
d}\sqrt{2\beta}/\xi_{0}}.$ (S49)

We have discussed the energy and attraction force of the vortex far away from
the vortex core, $R_{\rm d}\gg\xi_{0}/\sqrt{2\beta}$. Next, we will discuss
the energy and force of the vortex near the vortex core $R_{\rm
d}\ll\xi_{0}/\sqrt{2\beta}$. Let us rewrite Eq. (S39) in the polar
coordinates,

$\displaystyle F$ $\displaystyle=E_{\text{c}}|\alpha|\int^{R_{\rm
d}/2}_{0}rdr\int^{2\pi}_{0}d\phi\ \
\Big{[}\xi^{2}_{0}(\frac{\partial\theta}{\partial
r})^{2}+\frac{\xi^{2}_{0}}{r^{2}}(\frac{\partial\theta}{\partial\phi})^{2}+2\beta\sin^{2}\theta\Big{]}\equiv
E_{\text{c}}|\alpha|\int^{R_{\rm d}/2}_{0}dr\sigma(r).$ (S50)

For $r\ll\xi_{0}/\sqrt{2\beta}$, the vortex profile approaches the U(1)
symmetric vortex, where $\sigma(r)$ is dominated by the gradient term;

$\displaystyle\sigma(r)\equiv\int^{2\pi}_{0}d\phi\ \
\frac{\xi^{2}_{0}}{r}\Big{[}r^{2}(\frac{\partial\theta}{\partial
r})^{2}+(\frac{\partial\theta}{\partial\phi})^{2}+\frac{2\beta
r^{2}}{\xi^{2}_{0}}\sin^{2}\theta\Big{]}\approx\int^{2\pi}_{0}d\phi\ \
\frac{\xi^{2}_{0}}{r}\Big{[}r^{2}(\frac{\partial\theta}{\partial
r})^{2}+(\frac{\partial\theta}{\partial\phi})^{2}\Big{]}.$ (S51)

Thus, for $R_{\rm d}\ll\xi_{0}/\sqrt{2\beta}$, $E_{d}(R_{\rm d})\approx
E_{\text{c}}|\alpha|\int d^{2}\bm{r}\\!\
\xi_{0}^{2}|\nabla\theta|^{2}\simeq|\alpha|E_{\text{c}}\xi_{0}^{2}\ln(R_{\rm
d}/\xi)$. Here the coherence length $\xi$ plays the role of the core size of
the vortices. Thereby, the attractive force decays with a power law for
$R_{\rm d}\ll\xi_{0}/\sqrt{2\beta}$;

$\displaystyle F_{\text{d}}(R_{\rm d})=\frac{dE_{\rm d}}{dR_{\rm
d}}\sim\frac{\xi_{0}^{2}|\alpha|E_{\text{c}}}{R_{\rm d}},$ (S52)

which is the same as the $U(1)$-symmetric case.

The 2D spatial profile of the vortex discussed so far has a characteristic
length scale, $\xi_{0}/\sqrt{2\beta}$, which separates the long-length regime
with the exponential decay of the attractive force, and the short-length
regime with the power law of the force. We confirmed that numerical solutions
of the vortex are consistent with these analytical results in the two limits;
Fig. S5, verifying the conclusion above. In Fig. S5, we can see that the
defect behaves like a circular symmetric vortex in the shorter length regions
$|y|\ll\xi_{0}/\sqrt{2\beta}$, while in the long-length region
$|y|\gg\xi_{0}/\sqrt{2\beta}$, it behaves as a Sine-Gordon soliton with a
translation symmetry along the $y$-direction.

In the 2D interface of the 3D systems, the vortex forms a vortex string, and
the cross-section of the string contains the vortices discussed above. For a
vortex-string loop with diameter $\sim R_{\text{d}}$, its energy is estimated
as $\sim R_{\text{d}}E_{\rm d}(R_{\rm d})$ where $E_{\rm d}(R_{\rm d})$ takes
the expression obtained above.

Figure S5: Numerical solution of a vortex of the free energy Eq. (S3) for the
2D spatial profile of the phase with a vortex core at $(y,z)=(0,0)$. The
numerical solution is obtained with $|\alpha|=-1$, $\beta=0.02$. i.e.
$\xi_{0}/\sqrt{2\beta}=5\xi_{0}$.

###  III.2 Force balance of the defects

In this section, we start from the 2D systems and derive an equation of motion
for a single vortex from the “model A dynamics”. The equation of motion
becomes a force balance equation among the attractive force between vortex and
antivortex, frictional and Langevin forces. When a pair of vortex and
antivortex meet each other due to the attractive force, they annihilate.
Events of the pair annihilations keep enlarging the mean distance between
vortex and antivortex in the whole interface, leading to a coarsening
phenomenon [10, 11, 12]. From the force balance equation, we derive an
equation of motion of the mean distance and estimate a growth law of the mean
distance (correlation length) in the coarsening process. In 3D systems, a
similar thing happens and the topological defects are vortex strings. They
contract and annihilate driven by the tension that stems from the local
curvature [13, 14], leading to a coarsening phenomenon. The analytical methods
used in this section are mainly adopted from a review paper by Bray [8].

####  III.2.1 2D systems

Figure S6: A schematic picture of the force balance relation for a vortex
moving with a velocity $\bm{v}$ around the interface at $z=0$. An attractive
force $\bm{F}_{\text{d}}$ is along the grey dashed line, which connects the
vortex ($\otimes$) and its neighboring antivortex $(\odot)$ with a distance
$R_{\rm d}$. $\delta_{z}$ is a displacement of the vortex from $z=0$. We
consider the $\delta_{z}/R_{\rm d}\rightarrow 0$ limit. Thus, the direction of
$\bm{F}_{\rm d}$ is nearly along the interface.

Let us begin with the 2D system. Consider that a vortex core moves with a
velocity $\bm{v}(t)$, represented by the $\otimes$ symbol in Fig. (S6). The
spatial profile of such a moving vortex is given by
$\psi(\bm{r},t)=f(\bm{r}-\int^{t}\bm{v}(t^{\prime})dt^{\prime})$, where
$f(\bm{r})$ is given by a defect at ${\bm{r}}=(0,0)$ in Eq. (S43);
$f({\bm{r}})\equiv\sqrt{|\alpha|}{\rm exp}[i\theta({\bm{r}})]$.
$\psi(\bm{r},t)$ satisfies the dynamical equation,

$\displaystyle\frac{1}{\gamma^{2}_{0}}\partial^{2}_{t}\psi+\frac{1}{\gamma}\partial_{t}\psi=-\frac{\delta
F}{\delta\psi^{*}}+\eta,$ (S53)

where $F[\psi]$ takes the form of Eq. (S3). Here for the usage of the later
section, we also added the inertia term $\gamma_{0}^{-2}\partial^{2}_{t}\psi$
based on the “model A” dynamics in Eq. (S2) and derived a general equation of
motion for the vortex.

For a moving vortex with
$\psi(\bm{r},t)=f(\bm{r}-\int^{t}\bm{v}(t^{\prime})dt^{\prime})$, its distance
from the neighboring antivortex could be regarded as its size $R_{\rm d}$. In
the second-stage dynamics, the vortex moves toward the neighboring antivortex
under the attraction force and $R_{\rm d}$ gets smaller with time going by.
Thus, its energy $E_{\rm d}$ changes with the following energy-dissipation
rate

$\displaystyle\frac{dE_{\rm d}}{dt}=\frac{dE_{\rm d}}{dR_{\rm d}}\frac{dR_{\rm
d}}{dt}=-\bm{F}_{\rm d}\cdot\bm{v}.$ (S54)

Here $\bm{F}_{\rm d}$ is the attraction force from the neighboring antivortex
with a magnitude $F_{\rm d}=dE_{\rm d}/dR_{\rm d}$, as illustrated in Fig. S6.
We consider the vortices are nearly confined in the interface, i.e. the
distance between the vortices are much larger than their displacements from
the $z=0$ interface, which is illustrated by the $\delta_{z}/R_{\rm
d}\rightarrow 0$ limit in Fig. S6. Therefore, $\bm{F}_{\rm d}$ is along the
$y$-direction, i.e. ${F}_{{\rm d},z}=0$ and $F_{{\rm d},y}=F_{\rm d}$. Since
$\bm{v}$ is measured from the frame of the neighboring antivortex,
$v_{y}=-dR_{\rm d}/dt$ and $-\bm{F}_{\rm d}\cdot\bm{v}=-F_{\rm
d,y}v_{y}=F_{\rm d}dR_{\rm d}/dt$, yielding Eq. (S54). The asymptotic forms of
$F_{\rm d}$ in the two limits, $R_{\rm d}\gg\xi_{0}/\sqrt{2\beta}$ and $R_{\rm
d}\ll\xi_{0}/\sqrt{2\beta}$, are discussed in Eqs. (S49,S52).

According to Eq. (S53), the energy-dissipation rate can also be given by the
inertia term of the vortex, frictional force, and Langevin force acting on the
vortex,

$\displaystyle\frac{dE_{\text{d}}}{dt}$ $\displaystyle=\int d^{2}\bm{r}\ \
\frac{\delta F}{\delta\psi}\frac{\partial\psi}{\partial t}+\frac{\delta
F}{\delta\psi^{*}}\frac{\partial\psi^{*}}{\partial t}$ $\displaystyle=\int
d^{2}\bm{r}\ \ \Big{(}-\frac{E_{\rm
c}}{\gamma^{2}_{0}}\frac{\partial^{2}\psi^{*}}{\partial
t^{2}}-\frac{E_{\text{c}}}{\gamma}\frac{\partial\psi^{*}}{\partial
t}+E_{\text{c}}\eta^{*}\Big{)}\\!\ \frac{\partial\psi}{\partial
t}+\Big{(}-\frac{E_{\rm c}}{\gamma^{2}_{0}}\frac{\partial^{2}\psi}{\partial
t^{2}}-\frac{E_{\text{c}}}{\gamma}\frac{\partial\psi}{\partial
t}+E_{\text{c}}\eta\Big{)}\\!\ \frac{\partial\psi^{*}}{\partial t}$
$\displaystyle=-\sum_{i,j,k=y,z}\int d^{2}\bm{r}\ \
\Big{(}-\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}\partial_{k}f^{*}){\bm{v}}_{i}{\bm{v}}_{k}+\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}f^{*})\frac{\partial{\bm{v}}_{i}}{\partial
t}+\frac{E_{\text{c}}}{\gamma}(\partial_{i}f^{*})\bm{v}_{i}+E_{\text{c}}\eta^{*}\Big{)}\\!\
(\partial_{j}f)\bm{v}_{j}$ $\displaystyle\hskip
85.35826pt+\Big{(}-\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}\partial_{k}f){\bm{v}}_{i}{\bm{v}}_{k}+\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}f)\frac{\partial{\bm{v}}_{i}}{\partial
t}+\frac{E_{\text{c}}}{\gamma}(\partial_{i}f)\bm{v}_{i}+E_{\text{c}}\eta\Big{)}\\!\
(\partial_{j}f^{*})\bm{v}_{j}$
$\displaystyle\equiv(\bm{F}_{\lambda}-\bm{F}_{\eta}-{\bm{A}})\cdot\bm{v}.$
(S55)

Here ${\bm{A}}=\sum_{i}m_{i}a_{i}\hat{i}$ with ${\bm{a}}=d{\bm{v}}/dt$, is the
inertia term of the vortex [15],
$\bm{F}_{\lambda}=-\sum_{i}\lambda_{i}v_{i}\hat{i}$ is the frictional force,
and $\bm{F}_{\eta}$ is the Langevin force. $m_{i}$ and $\lambda_{i}$ are the
inertia mass and the friction coefficient in the $i$-direction, respectively.
The inertia mass and friction coefficients and the Langevin force are given as
follows;

$\displaystyle m_{i}$ $\displaystyle=\frac{E_{\text{c}}}{\gamma^{2}_{0}}\int
d^{2}\bm{r}\ \ 2|\partial_{i}f|^{2},$ (S56) $\displaystyle\lambda_{i}$
$\displaystyle=\frac{E_{\text{c}}}{\gamma}\int d^{2}\bm{r}\ \
2|\partial_{i}f|^{2},$ (S57) $\displaystyle\bm{F}_{\eta}$
$\displaystyle=E_{\text{c}}\int d^{2}\bm{r}\ \ (\eta^{*}\nabla f+\eta\nabla
f^{*}).$ (S58)

From the third line to the fourth line in Eq. (S55), we set $\int
d^{2}{\bm{r}}(\partial_{i}\partial_{k}f^{*})(\partial_{j}f)+{\rm c.c.}=0$,
because $f({\bm{r}})=-f(-{\bm{r}})$.

Equating Eq. (S55) with Eq. (S54), we obtain a general equation of motion for
the vortex under the force from its neighboring antivortex,

$\displaystyle m_{y}\frac{d^{2}y}{dt^{2}}+\lambda_{y}\frac{dy}{dt}$
$\displaystyle=F_{{\rm d},y}-F_{\eta,y},$ $\displaystyle
m_{z}\frac{d^{2}z}{dt^{2}}+\lambda_{z}\frac{dz}{dt}$
$\displaystyle=-F_{\eta,z}.$ (S59)

For the “model A” dynamics without the inertia term ($m_{i}=0$), the equation
of motion does not have the acceleration term, leading to the force balance
equation;

$\displaystyle\lambda_{y}\frac{dy}{dt}=F_{{\rm d},y}-F_{\eta,y},\ \
\lambda_{z}\frac{dz}{dt}=-F_{\eta,z}.$ (S60)

Using $F_{\rm d,z}=0$ and $F_{\rm d,y}=F_{\rm d}$, the force balance equation
is solved for the coordinate of the vortex core,

$\displaystyle z(t)$ $\displaystyle=\int dt\ \ \lambda_{z}^{-1}F_{\eta,z}$
$\displaystyle y(t)$ $\displaystyle=\int dt\ \ \lambda_{y}^{-1}(F_{\rm
d}-F_{\eta,y}).$ (S61)

The mean displacement is

$\displaystyle\langle z(t)\rangle$ $\displaystyle=0;\ \ \langle
y(t)\rangle=\int dt\ \ \lambda_{y}^{-1}F_{\text{d}}.$ (S62)

Eq. (S62) leads to a conclusion: the mean velocity $v(R_{\text{d}})$ of a
vortex in a contracting vortex-antivortex pair with distance $R_{\text{d}}$
reads

$\displaystyle v_{y}(R_{\text{d}})=\frac{d\langle
y(t)\rangle}{dt}=-\frac{dR_{\text{d}}}{dt}=\lambda_{y}^{-1}(R_{\text{d}})F_{\text{d}}(R_{\text{d}}).$
(S63)

For a coarsening system, a single typical time-dependent correlation length
$L(t)$, which is, physically, the typical inter-defect distance or the defect
size, characterizes the system. Thus, the typical velocity reads $v(R_{\rm
d}\sim L)$ [8, 11, 10, 10]. Equating it to the rate of change of $L$ gives the
following Eq. (S64), an equation to obtain the growth exponent $p$ of the
growing length $L(t)\propto t^{p}$ [8, 11, 10, 10].

$\displaystyle\frac{dL}{dt}=\lambda_{y}^{-1}(L)F_{\text{d}}(L).$ (S64)

In Sec. III.2.3, a formal of Eq. (S64) from Eq. (S63) is presented. Eq. (S64)
is referred to as Eq. (5) in the main text, where $\lambda$ denotes
$\lambda_{y}$ for simplicity.

To solve Eq. (S64), $\lambda_{i}(L)$ and $F_{\text{d}}(L)$ need to be
evaluated. In $L\ll\xi_{0}/\sqrt{2\beta}$ limit, the defects behave like
circular vortices in the $U(1)$-symmetric case. Therefore,
$E_{\text{d}}(L)\sim E_{\text{c}}|\alpha|\xi_{0}^{2}\ln(L/\xi)$,
$F_{\text{d}}(L)=|dE_{\text{d}}/dL|=E_{\text{c}}|\alpha|\xi_{0}^{2}/L$. For
$\lambda_{y}$, notice that for $L\ll\xi_{0}/\sqrt{2\beta}$,

$\displaystyle\lambda_{i=y,z}$
$\displaystyle\approx\frac{E_{\text{c}}|\alpha|}{\gamma}\int d^{2}\bm{r}\ \
2|\nabla\theta\cdot\hat{i}|^{2}=\frac{E_{\text{c}}|\alpha|}{\gamma}\int
d^{2}\bm{r}\ \
|\nabla\theta|^{2}=\frac{E_{\text{d}}}{\gamma\xi_{0}^{2}}\sim\frac{E_{\text{c}}|\alpha|}{\gamma}\ln(L/\xi).$
(S65)

Thus,

$\displaystyle L(t)\sim\xi_{0}\sqrt{\frac{\gamma t}{\ln(|\alpha|\gamma t)}}.$
(S66)

The scaling is the same as the 2D $U(1)$-symmetric systems after a global
thermal quench [8].

In the $L\gg\xi_{0}/\sqrt{2\beta}$ limit, the typical attraction force reads
$F_{\text{d}}(L)\sim E_{\text{c}}|\alpha|\xi_{0}e^{-L\sqrt{2\beta}/\xi_{0}}$
according to Eq. (S49). For $\lambda_{y}$,

$\displaystyle\lambda_{y}$
$\displaystyle=\frac{E_{\text{c}}|\alpha|}{\gamma}\int d^{2}\bm{r}\ \
2(\partial_{y}\theta)^{2}=\frac{2E_{\text{c}}|\alpha|}{\gamma}\int_{-L/2}^{L/2}dy\int
dz\ \ (\partial_{y}s)^{2}.$ (S67)

In the $|y|\rightarrow\infty$ limit, by Eq. (S46), the integrand in the
integral of $y$ approaches the following asymptotic form,

$\displaystyle\int dz\ (\partial_{y}s)^{2}=\frac{\sqrt{2\beta}}{\xi_{0}}\int
dk^{\prime}|c(k^{\prime})|^{2}(k^{\prime 2}+1)e^{-\sqrt{k^{\prime
2}+1}|y|\sqrt{2\beta}/\xi_{0}}.$ (S68)

As $|y|$ gets larger, the integrand show an exponential decay. Therefore, the
integral of $y$ is convergent for $L\rightarrow+\infty$. In other words, for
larger $L$, the leading-order contribution of $\lambda_{y}(L)$ is a
$L$-independent constant $\lambda_{y}(+\infty)\equiv\lambda_{0}\sim
E_{\text{c}}|\alpha|/\gamma$. Thus, from Eq. (S64, S49),

$\displaystyle\frac{dL}{dt}=\lambda_{y}^{-1}F_{\text{d}}\sim\gamma\xi_{0}\sqrt{2\beta}e^{-L\sqrt{2\beta}/\xi_{0}}$
(S69)

i.e.

$\displaystyle L\propto\frac{\xi_{0}}{\sqrt{2\beta}}\ln(2\beta\gamma t).$
(S70)

$\bullet$ Defect displacement in the $z$-direction and stability of the
interface: The analysis above assumes that the vortices are confined within
the interface. In fact, from Eq. (S60), the vortices acquire a weak Brownian-
like motion along the $z$-direction due to the thermal noise. Write the
equation of motion in the $z$-direction in the following dimensionless form,

$\displaystyle\frac{d(z/\xi_{0})}{d(\gamma t)}=H_{z}.$ (S71)

Here $H_{z}=-\lambda_{z}^{-1}F_{\eta,z}/(\xi_{0}\gamma)$, whose correlator
reads

$\displaystyle\langle
H_{z}(t)H_{z}(t^{\prime})\rangle=\frac{2T}{\gamma\xi_{0}^{2}}\lambda_{z}^{-1}(R_{\rm
d})\delta(\gamma t-\gamma t^{\prime}).$ (S72)

Here $R_{\rm d}$ is the vortex size. By the expression of $\lambda_{i}$ in Eq.
(S65), $\lambda_{i}=E_{\rm c}|\alpha|\gamma^{-1}h_{i}(R_{\rm d})$ where
$h_{i}(R_{\rm d})=\int d^{2}\bm{r}\ |\partial_{i}\theta|^{2}$ is
dimensionless. Therefore,

$\displaystyle\langle H_{z}(t)H_{z}(t^{\prime})\rangle=\frac{2T}{E_{\rm
c}\xi_{0}^{2}|\alpha|}h_{z}^{-1}(R_{\rm d})\delta(\gamma t-\gamma
t^{\prime})=2\zeta h_{z}^{-1}(R_{\rm d})\delta(\gamma t-\gamma t^{\prime}).$
(S73)

Replacing $R_{\rm d}$ in Eq. (S73) by the typical defect size $L$, we obtain
Eq. (9) in the main text. Eq. (S71) and Eq. (S73) imply the defect motion in
the $z$-direction is negligible because $\zeta$ is small in consideration of
this paper. Thus, the interface is stable.

For the $L\ll\xi_{0}/\sqrt{2\beta}$ case with a nonzero $\beta$, one has
$\lambda_{z}=\lambda_{y}$ and $h_{z}(L)=h_{y}(L)\sim\ln(L/\xi)$, while it is
hard to evaluate $h_{z}(L)$ for the other cases. In the
$L\ll\xi_{0}/\sqrt{2\beta}$ case, after replacing $R_{\rm d}$ in Eqs.(S72,S73)
by its typical value $L$ , we could solve the following scaling of $\langle
z^{2}(L)\rangle$

$\displaystyle\langle z^{2}(L(t))\rangle=$ $\displaystyle
2T\int^{t}\lambda_{z}^{-1}(L(t^{\prime}))\
dt^{\prime}=2T\int^{t}\lambda_{y}^{-1}(L(t^{\prime}))\ dt^{\prime}$
$\displaystyle\stackrel{{\scriptstyle\text{By
eq.}~{}(\ref{estimate_sm})}}{{=}}$ $\displaystyle 2T\int^{L}F_{\rm
d}^{-1}(L^{\prime})\ dL^{\prime}=2T\int^{L}\frac{L^{\prime}}{E_{\rm
c}|\alpha|\xi_{0}^{2}}\ dL^{\prime}\approx\zeta L^{2}.$ (S74)

i.e. $\sqrt{\langle z^{2}(L)\rangle}\approx\sqrt{\zeta}L\ll L$. $\sqrt{\langle
z^{2}(L)\rangle}$ may physically represent the typical value of the defects’
displacement in the $z$-direction when the typical inter-defect distance is
$L$. In that way, $\sqrt{\langle z^{2}(L)\rangle}=\sqrt{\zeta}L\ll L$ means
the interface is stable. For other cases, quantitative evaluation of
$\sqrt{\langle z^{2}(L)\rangle}$ is hard to achieve since $\lambda_{z}(L)$ and
$h_{z}(L)$ are hard to evaluate.

####  III.2.2 3D systems

In 3D systems, vortex string loops contract and $\psi(\bm{r},t)$ satisfies Eq.
(S53). For convenience, we separate the coordinates by $\bm{r}=(\bm{R}_{n},s)$
where $\bm{R}_{n}$ is the local coordinate of the intersecting plane normal to
the string and $s$ is the coordinate along the string. The motion of the
string is characterized by the ansatz
$\psi(\bm{r}=(\bm{R}_{n},s),t)=f(\bm{R}_{n}-\int^{t}\bm{v}(s,t^{\prime})dt^{\prime})$.
Here $f(\bm{R}_{n})$ is the order-parameter profile of a point defect in the
2D plane spanned by $\bm{R}_{n}$. As illustrated in Fig. S7, the ansatz
implies that $\bm{v}(s,t)$ is confined in the 2D plane spanned by
$\bm{R}_{n}$.

Figure S7: A circular vortex string loop contracts under the driving of the
force $\bm{F}_{\rm d}$, which points to the center of the circle. The local
motion of the string is confined in the local $R_{n,1}-R_{n,2}$ plane normal
to the string, which is marked by blue.

To further simplify the problem, we idealize the shape of the vortex-string
loop to be circular with a diameter $R_{\text{d}}$. The energy of such a
vortex-string loop is $E_{\text{s}}=R_{\rm d}E_{\text{d}}(R_{\text{d}})$,
where a constant $\pi$ factor is omitted. Here $E_{\text{d}}(R_{\rm d})$ is
the energy per unit length of the string. $E_{\text{d}}$ equals the energy of
point vortices in the 2D intersecting plane of the string, whose expression is
given in Sec. III.1. Since the vortex string loop is contracting, $R_{\rm d}$
as well as $E_{\text{s}}=R_{\rm d}E_{\text{d}}(R_{\text{d}})$ decrease with
time. The energy dissipation rate of the vortex string loop reads

$\displaystyle\frac{dE_{\text{s}}}{dt}=\frac{dE_{\text{s}}}{dR_{\text{d}}}\frac{dR_{\text{d}}}{dt}=\int_{0}^{R_{\rm
d}}ds\ \ -\bm{F_{\rm d}}(s)\cdot\bm{v}(s).$ (S75)

Here $\bm{F_{\rm d}}(s)$ is a force acting on per unit length of the string,
which comes from the string tension stem from the local curvature. $\bm{F_{\rm
d}}(s)$ points to the center of the circle and drives the loop to contract, as
shown in Fig. S7. Since the loop is circular, the magnitude of $\bm{F_{\rm
d}}$ is independent on $s$ with the following expression,

$\displaystyle F_{\text{d}}=\frac{1}{R_{\text{d}}}\frac{dE_{\rm
s}}{dR_{\text{d}}}=\frac{E_{\text{d}}(R_{\text{d}})}{R_{\text{d}}}+\frac{dE_{\text{d}}(R_{\text{d}})}{dR_{\text{d}}}.$
(S76)

In our theory, the strings are nearly confined in the interface, i.e. the
displacement of the strings from the $z=0$ interface is negligible compared to
$R_{\rm d}$. Therefore, the direction of $\bm{F_{\rm d}}(s)$ is nearly
confined in the $x-y$ plane.

Consider an infinitely small segment of the string with length $\epsilon$
located at $s=s_{0}$. Set its energy as $E_{\epsilon}$. Eq. (S75) impies that
the energy dissipation rate of the small segment reads

$\displaystyle\frac{dE_{\epsilon}}{dt}=\int_{s_{0}}^{s_{0}+\epsilon}ds\ \
-\bm{F_{\rm d}}(s)\cdot\bm{v}(s)=-\epsilon\bm{F}_{\rm
d}(s_{0})\cdot\bm{v}(s_{0}).$ (S77)

Moreover, by Eq. (S53) and the ansatz of
$\psi(\bm{r}=(\bm{R}_{n},s_{0}),t)=f(\bm{R}_{n}-\int^{t}\bm{v}(s_{0},t^{\prime})dt^{\prime})$,
the energy dissipation rate of the small segment can be given by the inertia
term, and friction and Langevin forces acting on the small segment.

$\displaystyle\frac{dE_{\epsilon}}{dt}$
$\displaystyle=\int_{s_{0}}^{s_{0}+\epsilon}ds\int d^{2}\bm{R}_{n}\ \
\frac{\delta F}{\delta\psi}\frac{\partial\psi}{\partial t}+\frac{\delta
F}{\delta\psi^{*}}\frac{\partial\psi^{*}}{\partial t}$
$\displaystyle=\int_{s_{0}}^{s_{0}+\epsilon}ds\int d^{2}\bm{R}_{n}\ \
\Big{(}-\frac{E_{\rm c}}{\gamma^{2}_{0}}\frac{\partial^{2}\psi^{*}}{\partial
t^{2}}-\frac{E_{\text{c}}}{\gamma}\frac{\partial\psi^{*}}{\partial
t}+E_{\text{c}}\eta^{*}\Big{)}\\!\ \frac{\partial\psi}{\partial
t}+\Big{(}-\frac{E_{\rm c}}{\gamma^{2}_{0}}\frac{\partial^{2}\psi}{\partial
t^{2}}-\frac{E_{\text{c}}}{\gamma}\frac{\partial\psi}{\partial
t}+E_{\text{c}}\eta\Big{)}\\!\ \frac{\partial\psi^{*}}{\partial t}$
$\displaystyle=-\sum_{i,j,k=R_{n,1},R_{n,2}}\int_{s_{0}}^{s_{0}+\epsilon}ds\int
d^{2}\bm{R}_{n}\ \
\Big{(}-\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}\partial_{k}f^{*}){\bm{v}}_{i}{\bm{v}}_{k}+\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}f^{*})\frac{\partial{\bm{v}}_{i}}{\partial
t}+\frac{E_{\text{c}}}{\gamma}(\partial_{i}f^{*})\bm{v}_{i}+E_{\text{c}}\eta^{*}\Big{)}\\!\
(\partial_{j}f)\bm{v}_{j}$ $\displaystyle\hskip
142.26378pt+\Big{(}-\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}\partial_{k}f){\bm{v}}_{i}{\bm{v}}_{k}+\frac{E_{c}}{\gamma^{2}_{0}}(\partial_{i}f)\frac{\partial{\bm{v}}_{i}}{\partial
t}+\frac{E_{\text{c}}}{\gamma}(\partial_{i}f)\bm{v}_{i}+E_{\text{c}}\eta\Big{)}\\!\
(\partial_{j}f^{*})\bm{v}_{j}$
$\displaystyle=(\epsilon\bm{F}_{\lambda}-\bm{F}_{\eta}-\epsilon{\bm{A}})\cdot\bm{v}(s_{0}).$
(S78)

Here $\epsilon{\bm{A}}=\sum_{i=y,z}\epsilon m_{i}a_{i}\hat{i}$ with
${\bm{a}}=d{\bm{v}}(s_{0})/dt$, is the inertia term for the small segment.
$\epsilon\bm{F}_{\lambda}=-\sum_{i=y,z}\epsilon\lambda_{i}v_{i}(s_{0})\hat{i}$
is the frictional force. $\bm{F}_{\eta}$ is the Langevin force. $m_{i}$ and
$\lambda_{i}$ are the inertia mass and the friction coefficient in the
$i$-direction, respectively. The inertia mass, friction coefficients, and the
Langevin force are given as follows;

$\displaystyle m_{i}$ $\displaystyle=\frac{E_{\text{c}}}{\gamma^{2}_{0}}\int
d^{2}\bm{R_{n}}\ \ 2|\partial_{i}f|^{2},$ (S79) $\displaystyle\lambda_{i}$
$\displaystyle=\frac{E_{\text{c}}}{\gamma}\int d^{2}\bm{R_{n}}\ \
2|\partial_{i}f|^{2},$ (S80) $\displaystyle\bm{F}_{\eta}$
$\displaystyle=E_{\text{c}}\int_{s_{0}}^{s_{0}+\epsilon}ds\int
d^{2}\bm{R}_{n}\ \ (\eta^{*}\nabla_{n}f+\eta\nabla_{n}f^{*})..$ (S81)

Here $\nabla_{n}=\partial/\partial\bm{R}_{n}$. Without loss of generality, we
set the tangent of the string at $s=s_{0}$ is along the $x$-direction and
$\bm{R}_{n}=(y,z)$. Therefore, $\bm{F}_{\rm d}=F_{\rm d}\hat{y}$. Equating Eq.
(S77) and Eq. (S78), we obtain an equation of motion for the small segment of
the string under the driving force $\bm{F}_{\rm d}$.

$\displaystyle m_{y}\frac{d^{2}y}{dt^{2}}+\lambda_{y}\frac{dy}{dt}$
$\displaystyle=F_{{\rm d},y}-\epsilon^{-1}F_{\eta,y},$ $\displaystyle
m_{z}\frac{d^{2}z}{dt^{2}}+\lambda_{z}\frac{dz}{dt}$
$\displaystyle=-\epsilon^{-1}F_{\eta,z}.$ (S82)

For the “model A” dynamics without the inertia term ($m_{i}=0$), we obtain the
force-balance equation

$\displaystyle\lambda_{y}\frac{dy}{dt}=F_{{\rm d},y}-\epsilon^{-1}F_{\eta,y},\
\ \lambda_{z}\frac{dz}{dt}=-\epsilon^{-1}F_{\eta,z}.$ (S83)

The force balance equation is solved for the coordinate of the vortex core,

$\displaystyle z(t)$ $\displaystyle=\int dt\ \
-\epsilon^{-1}\lambda_{z}^{-1}F_{\eta,z}$ $\displaystyle y(t)$
$\displaystyle=\int dt\ \
\lambda_{y}^{-1}(-\epsilon^{-1}F_{\eta,y}+F_{\text{d}}).$ (S84)

The mean displacement is

$\displaystyle\langle z(t)\rangle$ $\displaystyle=0;\ \ \langle
y(t)\rangle=\int dt\ \ \lambda_{y}^{-1}F_{\text{d}}.$ (S85)

Eq. (S85) leads to a conclusion: that the mean velocity of the string motion
reads $v(R_{\text{d}})=-dR_{\rm
d}/dt=\lambda_{y}^{-1}(R_{\text{d}})F_{\text{d}}(R_{\text{d}})$. Similar to
the analysis for the 2D systems, taking $R_{\text{d}}$ by its typical value
$L$ and equating $dL/dt=v(L)$ gives Eq. (S64) [8]. See Sec. III.2.3 for a more
formal derivation.

To solve Eq. (S64), $\lambda_{i}(L)$ and $F_{\text{d}}(L)$ need to be
evaluated. In the presence of the $U(1)$-symmetry ($\beta=0$) or in the
$L\ll\xi_{0}/\sqrt{2\beta}$ limit when $\beta\neq 0$, $E_{\text{d}}(L)\propto
E_{\text{c}}|\alpha|\xi_{0}^{2}\ln(L/\xi)$ and the leading order contribution
of $F_{\text{d}}(L)$ is the $E_{\text{d}}(L)/L$ term in Eq. (S76). Therefore,

$\displaystyle F_{\text{d}}$ $\displaystyle\sim
E_{\text{c}}|\alpha|\xi_{0}^{2}\frac{1}{L}\ln(\frac{L}{\xi}).$ (S86)

For $\lambda_{i=y,z}$, as shown in the analysis of 2D systems,

$\displaystyle\lambda_{z/y}$
$\displaystyle\sim\frac{E_{\text{c}}|\alpha|}{\gamma}\ln(\frac{L}{\xi}).$
(S87)

Thus, from Eq. (S64), we obtain

$\displaystyle L(t)\sim\xi_{0}\sqrt{\gamma t}.$ (S88)

In the $L\gg\xi_{0}/\sqrt{2\beta}$ limit, $E_{\rm d}(L)\sim E_{\rm
c}|\alpha|\xi_{0}^{2}$ according to the calculation in Sec. III.1. The leading
order contribution of $F_{\text{d}}(L)$ is the $E_{\text{d}}(L)/L$ term in Eq.
(S76). Therefore

$\displaystyle F_{\text{d}}$ $\displaystyle\sim
E_{\text{c}}|\alpha|\xi_{0}^{2}\frac{1}{L}.$ (S89)

For $\lambda_{y}$, according to the calculation in Eq. (S67),

$\displaystyle\lambda_{y}$ $\displaystyle=\lambda_{0}\propto
E_{\text{c}}|\alpha|/\gamma$ (S90)

Therefore, from Eq. (S64), we obtain

$\displaystyle L(t)\sim\xi_{0}\sqrt{\gamma t}.$ (S91)

The analysis reveals a discrepancy between 2D and 3D systems. In 2D systems
with $\beta\neq 0$, as $L$ gets larger, the scaling of $L(t)$ undergoes a
crossover from a diffusive behavior to a subdiffusive behavior. See Eq. (S66,
S70). However, in 3D systems, $L(t)$ always shows a diffusive behavior.

The discrepancy could be explained in the following way. In 2D systems, as
shown above, the $L$-dependence of the attraction force $F_{\text{d}}$ shows a
$1/L$ (long range) to $e^{-L\sqrt{2\beta}/\xi_{0}}$ (short range) crossover as
$L$ goes across the scale $\xi_{0}/\sqrt{2\beta}$. The crossover of
$F_{\text{d}}$ makes the motion and annihilation of vortices slower with time
going by, leading to the diffusive-to-subdiffusive crossover of the coarsening
dynamics. In 3D systems, the dominant contribution of the driving force
$F_{\text{d}}(L)$ always comes from $E_{\rm d}(L)/L$, the first term in Eq.
(S76), which is from the derivative of the string-loop diameter $L$ rather
than $E_{\rm d}(L)$. (Physically, it means that the driving force
$F_{\text{d}}(L)$ is dominated by the tension of the string.) Moreover, the
$L$-dependence of $E_{\text{d}}(L)$ is the same as $\lambda_{y}(L)$ in both
$L\ll\xi_{0}/\sqrt{2\beta}$ and $L\gg\xi_{0}/\sqrt{2\beta}$ limits. Therefore,
we always obtain $dL/dt\sim\gamma\xi_{0}^{2}/L$, leading to the result
$L(t)\sim\xi_{0}\sqrt{\gamma t}$ .

$\bullet$ Defect displacement in the $z$-direction and stability of the
interface: Similar to the analysis for 2D systems, the discussion above
assumes the strings are confined within the 2D interface. In fact, the string-
loop could acquire a weak Brownian-like motion in the $z$-direction due to the
thermal noise. Consider the center-of-mass motion in the $z$-direction of a
macroscopical segment of the string loop with diameter $R_{\rm d}$. The length
of the segment is $\sim R_{\rm d}$. The equation of motion is derived by
simply replacing $\epsilon$ by $R_{\rm d}$ in Eq. (S78, S81, S83), where
$R_{\rm d}$ is the diameter of the loop. Written in the dimensionless form,
the equation of motion reads

$\displaystyle\frac{d(z/\xi_{0})}{d(\gamma t)}=H_{z}$ (S92)

Here

$\displaystyle H_{z}=-(\xi_{0}\gamma)^{-1}(R_{\rm
d}\lambda_{z})^{-1}E_{\text{c}}\int_{0}^{R_{\rm d}}ds\int d^{2}\bm{R}_{n}\ \
(\eta^{*}\partial_{z}f+\eta\partial_{z}f^{*}).$ (S93)

$H_{z}$ satisfies

$\displaystyle\langle H_{z}(t)H_{z}(t^{\prime})\rangle=\frac{2T}{R_{\rm
d}\xi_{0}^{2}\gamma}\lambda_{z}^{-1}(R_{\rm d})\delta(\gamma t-\gamma
t^{\prime}).$ (S94)

By the expression of $\lambda_{i}$ in Eq. (S81), $\lambda_{i}=E_{\rm
c}|\alpha|\gamma^{-1}h_{i}(R_{\rm d})$ where $h_{i}(R_{\rm d})=\int
d^{2}\bm{R}_{n}\ |\partial_{i}\theta|^{2}$ is dimensionless. Therefore,

$\displaystyle\langle H_{z}(t)H_{z}(t^{\prime})\rangle$
$\displaystyle=\frac{2T}{R_{\rm d}\xi_{0}^{2}E_{\rm
c}|\alpha|}h_{z}^{-1}\delta(\gamma t-\gamma
t^{\prime})=\frac{2T}{\xi_{0}^{3}E_{\rm c}\sqrt{|\alpha|}}\Big{(}\frac{R_{\rm
d}}{\xi}h_{z}(R_{\rm d})\Big{)}^{-1}\delta(\gamma t-\gamma t^{\prime})$
$\displaystyle=2\zeta\Big{(}\frac{R_{\rm d}}{\xi}h_{z}(R_{\rm
d})\Big{)}^{-1}\delta(\gamma t-\gamma t^{\prime}).$ (S95)

Replacing $R_{\rm d}$ in Eq. (S95) by the typical defect size $L$, we obtain
Eq. (9) in the main text. Eq. (S92) and Eq. (S95) imply the defect motion in
the $z$-direction is negligible because $\zeta$ is small in consideration of
this paper. Thus, the interface is stable.

For the $U(1)$ symmetric case as well as the $L\ll\xi_{0}/\sqrt{2\beta}$ case
with a nonzero $\beta$, one has $\lambda_{z}=\lambda_{y}$ and
$h_{z}(L)=h_{y}(L)\sim\ln(L/\xi)$, while it is hard to evaluate $h_{z}(L)$ for
the other cases. In the $L\ll\xi_{0}/\sqrt{2\beta}$ case, after replacing
$R_{\rm d}$ in Eqs.(S94,S95) by its typical value $L$ , we could solve the
following scaling of $\langle z^{2}(L)\rangle$

$\displaystyle\langle z^{2}(L(t))\rangle=$ $\displaystyle
2T\int^{t}L^{-1}(t^{\prime})\lambda_{z}^{-1}(L(t^{\prime}))\
dt^{\prime}=2T\int^{t}L^{-1}(t^{\prime})\lambda_{y}^{-1}(L(t^{\prime}))\
dt^{\prime}$ $\displaystyle\stackrel{{\scriptstyle\text{By
eq.}~{}(\ref{estimate_sm})}}{{=}}$ $\displaystyle 2T\int^{L}L^{\prime-1}F_{\rm
d}^{-1}(L^{\prime})\ dL^{\prime}=2T\int^{L}_{2\xi}\frac{1}{E_{\rm
c}|\alpha|\xi_{0}^{2}\ln(L^{\prime}/\xi)}\ dL^{\prime}$ $\displaystyle=$
$\displaystyle
2\zeta\xi^{2}\int_{2}^{L/\xi}\frac{1}{\ln(L^{\prime}/\xi)}d(L^{\prime}/\xi)\stackrel{{\scriptstyle
L\gg\xi}}{{=}}2\zeta\frac{\xi L}{\ln(L/\xi)}$ (S96)

Here we set a lower limit of $2\xi$ in the integration over $L^{\prime}$ to
prevent the unphysical singularity. Eq. (S96) gives $\sqrt{\langle
z^{2}(L)\rangle}\ll L$. $\sqrt{\langle z^{2}(L)\rangle}$ may physically
represent the typical value of the defects’ displacement in the $z$-direction
when the typical inter-defect distance is $L$. In that way, $\sqrt{\langle
z^{2}(L)\rangle}\ll L$ means the interface is stable. For other cases,
quantitative evaluation of $\sqrt{\langle z^{2}(L)\rangle}$ is hard to achieve
since $\lambda_{z}(L)$ and $h_{z}(L)$ are hard to evaluate.

####  III.2.3 Derivation of Eq. (S64) from Eq. (S63).

This section discusses how Eq. (S64) is derived From Eq. (S63). The interface
contains random domains of size $R_{\text{d}}$. $R_{\text{d}}$ satisfies a
probability distribution function $p(R_{\text{d}})$ with
$\int_{0}^{+\infty}p(R_{\text{d}})dR_{\text{d}}=1$. Suppose the coarsening
dynamics is self-similar with a single length scale $L(t)$, we obtain
$p(R_{\text{d}})=g(R_{\text{d}}/L)/L$ where $g$ is a function satisfying
$\int_{0}^{+\infty}g(x)dx=1$.

Eq. (S63) implies that the typical velocity of the defect motion depends on
the defect size (domain size). Statistically, when the typical correlation
length is $L(t)$, the typical velocity of the defect’s motion reads
$v(L)=F_{\text{d}}(L)/\lambda_{y}(L)$. Therefore, within a time $dt$, the
number of the annihilated domains reads

$\displaystyle
dN_{\text{d}}=-N_{\text{d}}\int_{0}^{v(L)dt}p(R_{\text{d}})dR_{\text{d}}=-N_{\text{d}}\int_{0}^{v(L)dt/L}g(x)dx.$
(S97)

Here $N_{\text{d}}$ denotes the number of domains. Suppose the system size is
$l$. In the 1D interface of 2D systems, $N_{\text{d}}\approx l/L$. In the 2D
interface of 3D systems, $N_{\text{d}}\approx l^{2}/L^{2}$. Obviously, in both
cases, $dN_{\text{d}}/N_{\text{d}}\approx-dL/L$. Thus, Eq. (S97) implies

$\displaystyle\frac{dL}{L}\approx g(0)\frac{v(L)}{L}dt,$ (S98)

for $v(L)dt/L\ll 1$. Neglecting the dimensionless constant $g(0)$, we derive
Eq. (S64).

##  IV Effects of the $\partial_{t}^{2}\psi$ term

###  IV.1 Suppression of the $\partial_{t}^{2}\psi$ term when the temperature
is near $T_{\rm c}$.

The “model A” dynamics used in Eq. (S2) neglects the inertia term, the second-
order time-derivative $\partial_{t}^{2}\psi$. In fact, as mentioned in Eq.
(S1), the inertia term is crucial to form a domain wall. In this section, we
discuss the effect of the inertia term and analyze how it modifies our
results. With the inertia term, the dynamical equation reads

$\displaystyle\frac{1}{\gamma_{0}^{2}}\partial^{2}_{t}\psi+\frac{1}{\gamma}\partial_{t}\psi=-\frac{1}{E_{\text{c}}}\frac{\delta
F}{\delta\psi^{*}}+\eta=\xi_{0}^{2}\nabla^{2}\psi-2\alpha\psi-2\beta
i\psi_{2}-2|\psi|^{2}\psi+\eta.$ (S99)

Here $\gamma_{0}$ is a parameter controlling the strength of the inertia.
Typically, one has $\gamma,\,\gamma_{0}\sim T$,
$\alpha\sim(T-T_{\text{c}})/T_{\text{c}}$ for pure electronic orders arising
from Fermi surface nesting, see Ref. [16] and Appendix E of Ref. [17]. The
following scaling analysis shows that when the temperature $T$ is close to
$T_{\text{c}}$, the inertia term is effectively suppressed. Eq. (S99) is
invariant under the following rescaling;

$\displaystyle\alpha$ $\displaystyle=\bar{\alpha}b^{2},\ \
\beta=\bar{\beta}b^{2},\ \ E_{\text{c}}=\bar{E_{\text{c}}}b^{-4},\ \
\psi=\bar{\psi}b,$ $\displaystyle\ \ \eta$ $\displaystyle=\bar{\eta}b^{3},\ \
\xi_{0}=\bar{\xi_{0}}b,\ \ \gamma=\bar{\gamma}b^{-2},\ \
\gamma_{0}=\bar{\gamma_{0}}b^{-1}.$ (S100)

The rescaled equation reads,

$\displaystyle\frac{1}{\bar{\gamma_{0}}^{2}}\partial^{2}_{t}\bar{\psi}+\frac{1}{\bar{\gamma}}\partial_{t}\bar{\psi}=\bar{\xi_{0}}^{2}\nabla^{2}\bar{\psi}-2\bar{\alpha}\bar{\psi}-2\bar{\beta}i\bar{\psi_{2}}-2|\bar{\psi}|^{2}\bar{\psi}+\bar{\eta}.$
(S101)

with the noise correlator
$\langle\bar{\eta}_{i}(\bm{r},t)\bar{\eta}_{i^{\prime}}(\bm{r^{\prime}},t^{\prime})\rangle=2T(\bar{\gamma}\bar{E_{\text{c}}})^{-1}\delta_{ii^{\prime}}\delta(\bm{r}-\bm{r^{\prime}})\delta(t-t^{\prime})$.
We could set $\bar{\alpha}=-1$. If $T$ is close to $T_{\text{c}}$, $|\alpha|$
and $b=\sqrt{|\alpha|}$ will be small. Thus,
$\bar{\gamma}^{2}/\bar{\gamma_{0}}^{2}=b^{2}\gamma^{2}/\gamma_{0}^{2}$ is
suppressed to be small. Since we discuss the dynamics happening after $t_{0}$,
as long as $T(t>t_{0})=T_{0}$ is near $T_{\text{c}}$, the inertia term is much
weaker than the damping term and our analysis using the “model A” dynamics is
valid. The condition
$\bar{\gamma}^{2}/\bar{\gamma_{0}}^{2}=|\alpha|\gamma^{2}/\gamma_{0}^{2}=1$
gives a lower critical temperature $T_{\gamma}$, below which the inertia term
is non-negligible. Moreover, another condition $\zeta=G|\alpha|^{D/2-2}\geq 1$
sets an upper critical temperature $T_{\rm G}$, with $T<T_{\rm G}$ being the
Ginzburg criterion in the conventional theory of critical phenomenon [6]. If
$T_{\rm G}<T_{0}<T_{\rm c}$, the thermal fluctuation is strong and the mean-
field solution of $\psi=-\sqrt{\alpha}$ for $t<t_{\rm pump}$ is invalid and
the domain wall will not form. The successful formation of the domain wall in
the CDW materials proves those systems are inside the Ginzburg criterion. In
summary, the valid criterion of the analysis in this paper is
$T_{\gamma}<T_{0}<T_{\rm G}$, shown by Fig. S8.

Figure S8: Valid criteria of the analysis in this paper is marked by green.

Note that the inertia term is suppressed for $t>t_{0}$ by $T_{0}$. While
within $t_{\rm pump}<t<t_{\rm off}$, the pump lifts the temperature to $T_{\rm
H}$ and $\alpha\sim(T_{\rm H}-T_{\rm c})/T_{\rm c}$ may be larger, making the
inertia term significant enough to form the domain wall.

###  IV.2 Quantitative modifications of the first and second stage dynamics by
the $\partial_{t}^{2}\psi$ term.

In the following discussion, we discuss how the weak inertia modifies our
results quantitatively for $t>t_{0}$.

$\bullet$ First-stage dynamics. For the first-stage dynamics, we calculate the
correlation $C(\bm{R},t)=C_{1}(\bm{R},t)+C_{2}(\bm{R},t)$ defines in Eq. (S12,
S13, S14). A calculation is shown in Sec. IV.3. $C_{2}(\bm{R},t)$ is evaluated
in the $\gamma t\gg 1/|\alpha|$ limit,

$\displaystyle C_{2}(\bm{R},t)$
$\displaystyle=\frac{\zeta|\alpha|}{2\kappa^{2}c/|\alpha|}\frac{e^{2c\gamma
t}}{\left(\sqrt{8\pi|\alpha|\gamma
t/\kappa}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma t/\kappa}}.$ (S102)

Besides, $C_{1}(\bm{R},t)$ is given by the following integral.

$\displaystyle C_{1}(\bm{R},t)$ $\displaystyle=e^{2c\gamma
t}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K}(0)e^{-2\xi_{0}^{2}\gamma
tK^{2}/\kappa}e^{i\bm{K}\cdot\bm{R}}$ (S103)

Here
$\kappa=\sqrt{1+4(|\alpha|-2\beta)\gamma^{2}/\gamma_{0}^{2}},c=2^{-1}(\kappa-1)\gamma_{0}^{2}/\gamma^{2}$.

Compared with Eq. (S15, S13), Eq. (S102 ,S103) only contains corrections of
the coefficients. Espectially, the time-dependent correlation length
$\xi_{0}\sqrt{\gamma t}$ shown in the results of Eq. (S15, S34, S35, S36) is
renormalized to $\xi_{0}\sqrt{\gamma t}/\sqrt{\kappa}$.

$\bullet$ Second-stage dynamics. For the second-stage dynamics, the equations
of motion with the inertia term are given in Eqs. (S59,S82) for 2D and 3D
systems. have nonzero $m_{i}$. An average of Eqs. (S59,S82) over the thermal
noise $\eta$ gives the equation of motion

$\displaystyle m_{y}\frac{d^{2}\langle
y\rangle}{dt^{2}}+\lambda_{y}\frac{d\langle y\rangle}{dt}$
$\displaystyle=F_{{\rm d},y}$ $\displaystyle m_{z}\frac{d^{2}\langle
z\rangle}{dt^{2}}+\lambda_{z}\frac{d\langle z\rangle}{dt}$ $\displaystyle=0.$
(S104)

Therefore, Eq. (S63) is modified to

$\displaystyle m_{y}\frac{dv_{y}(R_{\rm d})}{dt}+\lambda_{y}v_{y}(R_{\rm d})$
$\displaystyle=F_{\text{d}}(R_{\rm d}).$ (S105)

Taking $v_{y}(R_{\rm d}\sim L)=dL/dt$, Eq. (S64) is modified to

$\displaystyle m_{y}\frac{d^{2}L}{dt^{2}}+\lambda_{y}\frac{dL}{dt}$
$\displaystyle=F_{\text{d}}(L).$ (S106)

Eq. (S106) is the Eq. (10) in the main text, where $\lambda$, $m$ in the main
text denote $\lambda_{y},m_{y}$ here for simplicity. Note that all the results
of $L(t)$ satisfying Eq. (S64) obtained in this work satisfy Eq. (S106) in the
long time limit, at which time the term $m_{y}(L)d^{2}L/dt^{2}$ is much
smaller than the other two terms.

###  IV.3 Derivation of Eq. (S102, S103).

In the derivation in this section, we use the rescaled notation used in the
dynamical equation Eq. (S133) with $\bar{\gamma}^{2}/\bar{\gamma}_{0}^{2}\ll
1$. Moreover, we drop the bar symbol on the top of the characters for
convenience.

In the first-stage dynamics, the dynamical equation Eq. (S4) with an inertia
term reads

$\displaystyle\frac{1}{\gamma_{0}^{2}}\partial^{2}_{t}\phi_{1/2}+\frac{1}{\gamma}\partial_{t}\phi_{1/2}$
$\displaystyle=-(\hat{L}_{1/2}-\xi_{0}^{2}\nabla_{\bm{R}}^{2})\phi_{1/2}+\eta_{1/2}.$
(S107)

Here $\gamma^{2}/\gamma_{0}^{2}\ll 1$. Only the modes
$\phi_{2}^{\bm{K}}(\bm{r})=\text{sech}(z/\xi)e^{i\bm{K}\cdot\bm{R}}$ with
$K<\sqrt{|\alpha|-2\beta}/\xi_{0}$ will be exponentially amplified with time.
Writing $\phi_{2}(\bm{r})=\varphi(\bm{R},t)\text{sech}(z/\xi)$, Eq. (S7, S9)
are modified to be

$\displaystyle\frac{1}{\gamma_{0}^{2}}\partial_{t}^{2}\varphi(\bm{R},t)+\frac{1}{\gamma}\partial_{t}\varphi(\bm{R},t)$
$\displaystyle=(|\alpha|-2\beta+\xi_{0}^{2}\nabla_{\bm{R}}^{2})\varphi(\bm{R},t)+\eta_{\varphi}(\bm{R},t),$
$\displaystyle\frac{1}{\gamma_{0}^{2}}\partial_{t}^{2}\tilde{\varphi}(\bm{K},t)+\frac{1}{\gamma}\partial_{t}\tilde{\varphi}(\bm{K},t)$
$\displaystyle=(|\alpha|-2\beta-\xi_{0}^{2}K^{2})\tilde{\varphi}(\bm{K},t)+\tilde{\eta}_{\varphi}(\bm{K},t).$
(S108)

Construct the following two functions,

$\displaystyle A_{1}(\bm{K},t)$
$\displaystyle=\frac{\gamma^{2}}{\gamma_{0}^{2}}\frac{1}{\gamma}\partial_{t}\tilde{\varphi}(\bm{K},t)+[\frac{\gamma^{2}}{\gamma_{0}^{2}}a_{1}(\bm{K})+1]\tilde{\varphi}(\bm{K},t),$
$\displaystyle A_{2}(\bm{K},t)$
$\displaystyle=\frac{\gamma^{2}}{\gamma_{0}^{2}}\frac{1}{\gamma}\partial_{t}\tilde{\varphi}(\bm{K},t)+[\frac{\gamma^{2}}{\gamma_{0}^{2}}a_{2}(\bm{K})+1]\tilde{\varphi}(\bm{K},t),$
(S109)

with

$\displaystyle
a_{1}(\bm{K})=\frac{-1+\Delta(K)}{2}\frac{\gamma_{0}^{2}}{\gamma^{2}},$
$\displaystyle
a_{2}(\bm{K})=\frac{-1-\Delta(K)}{2}\frac{\gamma_{0}^{2}}{\gamma^{2}}.$ (S110)

Here
$\Delta(K)=\sqrt{1+4(|\alpha|-2\beta-\xi_{0}^{2}K^{2})\gamma^{2}/\gamma_{0}^{2}}$.
From Eq. (S108), $A_{1/2}$ satisfies

$\displaystyle\frac{1}{\gamma}\partial_{t}A_{1/2}=a_{1/2}(\bm{K})A_{1/2}+\tilde{\eta}_{\varphi}(\bm{K},t).$
(S111)

The solution reads

$\displaystyle
A_{1/2}(\bm{K},t)=A_{1/2}(\bm{K},t_{0})e^{\gamma\int_{t_{0}}^{t}a_{1/2}(K)dt^{\prime}}+\gamma\int_{t_{0}}^{t}\tilde{\eta}_{\varphi}(\bm{K},t^{\prime})e^{\gamma\int_{t^{\prime}}^{t}a_{1/2}(K)dt^{\prime\prime}}dt^{\prime}.$
(S112)

Therefore, the function $A_{1}(\bm{K},t)$ will be exponentially amplified with
time if $K<\sqrt{|\alpha|-2\beta}/\xi_{0}$ since $a_{1}(\bm{K})>0$, while
$A_{2}(\bm{K},t)$ will exponentially decay with time since $a_{2}(\bm{K})<0$.
Notice that the decay of $A_{2}(\bm{K},t)$ is extremely rapid, since
$a_{2}(\bm{K})$ is large with $a_{2}(\bm{K})\approx\gamma_{0}^{2}/\gamma^{2}$
in the $\gamma^{2}/\gamma_{0}^{2}\ll 1$ limit. Therefore, by the relation

$\displaystyle
A_{1}(\bm{K},t)-A_{2}(\bm{K},t)=\Delta(K)\tilde{\varphi}(\bm{K},t),$ (S113)

$A_{1}(\bm{K},t)$ could be approximated by
$A_{1}(\bm{K},t)\approx\Delta(K)\tilde{\varphi}(\bm{K},t).$ for $t>t_{0}$.
Thereby, from Eq. (S111), we obtain

$\displaystyle\frac{1}{\gamma}\partial_{t}\tilde{\varphi}(\bm{K},t)=a_{1/2}(\bm{K})\tilde{\varphi}(\bm{K},t)+\frac{\tilde{\eta}_{\varphi}(\bm{K},t)}{\Delta(K)},$
(S114)
$\displaystyle\tilde{\varphi}(\bm{K},t)=\tilde{\varphi}(\bm{K},t_{1})e^{U_{K}(t,t_{0})}+\gamma\int_{t_{0}}^{t}\frac{\tilde{\eta}_{\varphi}(\bm{K},t^{\prime})}{\Delta(K)}e^{U_{K}(t,t^{\prime})}dt^{\prime}$
(S115)

Here $U_{K}(t_{1},t_{2})=\gamma\int_{t_{2}}^{t_{1}}a_{1}(K)dt$. Define the
amplitude
$D_{K}(t)=\langle\tilde{\varphi}(-\bm{K},t)\tilde{\varphi}(\bm{K},t)\rangle$.
We obtain

$\displaystyle D_{K}(t)$
$\displaystyle=\frac{1}{l^{D-1}}e^{2U_{K}(t,t_{0})}\Big{(}l^{D-1}D_{K}(t_{0})+\frac{T\gamma}{E_{\text{c}}\xi}\frac{1}{\Delta^{2}(K)}\int_{t_{0}}^{t}dt^{\prime}e^{-2U_{K}(t^{\prime},t_{0})}\Big{)}$
(S116)
$\displaystyle=\frac{1}{l^{D-1}}e^{2U_{K}(t,t_{0})}\Big{(}l^{D-1}D_{K}(t_{0})+\frac{T}{2E_{\text{c}}\xi}\frac{1}{\Delta^{2}(K)}\frac{1}{a_{1}(K)}(1-e^{-2U_{K}(t,t_{0})})\Big{)}.$

Here $l$ is the system size. In the following derivation, we set $t_{0}=0$ for
convenience. Apply the inverse Fourier transform to $D_{K}(t)$ to obtain the
function $C(\bm{R},t)$.

$\displaystyle C(\bm{R},t)$
$\displaystyle=\sum_{\bm{K}}D_{K}(t)e^{i\bm{K}\cdot\bm{R}}\equiv
l^{D-1}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}D_{K}(t)e^{i\bm{K}\cdot\bm{R}}$
$\displaystyle=C_{1}(\bm{R},t)+C_{2}(\bm{R},t),$ (S117)

with

$\displaystyle C_{1}(\bm{R},t)$
$\displaystyle=\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K}(0)e^{2U_{K}(t,0)}e^{i\bm{K}\cdot\bm{R}}$
(S118) $\displaystyle C_{2}(\bm{R},t)$
$\displaystyle=\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}\frac{T\sqrt{|\alpha|}}{2E_{\text{c}}\xi_{0}}\frac{1}{\Delta^{2}(K)}\frac{1-e^{-2U_{K}(t,0)}}{a_{1}(K)}e^{2U_{K}(t,0)}e^{i\bm{K}\cdot\bm{R}}.$
(S119)

To make the calculation of $C_{1}(\bm{R},t)$ and $C_{2}(\bm{R},t)$ practical,
we approximate $a_{1}(K)$ by the following Taylor expansion at $K=0$, which is
suitable for large $\gamma(t-t_{0})$ when the long-wavelength mode dominates.

$\displaystyle e^{2U_{K}(t,0)}\approx e^{2c\gamma t-2\xi_{0}^{2}\gamma
tK^{2}/\kappa}$ (S120)

Here
$\kappa=\Delta(0)=\sqrt{1+4(|\alpha|-2\beta)\gamma^{2}/\gamma_{0}^{2}},c=a_{1}(0)=2^{-1}(\kappa-1)\gamma_{0}^{2}/\gamma^{2}$.

We apply the approximation used in the derivation of Eq. (S15), in the $\gamma
t\gg 1/|\alpha|$ limit, we integrate over $\bm{K}$ with the integrand
$e^{-2\xi_{0}^{2}\gamma tK^{2}/b}e^{i\bm{K}\cdot\bm{R}}$ and omit the
$K$-dependence of other factors in Eq. (S119), i.e.

$\displaystyle C_{2}(\bm{R},t)$
$\displaystyle\approx\frac{T\sqrt{|\alpha|}}{2E_{\text{c}}\xi_{0}}\frac{1}{\kappa^{2}c}e^{2c\gamma
t}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}e^{-2\xi_{0}^{2}\gamma
tK^{2}/\kappa}e^{i\bm{K}\cdot\bm{R}}$ $\displaystyle=\frac{T}{E_{\rm
c}\xi_{0}^{D}}\frac{\sqrt{|\alpha|}}{2\kappa^{2}c}\frac{e^{2c\gamma
t}}{\left(\sqrt{8\pi\gamma
t/\kappa}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma t/\kappa}}$
$\displaystyle=\frac{\zeta|\alpha|}{2\kappa^{2}c/|\alpha|}\frac{e^{2c\gamma
t}}{\left(\sqrt{8\pi|\alpha|\gamma
t/\kappa}\right)^{D-1}}e^{-\frac{R^{2}}{8\xi_{0}^{2}\gamma t/\kappa}}.$ (S121)

$C_{1}$ is given by the following integration.

$\displaystyle C_{1}(\bm{R},t)$ $\displaystyle=e^{2c\gamma
t}\int\frac{d^{D-1}\bm{K}}{(2\pi)^{D-1}}l^{D-1}D_{K}(0)e^{-2\xi_{0}^{2}\gamma
tK^{2}/\kappa}e^{i\bm{K}\cdot\bm{R}}$ (S122)

Note that Eq. (S121, S122) are written by the notation used in Eq. (S133) with
the bar symbol on the top of the characters being left out. They could be
rescaled back to the notation used in Eq. (S99) to give Eq. (S102, S103).

##  V Coupled three-component field in 3D systems.

###  V.1 Free energy and first-stage dynamics

As mentioned in the last discussion part in the main text, for a theory of
coupled three-component field in 3D systems, similar dynamics could happen.
For example, we consider the following free energy and “model A” dynamics

$\displaystyle F=E_{\text{c}}\int d^{3}\bm{r}\ \
f_{0}(\psi)+f_{1}(\Psi,\Psi^{*})+2g|\Psi|^{2}\psi^{2}.$ (S123)
$\displaystyle\frac{1}{\gamma}\partial_{t}\Psi=-\frac{1}{E_{\text{c}}}\frac{\delta
F}{\delta\Psi^{*}}+\eta$
$\displaystyle\frac{1}{\gamma^{\prime}}\partial_{t}\psi=-\frac{1}{2E_{\text{c}}}\frac{\delta
F}{\delta\psi}+\eta^{\prime}$
$\displaystyle\langle\eta_{i}(\bm{r},t)\eta_{i^{\prime}}(\bm{r^{\prime}},t^{\prime})\rangle=\frac{2T}{\gamma
E_{\text{c}}}\delta_{ii^{\prime}}\delta(\bm{r}-\bm{r^{\prime}})\delta(t-t^{\prime})$
$\displaystyle\langle\eta^{\prime}(\bm{r},t)\eta^{\prime}(\bm{r^{\prime}},t^{\prime})\rangle=\frac{2T}{\gamma^{\prime}E_{\text{c}}}\delta(\bm{r}-\bm{r^{\prime}})\delta(t-t^{\prime})$

with

$\displaystyle f_{0}(\psi_{1})$
$\displaystyle=\xi_{0}^{2}|\nabla\psi|^{2}+(\psi^{2}+\alpha)^{2}$
$\displaystyle f_{1}(\Psi)$
$\displaystyle=\xi_{1}^{2}|\nabla\Psi|^{2}+(|\Psi|^{2}+\chi)^{2}.$

Here $\alpha<\chi<0$, $g>0$. $\psi$ is real. and $\Psi$ is complex with
$\Psi=\Psi_{1}+i\Psi_{2}$.

The theory describes a $Z_{2}$ real-scalar field competing with a $U(1)$
complex-scalar field. Physically, it could effectively describe some CDW-
superconductor competing systems [18, 19]. In this case, $\psi$ represents the
amplitude of the commensurate CDW order and its phase fluctuations can be
ignored. The CDW free energy in Eq. (S2) with $\beta\gg|\alpha|$ belongs to
this case. $\Psi$ can represent the order parameter of the superconducting
order, and $f_{1}(\Psi)$ has the $U(1)$ symmetry $\Psi\rightarrow\Psi
e^{i\theta}$.

The symmetry of the theory is $Z_{2}\otimes U(1)$. Its mean-field ground state
reads $\psi=\pm\sqrt{|\alpha|}$. By the optical pump, a domain wall
configuration $\psi_{0}(\bm{r})=\sqrt{|\alpha|}\tanh\frac{z}{\xi}$ could be
generated by the mechanism mentioned in Sec. I. Similar to the analysis
applied in the first-stage dynamics, we expand the fluctuation of $\Psi$ to
the linear term and obtain

$\displaystyle\frac{1}{\gamma}\partial_{t}\Psi=-(\hat{M}-\xi_{1}^{2}\nabla_{\bm{R}}^{2})\Psi+\eta.$
(S125)

with

$\displaystyle\hat{M}$
$\displaystyle=-\xi_{1}^{2}\partial_{z}^{2}+2\chi+2g\psi_{0}^{2}$
$\displaystyle=-\xi_{1}^{2}\partial_{z}^{2}-2g|\alpha|\frac{1}{\cosh^{2}(z/\xi)}+2g|\alpha|+2\chi$
$\displaystyle=\xi_{1}^{2}[-\partial_{z}^{2}-\frac{1}{\xi^{2}}\frac{2g\xi_{0}^{2}/\xi_{1}^{2}}{\cosh^{2}(x/\xi)}]+2g|\alpha|+2\chi.$
(S126)

Based on the result of the eigenvalue problem in Eq. (S5), the lowest
eigenvalue of $\hat{M}$, which is denoted by $m_{0}$, reads

$\displaystyle
m_{0}=-\frac{\xi_{1}^{2}}{\xi^{2}}(\frac{\sqrt{1+8g\xi_{0}^{2}/\xi_{1}^{2}}-1}{2})^{2}+2g|\alpha|+2\chi.$
(S127)

If $m_{0}<0$, from Eq. (S125), the domain wall is unstable and the first-stage
dynamics, the growth of $\Psi$ around the interface, could happen.

###  V.2 The second-stage dynamics

After the first-stage dynamics, topological defects are generated and a
second-stage dynamics will happen. Basically, the following two situations may
appear:

1\. Only one component of $\Psi$ grows up around the interface. In that case,
the systems go back to the case of the two-component parameter theory we
discuss above. The defects are string defects. The coarsening should be
diffusive.

2\. Both two components of $\Psi$ grows up around the interface. In that case,
the defects are point defects called monopoles.

For the last case, we present a brief discussion of the monopole spatial
profile and the corresponding coarsening dynamics. We separate the coordinate
by $\bm{r}=(z,\bm{R})=(z,R\cos\phi,R\sin\phi)$ and write $\Psi=\rho
e^{i\theta}$. Suppose a monopole solution with a rotational symmetry around
the line of $\bm{R}=0$, i.e. $\rho(\bm{R}=0)=0$. The monopole solution
satisfies

$\displaystyle\frac{\delta F}{\delta\psi}$
$\displaystyle=-2\xi_{0}^{2}\nabla^{2}\psi+4\alpha\psi+4\psi^{3}+4g\rho^{2}\psi=0$
$\displaystyle\frac{\delta F}{\delta\rho}$
$\displaystyle=-2\xi_{1}^{2}\nabla^{2}\rho+2\xi_{1}^{2}\rho|\nabla\theta|^{2}+4\chi\rho+4\rho^{3}+4g\psi^{2}\rho=0$
$\displaystyle\frac{\delta F}{\delta\theta}$
$\displaystyle=-2\xi_{1}^{2}\rho^{2}\nabla^{2}\theta=0,$ (S128)

with boundary conditions

$\displaystyle\lim_{z\rightarrow\pm\infty}\psi$
$\displaystyle=\pm\sqrt{|\alpha|}$
$\displaystyle\lim_{z\rightarrow\pm\infty}\rho$ $\displaystyle=0$
$\displaystyle\oint d\bm{R}\ \ \theta(z,\bm{R})$ $\displaystyle=\pm 2\pi.$
(S129)

In the asymptotic regions with $|\bm{R}-\bm{R}_{0}|\rightarrow+\infty$, the
following solution

$\displaystyle\rho(\bm{r})$ $\displaystyle=\rho_{0}(z),$ $\displaystyle\theta$
$\displaystyle=\pm\phi,$ $\displaystyle\psi(\bm{r})$
$\displaystyle=\psi_{0}(z),$ (S130)

satisfies Eq. (S128, S129), as long as the the solution of $\psi_{0}(z)$ and
$\rho_{0}(z)$ with respect to the following equation exists.

$\displaystyle-2\xi_{0}^{2}\partial_{z}^{2}\psi_{0}+4\alpha\psi_{0}+4\psi_{0}^{3}+4g\rho_{0}^{2}\psi_{0}$
$\displaystyle=0,$
$\displaystyle-2\xi_{1}^{2}\partial_{z}^{2}\rho_{0}+4\chi\rho_{0}+4\rho_{0}^{3}+4g\psi_{0}^{2}\rho_{0}$
$\displaystyle=0,$ (S131)

with the boundary condition
$\lim_{z\rightarrow\pm\infty}\psi_{0}(z)=\pm\sqrt{|\alpha|}$,
$\lim_{z\rightarrow\pm\infty}\rho_{0}(z)=0$.
Substituting $\rho_{0}(z)$, $\psi_{0}(z)$ into Eq. (S123), we can separate the
free energy by

$\displaystyle F=E_{\rm d}+E_{0}$

with

$\displaystyle E_{\rm d}$ $\displaystyle=E_{\text{c}}\int d^{2}\bm{R}\ \
\xi_{1}^{2}\Gamma|\nabla_{\bm{R}}\theta|^{2},\ \ \ \ \ \Gamma=\int dz\ \
\rho_{0}^{2}(z)$ $\displaystyle E_{0}$ $\displaystyle=E_{\text{c}}\int
d^{2}\bm{R}\int dz\ \
\xi_{0}^{2}(\psi_{0}^{\prime}(z))^{2}+(\psi_{0}^{2}+\alpha)^{2}+\xi_{1}^{2}(\rho_{0}^{\prime}(z))^{2}+(\rho_{0}^{2}+\chi)^{2}+2g\rho_{0}^{2}\psi_{0}^{2}.$
(S132)

Here $E_{0}$ is the energy of the solution $\rho_{0}(z)$ and $\psi_{0}(z)$,
$E_{\rm d}$ is the energy per monopole, which has the same form as the energy
per vortex of a 2D system with the $U(1)$ symmetry. Therefore, in this case,
the coarsening dynamics has the same scaling behavior as the 2D vortex system,
where the typical length scale $L(t)$ scales as
$L(t)\sim\xi_{0}\sqrt{\frac{\gamma t}{\ln(|\alpha|\gamma t)}}$ [8].

##  VI Details of the numerical simulation

In Fig. 3 in the main text, the result of the numerical simulation of the
scaling of $L(t)$ in the second stage dynamics was shown. In the simulation,
the “model A” dynamics were rescaled in the same way as in Eq. (S100, S101)
with $b=\sqrt{|\alpha|}$, $\bar{\alpha}=-1$. We further define
$\tau=|\alpha|\gamma t$, $\tilde{\bm{r}}=\bm{r}/\xi$ and
$\tilde{\nabla}=\partial/\partial\tilde{\bm{r}}$ and write the dynamical
equation in the following form;

$\displaystyle\frac{\partial\bar{\psi}}{\partial\tau}=\tilde{\nabla}^{2}\bar{\psi}+2\bar{\psi}-2\frac{\beta}{|\alpha|}i\bar{\psi_{2}}-2|\bar{\psi}|^{2}\bar{\psi}+\bar{\eta},$
(S133)

Here $\bar{\psi}=\psi/\sqrt{|\alpha|}$, $\bar{\eta}=\eta|\alpha|^{-3/2}$. The
rescaled noise correaltor reads
$\langle\bar{\eta}_{i}(\tilde{\bm{r}},\tau)\bar{\eta}_{i^{\prime}}(\tilde{\bm{r}}^{\prime},\tau^{\prime})\rangle=2\zeta\delta_{ii^{\prime}}\delta(\tilde{\bm{r}}-\tilde{\bm{r}}^{\prime})\delta(\tau-\tau^{\prime})$.
In this equation, there are only two dimensionless parameters, $\zeta$ and
$\beta/|\alpha|$. We discretize the space and time by the steps
$d\tilde{r}_{j=x,y,z}=0.5$ and $d\tau=0.03$ in the program and use the
following Euler method to update $\bar{\psi}$ in the time evolution;

$\displaystyle\bar{\psi}(\tilde{\bm{r}},\tau+d\tau)=\Big{(}\tilde{\nabla}^{2}\bar{\psi}(\tilde{\bm{r}},\tau)+2\bar{\psi}(\tilde{\bm{r}},\tau)-2\frac{\beta}{|\alpha|}i\bar{\psi_{2}}(\tilde{\bm{r}},\tau)-2|\bar{\psi}(\tilde{\bm{r}},\tau)|^{2}\bar{\psi}(\tilde{\bm{r}},\tau)\Big{)}d\tau+\bar{\eta}(\tilde{\bm{r}},\tau)d\tau.$
(S134)

Here $\bar{\eta}(\tilde{\bm{r}},\tau)d\tau$ are generated by
$\bar{\eta}_{i}(\tilde{\bm{r}},\tau)d\tau=\sqrt{\frac{2\zeta
d\tau}{\Pi_{j}dr_{j}}}n_{i}(\tilde{r},\tau)$ with $2\zeta=0.0001$ with
$n_{i}(\tilde{r},\tau)$ being generated from independent standard normal
distribution for each discretized $\tilde{\bm{r}}$ and $\tau$.

The initial condition of Eq. (S133) in the program is simply set as
$\bar{\psi}(\tilde{\bm{r}},\tau=0)=\tanh{\tilde{z}}$. Note that the initial
fluctuation $\delta\bar{\psi}(\tilde{\bm{r}})$ is neglected in the program
because it is hard to generate the initial fluctuation, which is a random
function of $\bm{r}$ satisfying a specific statistical correlation function.
The simplification is appropriate since the simulation is to fit the scaling
of $L(t)$ in the second-stage dynamics, which is unrelated to the initial
condition information.

The following method evaluates the length scale $L(t)$ in the simulation. In
2D systems, we count a number $N$ of vortices in the $y$-$z$ plane and define
$L=L_{y}/N$. In 3D systems, we count a number $N$ of vortices in the $y$-$z$
planes for each discrete $x$, average $N$ over different $x$, and define
$L=L_{y}/\langle N\rangle$, where $L_{y}$ is the system size along $y$. To
count the number of vortices, we extract the function
$P(y)=1-|\bar{\psi}(y,z=0)|^{2}$. A peak of $P(y)$ with a height larger than
0.6 is identified as a vortex(or antivortex). $N$ equals the number of peaks.

Due to the random noise, the fitted value of
$p=\frac{\log_{10}(L/\xi)}{\log_{10}(|\alpha|\gamma t)}$ gets a statistical
error for different runs of the program. Table. 1 lists the results of the
fitted value of $p$ for ten program runs, which shows that the statistical
error is small.

$D=3$. Run number: | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$\beta/|\alpha|=0$ | 0.45365 | 0.45123 | 0.45414 | 0.44534 | 0.45779 | 0.46566 | 0.46025 | 0.44490 | 0.44960 | 0.45672
$\beta/|\alpha|=0.1$ | 0.45910 | 0.44516 | 0.45337 | 0.46073 | 0.45576 | 0.45869 | 0.46217 | 0.45613 | 0.45594 | 0.44143

$D=2$. Run number: | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$\beta/|\alpha|=0.001$ | 0.27497 | 0.26555 | 0.30045 | 0.30077 | 0.27175 | 0.29875 | 0.28527 | 0.27797 | 0.28693 | 0.27962
$\beta/|\alpha|=0.1$ | 0.19070 | 0.16815 | 0.18172 | 0.17463 | 0.17473 | 0.16713 | 0.17261 | 0.16296 | 0.18281 | 0.17346

Table 1:  Linear fit results of the values of
$p=\frac{\log_{10}(L/\xi)}{\log_{10}(|\alpha|\gamma t)}$ for ten runs of the
program in $D=3$ and $D=2$. The parameter set is the same as mentioned in the
caption of Fig. (3) in the main text.

## References

  * Yusupov  _et al._ [2010] R. Yusupov, T. Mertelj, V. V. Kabanov, S. Brazovskii, P. Kusar, J.-H. Chu, I. R. Fisher, and D. Mihailovic, Coherent dynamics of macroscopic electronic order through a symmetry breaking transition, Nature Physics 6, 681 (2010).
  * 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._ , Optical manipulation of electronic dimensionality in a quantum material, Nature 595, 239 (2021).
  * 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, Ultrafast formation of domain walls of a charge density wave in SmTe3, Phys. Rev. B 103, 054109 (2021).
  * Flügge [2012] S. Flügge,  _Practical quantum mechanics_ (Springer Science & Business Media, 2012).
  * Sun and Millis [2020] Z. Sun and A. J. Millis, Transient trapping into metastable states in systems with competing orders, Phys. Rev. X 10, 021028 (2020).
  * Cardy [1996] J. Cardy,  _Scaling and Renormalization in Statistical Physics_, Cambridge Lecture Notes in Physics (Cambridge University Press, 1996).
  * Kardar [2007] M. Kardar,  _Statistical physics of fields_ (Cambridge University Press, 2007).
  * Bray [2002] A. J. Bray, Theory of phase-ordering kinetics, Advances in Physics 51, 481 (2002).
  * Eto  _et al._ [2023] M. Eto, Y. Hamada, and M. Nitta, Composite topological solitons consisting of domain walls, strings, and monopoles in O(N) models, Journal of High Energy Physics 2023, 150 (2023).
  * Yurke  _et al._ [1993] B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse, Coarsening dynamics of the XY model, Phys. Rev. E 47, 1525 (1993).
  * Pargellis  _et al._ [1992] A. N. Pargellis, P. Finn, J. W. Goodby, P. Panizza, B. Yurke, and P. E. Cladis, Defect dynamics and coarsening dynamics in smectic-C films, Phys. Rev. A 46, 7765 (1992).
  * Mondello and Goldenfeld [1990] M. Mondello and N. Goldenfeld, Scaling and vortex dynamics after the quench of a system with a continuous symmetry, Phys. Rev. A 42, 5865 (1990).
  * Toyoki and Honda [1987] H. Toyoki and K. Honda, Ordering Dynamics of a Deeply Quenched Complex Field, Progress of Theoretical Physics 78, 237 (1987).
  * Mondello and Goldenfeld [1992] M. Mondello and N. Goldenfeld, Scaling and vortex-string dynamics in a three-dimensional system with a continuous symmetry, Phys. Rev. A 45, 657 (1992).
  * Duan and Leggett [1992] J.-M. Duan and A. J. Leggett, Inertial mass of a moving singularity in a fermi superfluid, Phys. Rev. Lett. 68, 1216 (1992).
  * Lemonik and Mitra [2018] Y. Lemonik and A. Mitra, Quench dynamics of superconducting fluctuations and optical conductivity in a disordered system, Phys. Rev. B 98, 214514 (2018).
  * Sun  _et al._ [2020] Z. Sun, M. M. Fogler, D. N. Basov, and A. J. Millis, Collective modes and terahertz near-field response of superconductors, Phys. Rev. Res. 2, 023413 (2020).
  * Yan  _et al._ [2017] S. Yan, D. Iaia, E. Morosan, E. Fradkin, P. Abbamonte, and V. Madhavan, Influence of domain walls in the incommensurate charge density wave state of cu intercalated 1T-TiSe2, Phys. Rev. Lett. 118, 106405 (2017).
  * Joe  _et al._ [2014] Y. I. Joe, X. Chen, P. Ghaemi, K. Finkelstein, G. de La Peña, Y. Gan, J. Lee, S. Yuan, J. Geck, G. MacDougall, _et al._ , Emergence of charge density wave domain walls above the superconducting dome in 1T-TiSe2, Nature Physics 10, 421 (2014).