# Cosmological scaling of precursor domain walls

Mainak Mukhopadhyay ID<EMAIL_ADDRESS>Department of Physics; Department of
Astronomy & Astrophysics; Center for Multimessenger Astrophysics, Institute
for Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA  Oriol Pujolas ID<EMAIL_ADDRESS>Institut de Física
d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology
(BIST), Campus UAB, 08193 Bellaterra, Barcelona, Spain  George Zahariade ID
<EMAIL_ADDRESS>Instytut Fizyki Teoretycznej, Uniwersytet
Jagielloński, Lojasiewicza 11, 30-348 Kraków, Poland

###### Abstract

Domain wall (DW) networks have a large impact on cosmology and present
interesting dynamics that can be controlled by various scaling regimes. In the
first stage after spontaneous breaking of the discrete symmetry, the network
is seeded with ‘DW precursors’, the zeros of a tachyonic field. At
sufficiently weak coupling, this stage can be quite long. The network is then
driven to a non-relativistic scaling regime: in flat spacetime the correlation
length grows like $L\sim t^{\kappa}$ with $\kappa=1/2$. We focus on the
precursor regime in cosmology, assuming a power-law scale factor $a\propto
t^{\alpha}$. We obtain the scaling exponent as a function of the external
parameter, $\kappa(\alpha)$, by explicit computation in $1+1$ and $2+1$
dimensions, and find a smooth transition from nonrelativistic scaling with
$\kappa\simeq 1/2$ for $\alpha\lesssim 1/2$ to DW gas regime
$\kappa\simeq\alpha$ for $\alpha\gtrsim 1/2$, confirming previous arguments.
The precise form of the transition $\kappa(\alpha)$ is surprisingly
independent of dimension, suggesting that similar results should also be valid
in $3+1$ dimensions.

##  I Introduction

Recent results from pulsar timing arrays [1, 2, 3, 4, 5, 6, 7] reported strong
evidence for the existence of a nanohertz (nHz) stochastic gravitational wave
(GW) background. Although the most widely accepted interpretation of this
evidence is a GW background generated by supermassive black hole binary
systems, other possibilities and sources are reasonably plausible, including
primordial GW backgrounds originating from beyond standard model (BSM)
physics. In particular, a GW background from cosmic domain walls (DWs) could
be an alternative interpretation of the signal [8, 9, 10, 11, 12, 13, 14, 15].

DWs are codimension one topological defects that are generically formed as a
result of the spontaneous breaking of a discrete symmetry, see [16, 17, 18,
19] for reviews. They are stable solutions of classical field theories with a
discrete vacuum manifold, that smoothly interpolate between spatial regions
with different vacua. Various BSM scenarios, including QCD axion models [20,
21], lead to the formation of such defects in the early universe. The energy
density stored in such objects decays more slowly than matter leading to a so-
called DW problem [22, 20, 17]. There are natural mechanisms to avoid DW
domination so the problem in reality is a virtue: the tendency of DWs to
dominate the energy density implies that models with DWs leave a large impact
on cosmology. Phenomenologically, this can be in the form of primordial black
holes, GWs, and/or other relics.

In order to confront the DW models with observations it is important to know
in detail their dynamics, abundance etc. As is well known, the evolution of a
DW network is dominated by a self-similar ‘scaling’ attractor, characterized
by having about one Hubble-sized DW per Hubble volume at all times. Most of
the previous work on DW evolution result from numerical simulations in
classical field theory simulations [23, 24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34, 35, 15]. In the analytical front, results derived using the thin wall
limit (combined with a mean field approximation) [36, 37] are able to
reproduce some scaling properties of DW networks in cosmology.

Further analytic progress was made in [38, 39] and [40] by addressing the
problem using quantum field theory methods. The formation of a DW network upon
$\mathbb{Z}_{2}$ symmetry breaking transition was analyzed exploiting the DW
precursor limit, where the symmetry is broken but the field hasn’t stabilized
yet. This already leads to well defined dynamical DW objects (the seeds of the
actual DWs formed when the field is nonliear) that arise in the linear theory
and so their behaviour can be computed.

Interestingly, a network of precursor DWs formed after the $\mathbb{Z}_{2}$
quench reaches a scaling regime characterized by a correlation length that
grows in time like $L\sim t^{1/2}$ [38, 39], see also [40]. These works
assumed the DW precursor limit and a Minkowski background. In this paper we
aim to relax the latter assumption and extend the previous analysis to an
expanding spacetime more directly relevant to cosmology.

In Ref. [40], a simple argument was given on how the $t^{1/2}$ scaling of
precursor DW networks extends to cosmology. The origin of the $t^{1/2}$
scaling resides in the fact that precursor DWs are non-relativistic at late
times. This is expected to hold in a slowly expanding background, and the
scaling should continue to be $t^{1/2}$ in terms of cosmic time. However, this
cannot remain true in cosmologies with a sufficiently fast expansion rate. For
power law cosmologies, $a(t)\sim t^{\alpha}$, with $\alpha\gtrsim 1/2$ one
expects a change of behaviour to a DW gas regime where $L(t)\propto a(t)$.
From this argument one cannot tell whether this transition is sharp, occurring
at the precise value $\alpha=1/2$, or rather smooth, and to what extent. Also,
$\alpha=1/2$ should mark the onset transition independently of $d$, the number
of space dimensions, although the details of the transition might depend on
$d$. The purpose of this work is to exhibit this transition in detail by
direct calculation.

For simplcity we first work in $1+1$ dimensions, and thus limit ourselves to
studying the evolution of the number density of kinks (and antikinks)
following a $\mathbb{Z}_{2}$-symmetry breaking phase transition in an
expanding background. We assume a spatially flat Friedmann-Lemaître-Robertson-
Walker (FLRW) background metric with a scale factor that grows as a power law
with respect to cosmological time. In particular, we aim to understand how the
growth of the scale factor affects the late time scaling of the kink number
density. We then extend the results to domain walls in $2+1$ dimensions.

The paper is organized as follows. We begin by reviewing the formalism and
setting up the main equations for the $1+1$ dimensional case in Sec. II (with
more details included in the appendices). The extension to the $2+1$
dimensional case is briefly outlined in Sec. III. Our main results are
showcased in Sec. IV. We discuss the effective dynamics of precursor DWs in
the velocity-one-scale (VOS) model in Sec. V. We conclude with a brief
discussion in Sec. VI.

##  II Kinks in an expanding background

In this section we extend the formalism of Refs. [38, 39, 40] to a $1+1$
dimensional expanding spacetime. As usual, dotted quantities denote time
derivatives while primed quantities denote spatial derivatives. We work in
units where $\hbar=c=1$.

The aim is to investigate the dynamics of the average number density of kinks
formed as a scalar field $\phi(t,x)$ undergoes a quantum phase transition111We
consider a quantum phase transition implying a zero temperature vacuum state
for the field. However, we expect that our results on the late time scaling
would also apply to a thermal transition [41, 42, 43, 44]. This is because the
precursor DWs are over-damped and hence their late time behaviour follows non-
relativistic scaling [40] (refer to Sec. V for more details). . The background
metric is taken to be $ds^{2}=dt^{2}-a(t)^{2}dx^{2}$, with scale factor

$a(t)\propto t^{\alpha}\,,$ (1)

where $0\leq\alpha<1$ to focus on the case of decelerating expansion. The
action for the scalar field is thus given by

$S=\int
dt\,dx\,a(t)\left[\frac{1}{2}\dot{\phi}^{2}-\frac{1}{2}a(t)^{-2}\phi^{\prime
2}-\frac{1}{2}m^{2}(t)\phi^{2}+\cdots\right],$ (2)

where

$m^{2}(t)=-m^{2}{\rm tanh}\bigg{(}\frac{t-t_{\rm pt}}{\tau}\bigg{)},$ (3)

and the dots denote stabilizing terms for the potential e.g.
$\lambda\phi^{4}$. In this work (as in [38, 39, 40]) we neglect such terms by
focusing on the so-called spinodal instability phase which can last for
arbitrary long time if the self-coupling $\lambda$ is weak enough. The phase
transition is triggered at $t=t_{\rm pt}>0$ when the mass squared becomes
negative and the potential turns tachyonic. The quench time scale $\tau<t_{\rm
pt}$ parametrizes the rapidity of the phase transition. Notice that we expect
kinks (and antikinks) to form as a consequence of this phase transition, but
since we neglect self-coupling terms we can only strictly talk about kink (and
antikink) precursors. Therefore, and in order to not overburden the reader, we
will assume the word precursor to be implied even when not explicitly stated.

Let us now summarize the task at hand. We want to determine the average number
density of kinks formed during the phase transition as a function of time,
$n_{K}(t)$. This roughly boils down to counting _zeros_ of the scalar field
$\phi$. Since we are interested in finding the quantum dynamics of $n_{K}(t)$
we will need to compute the quantum average of a suitably constructed number
density of zeros operator in a state $|\Psi\rangle$ of the quantized field
$\phi$. We will work in the Schrödinger picture and choose the state
$|\Psi(t)\rangle$ to coincide with the zeroth order adiabatic vacuum state of
the quadratic theory at some time $t_{0}$ verifying $0<t_{0}<t_{\rm pt}$. This
will depend only on the values $m^{2}(t_{0})$ and $a_{0}\equiv a(t_{0})$. If
$t_{\rm pt}-t_{0}\gg\tau$, our physical results will not depend on the choice
of initial time since $m^{2}(t_{0})\approx m^{2}$ does not depend on $t_{0}$
while $a_{0}$ only serves to rescale $L$ and thus has no impact in the
infinite volume limit.

By discretizing the field theory on a periodic lattice of size $L$ with $N$
evenly spaced points (such that the lattice spacing is given by $\ell\equiv
L/N$) and deploying the methods of Refs. [38, 39, 40] (see Appendix A) we find
the probability density functional for a given discretized field configuration
$\bm{\phi}\equiv(\phi_{1},\phi_{2},\dots,\phi_{N})^{T}$ at time $t$ to
be222Here we use the notation $\phi_{i}\equiv\phi(t,i\ell)$ for the
discretized field values, and we denote the transpose using T.

$\mathcal{P}(t;\bm{\phi})=\frac{1}{\sqrt{{\rm
det}\big{(}2\pi\bm{K}(t)\big{)}}}e^{-\bm{\phi}^{T}\bm{K}(t)^{-1}\bm{\phi}/2},$
(4)

where the covariance matrix $\bm{K}(t)$ is given by

$[\bm{K}(t)]_{jk}=\frac{1}{N}\sum_{n=0}^{N-1}|c_{n}(t)|^{2}{\rm
cos}\big{(}2\pi n(j-k)/N\big{)},$ (5)

and the mode functions $c_{n}(t)$ (for $0\leq n<N$) obey

$\ddot{c}_{n}+H(t)\dot{c}_{n}+\left[\frac{4}{(a(t)\ell)^{2}}{\rm
sin}^{2}\left(\frac{\pi n}{N}\right)+m^{2}(t)\right]c_{n}=0,$ (6)

$H\equiv\dot{a}/a$ being the Hubble rate. Our particular choice of quantum
state imposes the specific initial conditions

$\displaystyle c_{n}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2a_{0}\ell}}\left[\frac{4}{(a_{0}\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}(t_{0})\right]^{-1/4},$ (7)
$\displaystyle{\dot{c}}_{n}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{i}{\sqrt{2a_{0}\ell}}\left[\frac{4}{(a_{0}\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}(t_{0})\right]^{1/4}.$ (8)

Since we effectively work in a regime where the Gaussian approximation holds,
the mode functions are seen to contain all the information about the quantum
dynamics of the system.

The quantum operator used to estimate the average kink number density is the
direct generalization of the flat spacetime one described at length in [38,
39, 40]. Essentially, the idea is to look for kinks and antikinks among the
zeros of the field configuration. Of course quantum fluctuations introduce
spurious zeros that need to be discarded and we will do so through a
prescription which will be briefly described later on. We define a number
density of zeros operator by

$n_{\rm
Z}\left(\hat{\bm{\phi}}\right)\equiv\frac{1}{a(t)L}\sum_{i=1}^{N}\frac{1}{4}\left[{\rm
sgn}\left(\hat{\phi}_{i}\right)-{\rm
sgn}\left(\hat{\phi}_{i+1}\right)\right]^{2},$ (9)

where the symbol ${\rm sgn}$ denotes the signum function. This counts the
number density of sign changes of the $\phi$ field profile within the
(comoving) size $a(t)L$ of the lattice, to a level of accuracy determined by
its coarseness $\ell$.

Using the translational invariance of the problem, the vacuum expectation
value of this operator in the state $\left|\Psi\right\rangle$ can readily be
computed (see [38, 39, 40] for details). It reads

$\displaystyle\left\langle n_{\rm Z}\left(\hat{\bm{\phi}}\right)\right\rangle$
$\displaystyle\equiv$ $\displaystyle\int d\phi_{1}\cdots d\phi_{n}\,n_{\rm
Z}\left(\hat{\bm{\phi}}\right)\mathcal{P}(t;\bm{\phi})$ (10) $\displaystyle=$
$\displaystyle\frac{N}{\pi a(t)L}\cos^{-1}\left(\frac{\beta}{\alpha}\right),$

where $\alpha$ and $\beta$ are given by

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}|c_{n}|^{2},$ (11) $\displaystyle\beta$
$\displaystyle=$ $\displaystyle\frac{1}{N}\sum_{n=1}^{N}|c_{n}|^{2}\cos(2\pi
n/N).$ (12)

As we have already mentioned, we expect this quantity to not be an accurate
estimator of the number density of kinks and antikinks (even for very fine
lattices) because of the random fluctuations of the field caused by
oscillating modes. This can remediated by restricting the sums in Eqs. (11)
and (12) to only those modes $c_{n}(t)$ that are _non-oscillating_ or in other
words _unstable_. Such modes verify

$\omega_{n}^{2}(t)\equiv\frac{4}{(a(t)\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}(t)\leq 0.$ (13)

Before the phase transition, that is, for $t<t_{\rm pt}$, all the modes are
stable and there are no kinks or antikinks. For $t>t_{\rm pt}$ the time-
dependent critical value of $n$, $n_{c}(t)$, separating the stable modes from
the unstable ones is given by

$n_{c}(t)\equiv\left\lfloor\frac{N}{\pi}\sin^{-1}\left(\frac{a(t)\ell\sqrt{|m^{2}(t)|}}{2}\right)\right\rfloor,$
(14)

where $\lfloor\cdot\rfloor$ denotes the integer part function. Thus we find
the average number density of kinks $n_{K}$ to be

$n_{K}=\frac{N}{\pi a(t)L}{\rm
cos}^{-1}\left(\frac{\bar{\beta}}{\bar{\alpha}}\right),$ (15)

where

$\displaystyle\bar{\alpha}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\sum_{\omega_{n}^{2}(t)\leq
0}|c_{n}|^{2}=\frac{1}{N}\left[|c_{0}|^{2}+2\sum_{n=1}^{n_{c}(t)}|c_{n}|^{2}\right],$
(16) $\displaystyle\bar{\beta}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N}\sum_{\omega_{n}^{2}(t)\leq 0}|c_{n}|^{2}\cos(2\pi
n/N)$ (17) $\displaystyle=$
$\displaystyle\frac{1}{N}\left[|c_{0}|^{2}+2\sum_{n=1}^{n_{c}(t)}|c_{n}|^{2}\cos(2\pi
n/N)\right].$

Notice that the number of unstable modes grows with time (after the phase
transition). This is due to the physical mode wavelengths getting stretched by
the (cosmological) expansion of the lattice or in other words to UV modes
becoming IR. In fact, for fixed lattice spacing $\ell$, after a time of order
$t/t_{0}\sim(a_{0}\ell)^{-1/\alpha}m^{-1/\alpha}$ all modes are unstable and
the estimate of $n_{K}(t)$ ceases to be accurate: the longer we want to track
its dynamics, the more modes we will need to take into consideration (and the
finer our lattice needs to be).

In general, Eq. (6) will need to be solved numerically allowing us to compute
$\bar{\alpha}$, $\bar{\beta}$, and finally $n_{K}$ via Eq. (15). Our results
will be showcased in Sec. IV and we leave the more precise algorithmic details
to Appendix B.1.

Before ending this section, let us quickly mention that although it is not
possible to estimate $n_{K}$ analytically in general, or for arbitrary rates
of expansion $\alpha$ and phase transition time scales $\tau$, the case of the
sudden phase transition $\tau=0$ constitutes a notable exception. In this
scenario $m^{2}(t)=-m^{2}\Theta(t-t_{\rm pt})$, where $\Theta$ denotes the
Heaviside step function, and one can choose $t_{0}=t_{\rm pt}^{-}$. Equation
(6) then reduces to

$\ddot{c}_{n}+\frac{\alpha}{t}\dot{c}_{n}+\left[\frac{4}{(a_{0}\ell)^{2}(t/t_{\rm
pt})^{2\alpha}}\sin^{2}\left(\frac{\pi n}{N}\right)-m^{2}\right]c_{n}=0,$ (18)

with initial conditions,

$\displaystyle c_{n}(t_{\rm pt})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2a_{0}\ell}}\left[\frac{4}{(a_{0}\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}\right]^{-1/4},$ (19) $\displaystyle{\dot{c}}_{n}(t_{\rm
pt})$ $\displaystyle=$
$\displaystyle\frac{i}{\sqrt{2a_{0}\ell}}\left[\frac{4}{(a_{0}\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}\right]^{1/4}.$ (20)

It is still non-trivial to find a closed form expression for $c_{n}(t)$ for
general $\alpha$. However, the particular cases of $\alpha=0$ and $1$ can be
solved analytically in the continuum ($\ell\rightarrow 0$) and infinite volume
($L\rightarrow\infty$) limit333In the $\alpha=1/2$ case, Eq. (18) takes the
form of a hypergeometric differential equation which can be solved
analytically as well but it is no easy task to use the resulting closed-form
expressions for a useful estimate of $n_{K}(t)$., thus giving access to closed
form expressions for $n_{K}$ and providing in principle a useful benchmark
against which to gauge the accuracy of our numerics. The case of Minkowski
spacetime, $\alpha=0$, has already been studied in detail in Refs. [38, 39,
40] and it was found that $n_{K}(t)\sim t^{-1/2}$ for late times. In the case
of linear expansion, $\alpha=1$, Eq. (18) takes the form of a Bessel
differential equation whose solutions are given in terms of modified Bessel
functions $I_{n}$ and $K_{n}$,

$\displaystyle c_{n}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{2\ell}(k_{n}^{2}+m^{2})^{1/4}}\left[I_{ik_{n}}(mt)\left(m\left(K_{ik_{n}-1}(m)+\right.\right.\right.$
(21)
$\displaystyle\left.\left.\left.K_{ik_{n}+1}(m)\right)+2i\sqrt{k_{n}^{2}+m^{2}}K_{ik_{n}}(m)\right)+2K_{ik_{n}}(mt)\right.$
$\displaystyle\left.\left(mI_{ik_{n}+1}(m)-i\left(\sqrt{k_{n}^{2}+m^{2}}-k_{n}\right)I_{ik_{n}}(m)\right)\right],$

where $k_{n}\equiv 2\sin(\pi n/N)/\ell$ and we chose $a_{0}=1$, $t_{\rm pt}=1$
to simplify the expressions. In the late time limit $I_{z}(mt)\sim
e^{mt}/\sqrt{2\pi mt}$ while $K_{z}(mt)\sim\sqrt{2}\,e^{-mt}/\sqrt{\pi mt}$,
for all $z\in\mathbb{C}$ (see e.g. p.203 in [45], (10.40.5) in [46]), so it is
easy to see that the time dependence in $|c_{n}(t)|^{2}$ factors out and is
equal to $e^{2mt}/t$. Therefore the ratio $\bar{\beta}/\bar{\alpha}$ becomes
time independent at large times and $n_{K}(t)\sim 1/a(t)\sim t^{-1}$. This
corresponds to the domain wall gas limit and we will discuss it more in Sec.
IV.

The previous discussion applied to kinks and antikinks in one spatial
dimension. Such objects are expected to behave like point particles and do not
have any internal structure, so we now proceed to briefly extend the formalism
to similar but more complex objects: domain walls in two (and higher) spatial
dimensions.

##  III Domain walls

The formalism outlined in Sec. II for kinks in $1+1$ dimensions can be
trivially extended to study the dynamics of the domain wall length density in
a $2+1$ dimensional theory with a real scalar field $\phi(t,x,y)$ undergoing a
$\mathbb{Z}_{2}$-symmetry breaking phase transition. We begin with the
Lagrangian density given by

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
dt\,dx\,dy\,a^{2}(t)\left[\frac{1}{2}(\partial_{t}\phi)^{2}-\frac{1}{2}a^{-2}(t)(\partial_{x}\phi)^{2}\right.$
(22) $\displaystyle\hskip
28.45274pt\left.-\frac{1}{2}a^{-2}(t)(\partial_{y}\phi)^{2}-\frac{1}{2}m^{2}(t)\phi^{2}+\cdots\right],$

where the definitions of $m^{2}(t)$ and $a(t)$ are unchanged, and we discard
any non-linear stabilizing terms. Repeating the steps outlined in the previous
section (and carefully dealing with formal complications due to the presence
of two spatial indices), it is not hard to convince oneself that the quantum
dynamics of the system is fully determined by the mode functions
$c_{n,n^{\prime}}(t)$ (for $0\leq n,n^{\prime}<N$) verifying

$\displaystyle\ddot{c}_{n,n^{\prime}}$ $\displaystyle+$ $\displaystyle
2H(t)\dot{c}_{n,n^{\prime}}+\left[\frac{4}{(a(t)\ell)^{2}}\left\\{{\rm
sin}^{2}\left(\frac{\pi n}{N}\right)\right.\right.$ (23)
$\displaystyle+\left.\left.{\rm sin}^{2}\left(\frac{\pi
n^{\prime}}{N}\right)\right\\}+m^{2}(t)\right]c_{n,n^{\prime}}=0,$

with initial conditions

$\displaystyle c_{n,n^{\prime}}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{a_{0}\ell\sqrt{2}}\left[\frac{4}{(a_{0}\ell)^{2}}\left\\{\sin^{2}\left(\frac{\pi
n}{N}\right)\right.\right.$ (24)
$\displaystyle+\left.\left.\sin^{2}\left(\frac{\pi
n^{\prime}}{N}\right)\right\\}+m^{2}(t_{0})\right]^{-1/4},$
$\displaystyle{\dot{c}}_{n,n^{\prime}}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{i}{a_{0}\ell\sqrt{2}}\left[\frac{4}{(a_{0}\ell)^{2}}\left\\{\sin^{2}\left(\frac{\pi
n}{N}\right)\right.\right.$ (25)
$\displaystyle+\left.\left.\sin^{2}\left(\frac{\pi
n^{\prime}}{N}\right)\right\\}+m^{2}(t_{0})\right]^{1/4}.$

As mentioned above the relevant quantity that needs to be computed here is the
domain wall length density (in three spatial dimensions this corresponds to
the more familiar domain wall area density). As discussed in detail in Ref.
[40], this is obtained from the quantum average of a length density of zeros
operator444Again we use the notation $\phi_{ij}\equiv\phi(t,i\ell,j\ell)$ for
the discretized field values.

$\displaystyle\mu_{Z}\left(\hat{\bm{\phi}}\right)$ $\displaystyle\equiv$
$\displaystyle\frac{a(t)\ell}{(a(t)L)^{2}}\sum_{i,j=1}^{N}\frac{1}{4}\left[{\rm
sgn}(\hat{\phi}_{ij})-{\rm sgn}(\hat{\phi}_{i+1,j})\right]^{2}$ (26)
$\displaystyle+\frac{a(t)\ell}{(a(t)L)^{2}}\sum_{i,j=1}^{N}\frac{1}{4}\left[{\rm
sgn}(\hat{\phi}_{ij})-{\rm sgn}(\hat{\phi}_{i,j+1})\right]^{2}.$

Notice that, for small enough $\ell$, a domain wall can cross a unit cell in
the $x$ or in the $y$ direction, hence the presence of two distinct terms in
the above definition. While this is of course an idealization, since domain
walls can have arbitrary orientations, we expect the error to be
$\mathcal{O}(1)$ on average.

Using the translational and rotational invariance of the problem we
find555Note that in Ref. [40] there is a typo in Eq. (59), where a factor of
$2$ appears in the denominator of formula for the average DW area density
(instead of the numerator). This however does not alter the main results.

$\left\langle\mu_{Z}\left(\hat{\bm{\phi}}\right)\right\rangle=\frac{2N}{\pi
a(t)L}{\rm cos}^{-1}\left(\frac{{\beta}}{{\alpha}}\right)$ (27)

where

$\displaystyle\alpha(t)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N^{2}}\sum_{n,n^{\prime}=0}^{N-1}|c_{n,n^{\prime}}(t)|^{2},$
(28) $\displaystyle\beta(t)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N^{2}}\sum_{n,n^{\prime}=0}^{N-1}|c_{n,n^{\prime}}(t)|^{2}\cos(2\pi
n^{\prime}/N).$ (29)

The factor $2$ appearing in (27) is a consequence of the higher spatial
dimension. Once again the actual domain wall length density, $\mu_{DW}$, is
obtained from (27) by requiring the sums in $\alpha$, $\beta$ to run only over
those modes that are unstable. Numerical results for the dynamics of the
domain wall length density in two spatial dimensions will be presented in the
next section.

It is also easy to see how the formalism can be extended to higher spatial
dimensions $d$, where the quantum average of the area density of zeros of the
field $\phi$ is obtained by replacing the factor $2$ in (27) by $d$ and
$\alpha$, $\beta$ are defined in an analogous way with the help of the higher
dimensional mode functions.

##  IV Results

We are now in a position to discuss the main results of this work. We begin
with the one spatial dimension case of Sec. II. The dynamics of the kink
number density $n_{K}(t)$ is obtained by determining the mode functions
$c_{n}(t)$ and plugging them in Eq. (15). More precisely we numerically solve
the differential equations (6) with initial conditions given in (7) and
(8)666For numerical reasons we solve a slightly different but equivalent set
of equations, see Appendix B.1 for details. using an explicit Crank-Nicholson
method. The computations were performed using a _parallelized CPU code_.
Working in units where $m=1$ (and unless specified otherwise) we set the time
of the phase transition to $t_{\rm pt}=1$ and start evolving the equations at
$t_{0}=0.5$. To have an accurate estimate of the average kink number density
up to $\mathcal{O}(10^{3})$ times for all values of $\alpha$ we need to choose
the parameters $L,N,$ and the time increment of the numerical evolution
method, $dt$, carefully. We discuss this in detail in Appendix B.2.

Figure 1:  Log-log plot of the average number density of kinks $n_{K}$ as a
function of time for different values of the power law index $\alpha$. The
power law indices are color coded as follows: $\alpha=0$ (purple),
$\alpha=0.1$ (red), $\alpha=0.2$ (darker green), $\alpha=0.3$ (orange),
$\alpha=0.4$ (blue), $\alpha=0.5$ (darker yellow), $\alpha=0.6$ (brown),
$\alpha=0.7$ (pink), $\alpha=0.8$ (darker blue), and $\alpha=0.9$ (green). The
parameters used are $t_{\rm pt}=1$, $t_{0}=0.5$, $\tau=0.01$, $a_{0}=1$ (in
units where $m=1$). For details on the choice of $L$, $N$, and $dt$, see
Appendix B.2.

In Fig. 1, we show the average number density of kinks $n_{K}(t)$ as a
function of time $t$ for different values of the index $\alpha$ associated
with the power law growth of the scale factor (1). The phase transition occurs
at $t_{\rm pt}=1$, before which $n_{K}=0$. The quench timescale is taken to be
$\tau=0.01$, which makes the phase transition quasi-instantaneous. Thus
immediately after the phase transition kinks and antikinks are quickly
produced with random positions and velocities and $n_{K}$ reaches a maximum
value almost instantaneously. Then, as kink/antikink annihilation becomes
important, $n_{K}$ begins to decrease. At late times, similarly to Refs. [38,
39], we observe a power law decline with $n_{K}(t)\propto t^{-\kappa}$. For
$\alpha=0$ (Minkowski spacetime) we reproduce the results from [38, 39] where
we see that at late times $n_{K}(t)\propto t^{-1/2}$, or $\kappa=0.5$. As the
value of $\alpha$ increases, the slope associated with the late time power law
decay of $n_{K}$ also increases. This is due to the fact that the expansion
helps dilute the leftover kinks and antikinks.

Figure 2:  Log-log plot of the average kink number density $n_{K}$ as a
function of time for different values of the quench timescale $\tau=0.01$
(solid purple), $0.05$ (dashed red), and $0.1$ (dot-dashed orange). The inset
provides a zoomed in view of the plot at early times after the phase
transition. The parameters used are $t_{\rm pt}=1$, $t_{0}=0.5$, $\alpha=0.5$,
$L=32,\ N=12800,\ dt=0.25$, $a_{0}=1$ (in units where $m=1$). The inset is
made using a finer temporal step, $dt=0.0025$ (see Appendix B.2).

In Fig. 2 we show $n_{K}$ as a function of time for different values of the
quench timescale $\tau$ ($0.01$, $0.05$ and $0.1$) and $\alpha=0.5$. We can
see that at late times the curves for $n_{K}(t)$ overlap, which proves the
universality of the power law decay (or its insensitivity to the details of
the phase transition). However the growth of $n_{K}$ immediately after the
phase transition is different for the three cases as can be seen from the
inset in Fig. 2. We note that as $\tau$ increases the growth rate of $n_{K}$
decreases. (In fact the maximum value of $n_{K}$ also decreases but this is
not clearly visible on the plot.) This implies that the quicker the phase
transition (shorter quench timescales) the more efficient the kink production
rate is.

Figure 3:  The late time power law index $n_{K}(t)$, $\kappa$, given by
$n_{K}(t)\sim t^{-\kappa}$ as a function of the power law index, $\alpha$,
associated with the growth of the scale factor (see Eq. (1)). The dashed blue
line shows the broken linear function fit to the data points. For
$\alpha\lesssim 1$ the dashed line has a slope of $\sim 1$, that is,
$\kappa\sim\alpha$. The solid green line shows $\kappa=\alpha$. For details on
the choice of $L$, $N$, and $dt$, see Appendix B.2.

To make the relationship between the expansion rate and the power law followed
by $n_{K}$ at late times more explicit, we also plot the power law index
$\kappa$ versus $\alpha$ in Fig. 3. As $\alpha$ varies from $0$ to $1$,
$\kappa$ smoothly interpolates between $0.5$ and a linear function with slope
$1$. This agrees with the analytic estimate of the previous section. For
$\alpha\lesssim 1$, $\kappa\sim\alpha$ which means that the average number of
kinks and antikinks per comoving volume remains approximately constant and
they are unable to mutually annihilate because of the Hubble flow: they are
frozen and behave like a dilute non-interacting gas. Any further dilution is
entirely due to the increase of the physical size of the lattice $a(t)L$. For
the sake of completeness we fit a broken linear function
$0.5(1+(2\alpha)^{n})^{1/n}$ to $\kappa(\alpha)$ and find a best fit value of
the parameter $n=5.55$.

Similar results are obtained for the domain wall length density in two spatial
dimensions. In particular at late times and for a given $\alpha$, the length
density of domain walls, $\mu_{DW}$, also decays following a power law
$t^{-\kappa}$. The associated annihilation process is harder to visualize
since domain walls can have complex spatial structure in two dimensions and
are no longer point-like. But it is not difficult to convince oneself that
because they are created with random shapes, positions and velocities during
the phase transition, and they percolate the lattice, at late times they are
expected to disappear, and the whole lattice is expected to end up in one
domain. We plot $\kappa$ with respect to $\alpha$ in Fig. 4, and again notice
a smooth interpolating behavior between $\kappa=0.5$ at $\alpha=0$ (Minkowski
spacetime) and $\kappa\sim\alpha$ for $\alpha\lesssim 1$. The latter case
corresponds to a frozen network of domain walls blown away by the Hubble flow
and whose length scales as $a(t)$. A broken linear function fit to
$\kappa(\alpha)$ this time yields a best fit value of the parameter $n=4.35$.
Because of our limited numerical precision, it is unclear whether the
difference in parameters $n$ between the one- and two-dimensional cases is
significant.

Figure 4:  The late time power law index for $\mu_{\rm DW}(t)$, $\kappa$,
given by $\mu_{\rm DW}(t)\sim t^{-\kappa}$, as a function of the power law
index, $\alpha$, associated with the growth of the scale factor (see Eq. (1)).
The dashed orange line shows the broken linear function fit to the data
points. For $\alpha\lesssim 1$ the dashed line has a slope of $\sim 1$, that
is, $\kappa\sim\alpha$. The solid purple line shows $\kappa=\alpha$. The error
bars and the choice of parameters $L$, $N$, and $dt$ are discussed in Appendix
B.2.

##  V VOS model

It is worth looking at our results from the perspective of the velocity-one-
scale (VOS) model – an effective description of the DW network in terms of two
quantities: $L(t)$, the correlation length/average separation between walls,
and $v(t)$ the typical (root-mean-squared) wall velocity. For standard DW
networks the VOS model was introduced in [47] (see also [48] and references
therein).

As argued in [40], VOS equations for a network of precursor DWs in an
expanding universe should reduce to

$\displaystyle\dot{L}$ $\displaystyle=$ $\displaystyle H\,L+c\,v\,,$ (30)
$\displaystyle\dot{v}$ $\displaystyle=$
$\displaystyle-k\frac{v^{2}}{L}-d\,H\,v$ (31)

where $d$ is the number of space dimensions, $H=\dot{a}/a=\alpha/t$ and
$c,\,k$ the VOS parameters. The latter are expected to depend on $\alpha$ in
general.

Equations (30)-(31) have 2 special solutions. The trivial one,

$L(t)\propto a(t)\,,\qquad v=0\,,$ (32)

corresponds to a ‘gas’ of noninteracting DWs. Wall-wall interactions are
frozen, and the walls are carried and blown away by the expansion.

The other solution, for a power-law model
$a(t)=\left(t/t_{0}\right)^{\alpha}$, is

$\begin{split}L(t)&=L_{0}\left(t/t_{0}\right)^{\kappa}\,,\\\
v(t)&=v_{0}\,\left(t/t_{0}\right)^{\kappa-1}\end{split}$ (33)

with

$\displaystyle\kappa$ $\displaystyle=$
$\displaystyle\alpha+{\scriptstyle\mathcal{V}}\,,$ (34)
$\displaystyle{\scriptstyle\mathcal{V}}$ $\displaystyle=$
$\displaystyle\frac{1-(d+1)\,\alpha}{1+k/c}~{},$ (35)

where we have introduced the ‘velocity parameter’ as the combination

${\scriptstyle\mathcal{V}}\equiv\frac{v(t)}{L(t)/c\,t}=\frac{v_{0}}{L_{0}/c\,t_{0}}~{},$
(36)

which measures the typical velocity of the DWs in units of the characteristic
velocity $L/c\,t$. (We introduce a conventional factor $c$ in (36) for
convenience). In contrast with the standard cosmological attractor of
nonlinear DWs, here $L_{0}$ is a free parameter. The scaling (33) is not
characterized by a specific number of Hubble-sized DWs per Hubble patch,
because this number increases steadily.

The parameter ${\scriptstyle\mathcal{V}}$ is more relevant. Physically it
plays the role of a ‘temperature’, the amount of DW motion. For any given
$\alpha$ and $d$, it is in one-to one correspondence with (the ratio of) the
VOS parameters $c,\,k$, see (34), so in a way it is a more physical
parameterization. Importantly, it is clear from (34) that
${\scriptstyle\mathcal{V}}$ also determines the scaling exponent $\kappa$,
that is, how $1/L$ and $v$ decay in time. This makes manifest that the non-
relativistic scaling of precursor DWs is tied to the velocity distribution (in
particular the velocity parameter ${\scriptstyle\mathcal{V}}$), as noted
before in flat space [40]. From the VOS perspective, there is little more to
be said without knowing ${\scriptstyle\mathcal{V}}$, that is, specifying the
typical velocity $v$ in addition to the typical correlation length $L$.

It would require considerable additional work to compute directly the
distribution of velocities of kinks/DWs, or even their characteristic
velocity. However, assuming the VOS equations (30), (31), we can already
extract the velocity parameter ${\scriptstyle\mathcal{V}}$ from (34) together
with the results shown in Figs. 3 and 4777Similarly, using (34) and (35) we
could read off how the ratio of VOS parameters $k/c$ depends on $\alpha$.. The
result is shown in Fig. 5. The nonrelativistic scaling ($\alpha\lesssim 1/2$)
- DW gas ($\alpha\gtrsim 1/2$) transition is manifest: $v(t)$ goes from an
$O(1)$ fraction of the characteristic velocity for $\alpha\to 0$ and drops
almost linearly in $\alpha$ until $\alpha\sim 1/2$ and then remains
suppressed.

Figure 5:  Typical velocity in function of the expansion rate parameter
$\alpha$ inferred from the VOS model. The vertical axis shows $v$ in units of
$L/(c\,t)$. See Eq. (36).

The interpretation of these results is clear. In flat space, $\alpha=0$,
${\scriptstyle\mathcal{V}}$ coincides with the scaling exponent $\kappa$ and
it takes the ‘universal’ value888Equivalently, $k/c\to 1$ for $\alpha\to 0$
${\scriptstyle\mathcal{V}}=1/2$, forced upon by non-the relativistic scaling
that emerges at late times [40]. The calculations in Sections II and III
extend this to $\alpha\neq 0$. The expansion forces the velocity parameter to
decrease. However, for $\alpha\lesssim 1/2$, the correlation length exponent
must remain close to $1/2$ because the correlation length is sub-Hubble – the
physical cosmic horizon grows faster than $L$ and so the effects of expansion
must be small. In order for $\kappa(\alpha)$ to be approximately constant (as
seen in Figs. 3 and 4), the velocity parameter ${\scriptstyle\mathcal{V}}$
must be linear in $\alpha$. This explains the dependence of $\kappa$ and
${\scriptstyle\mathcal{V}}$ on $\alpha$. At the core, for slow enough
expansion ($\alpha\lesssim 1/2$) the scaling is very close to the standard
nonrelativistic one $L\sim t^{1/2}$. For faster expansion ($\alpha\sim 1/2$),
this scaling is completely lost and the DWs are basically blown away with the
Hubble flow.

##  VI Discussion

We studied the formation of domain walls (DWs) as a result of a cosmological
$\mathbb{Z}_{2}$-symmetry breaking phase transition in an expanding
background. Our description includes full quantum field theory effects but is
limited to the regime of small fluctuations, that is to the precursor DWs
regime. We find that the domain wall area density scales universally as
$t^{-\kappa}$ at late times, independently of the details of the quench, with
$\kappa$ depending only on the growth rate of the scale factor $\alpha=t\,H$
and the number of dimensions. Just like in the flat spacetime case, the quench
time scale only impacts the domain wall area density at early times after the
phase transition. The value of $\kappa$ in that case, $\alpha=0$, is known to
be $1/2$ [38, 39], which is interpreted as the wall network entering a
nonrelativistic scaling regime [40]. In presence of expansion, $\kappa$
depends on $\alpha$ in a distinctive way, interpolating smoothly between
$\kappa\simeq 0.5$ for small growth rates $\alpha\lesssim 0.5$, and
$\kappa\sim\alpha$ for $\alpha\gtrsim 0.5$. This behaviour is understood as
follows. For $\alpha\lesssim 0.5$, the nonrelativistic scaling prevails (with
small corrections) because the correlation length is always sub-Hubble and so
one is ‘close’ to the flat space case. However, for $\alpha>1/2$ the scaling
$L\sim t^{1/2}$ cannot be sustained, because the network is dynamical so $L$
is bound to grow at least as fast as the scale factor, $t^{\alpha}$. Thus, the
expectation from Ref. [40] that there is a transition from nonrelativistic
scaling to DW gas behaviour around $\alpha\simeq 0.5$ is confirmed. Moreover,
we find the that the transition (the shape of $\kappa(\alpha)$) is smooth.

Even though, for numerical reasons, our results are limited to low spatial
dimensions, we expect qualitatively similar results to hold in the physically
relevant three-dimensional case (and even in higher $d$). Indeed while the
$1+1$ dimensional case is special because domain walls there have no internal
degrees of freedom (they are point-like kinks and antikinks), already the two-
dimensional case captures the richness of the continuously infinite number of
degrees of freedom expected of domain walls in higher dimensions.

In fact, a somewhat unexpected outcome of our work is that the details of the
(nonrelativistic scaling) - (DW gas) transition depend quite weakly on
dimensionality, as is evident from the similitude of Figs. 3 (valid for $d=1$)
and 4 (valid for $d=2$). At first glance this seems very reasonable because
DWs are codimension-1 objects, so a DW network in any number of dimensions
should retain features of a bath of particles in $1+1$. However, DWs in $d+1$
interact/interconnect in much richer ways, so $O(1)$ discrepancies could have
been present. In this sense, the most dramatic comparison is with the $d=1$
case where DWs are point-like. In other words, the closeness between the $d=1$
and $d=2$ cases suggests that $d=3$ should not differ greatly, because the
major qualitative difference – the extended character of the DWs – is already
present in $d=2$.

This suggests the following extrapolation of our results to $d=3$ as a
reasonable expectation. One would expect $\kappa\to 1/2$ for $\alpha\to 0$,
$\kappa\to\alpha$ for $\alpha\gtrsim 1/2$. Concerning the most interesting
case phenomenologically, $\alpha=1/2$ or radiation domination, the tendency
seems to be that $\kappa$ should depart from $1/2$, slightly increasing its
value between $d=1$ and $d=2$. Assuming that this tendency is preserved, this
points to a ballpark value of $\kappa\sim 0.6$. It would be interesting to see
if this is confirmed.

Another interesting extension of our work is to include temperature. Thermal
phase transitions have been considered in the instantaneous quench limit [41,
42, 43, 44]. A network of precursor DWs is expected to display the same
‘universal’ late-time scaling behaviour derived here. It would be interesting
to extend our calculation and check that this is indeed the case.

An important question which we haven’t addressed is of course how feasible it
is to have a precursor DW network epoch in realistic models. Is it possible
that a DW network would spend a noticeable amount of time in an early scaling
regime such as the one described here? This requires very weak coupling or
that the initial field amplitude at the quench be very small, which do not
seem unreasonable conditions. Another possibility is that the precursor limit
captures properties of standard, nonlinear DWs in friction-dominated models,
as these fall in a non-relativistic regime too. We leave this for future
research.

For the moment, let us point out that even if it turns out that a long period
with scaling precursor DW networks cannot be realized in phenomenologically
interesting models, precursor DWs still provide a nice and simple theoretical
laboratory to understand the behaviour of extended objects in cosmology.
Indeed, it is quite intriguing that interacting DW-like objects with rich
‘dynamics’ emerge from a free theory.

###### Acknowledgements.

We thank Tanmay Vachaspati for useful discussions and comments. Computations
for this research were performed on the Pennsylvania State University’s
Institute for Computational and Data Sciences’ Roar supercomputer. We also
acknowledge the use of the Deep Learning for Statistics, Astrophysics,
Geoscience, Engineering, Meteorology and Atmospheric Science, Physical
Sciences and Psychology (DL-SAGEMAPP) computing cluster at the Institute for
Computational and Data Sciences (ICDS) at the Pennsylvania State University.
M. M. and G. Z. wish to thank the Yukawa Institute for Theoretical Physics,
Kyoto University for hospitality during the period when this project was
initiated. M. M. is supported by National Science Foundation Grant No.
AST-2108466. M. M. also acknowledges support from the Institute for
Gravitation and the Cosmos Postdoctoral Fellowship at Pennsylvania State
University. OP acknowledges the support from the Departament de Recerca i
Universitats from Generalitat de Catalunya to the Grup de Recerca ‘Grup de
Física Teòrica UAB/IFAE’ (2021 SGR 00649) and the Spanish Ministry of Science
and Innovation (PID2020-115845GB-I00/AEI/10.13039/501100011033). IFAE is
partially funded by the CERCA program of the Generalitat de Catalunya. G. Z.
is supported by SONATA BIS grant 2023/50/E/ST2/00231 from the Polish National
Science Centre.

##  Appendix A Schrödinger picture quantization

Below we outline the steps leading to Eq. (4). In other words we derive the
Schrödinger picture quantum evolution of the wavefunctional describing the
quantum state of the field $\phi$. In general this can only be done in a
consistent way by regularizing the model on a discrete lattice.

Assuming periodic boundary conditions at spatial infinity, we start by doing
an integration by parts on (2) thus rewriting the Lagrangian density as

$\mathcal{L}=\int dx\
a(t)\left[\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2}a(t)^{-2}\phi\phi^{\prime\prime}-\frac{1}{2}m^{2}(t)\phi^{2}\right].$
(37)

As described in the main text, we then discretize it on an $N$ point periodic
lattice of size $L$ and lattice spacing $\ell=L/N$. For any given lattice
point $i$ the value of the field $\phi$ at that point is denoted $\phi_{i}$.
The periodicity of the lattice dictates, $\phi_{N+1}=\phi_{1}$ and
$\phi_{0}=\phi_{N}$. We use the central finite difference approximation to
lowest order for the spatial double derivative, that is

$\phi^{\prime\prime}_{i}\rightarrow\frac{\phi_{i-1}-2\phi_{i}+\phi_{i+1}}{\ell^{2}}.$
(38)

With these prescriptions, the discretized Lagrangian density reads

$\displaystyle\mathcal{L}_{\rm disc.}$ $\displaystyle=$ $\displaystyle
a(t)\ell\sum_{i=1}^{N}\left[\frac{1}{2}\dot{\phi}_{i}^{2}-\frac{1}{2(a(t)\ell)^{2}}\phi_{i}\right.$
(39)
$\displaystyle\times\big{(}\phi_{i-1}-\left.2\phi_{i}+\phi_{i+1}\big{)}-\frac{1}{2}m^{2}(t)\phi_{i}^{2}\right].$

In matrix notation the above equation can be rewritten more compactly as

$\mathcal{L}_{\rm
disc.}=\frac{a(t)\ell}{2}\dot{\bm{\phi}}^{T}\dot{\bm{\phi}}-\frac{a(t)\ell}{2}\bm{\phi}^{T}\bm{\Omega}^{2}(t)\bm{\phi},$
(40)

where, $\bm{\phi}\equiv\big{(}\phi_{1},\phi_{2},\dots,\phi_{N}\big{)}^{T}$,
and $\bm{\Omega}^{2}(t)$ is given by

$[\bm{\Omega}^{2}(t)]_{ij}=\begin{cases}+{2}/(a(t)\ell)^{2}+m^{2}(t),&i=j\\\
-{1}/(a(t)\ell)^{2},&i=j\pm 1\ (\text{mod}\ N)\\\
0,&\text{otherwise}.\end{cases}$ (41)

The full field theory has thus effectively been reduced to an $N$ particle
quantum mechanics model whose quantum state (in the Schrödinger picture) is
fully determined by a wavefunctional $\Psi(t;\bm{\phi})$ obeying the
Schrödinger equation

$i\frac{\partial\Psi}{\partial
t}=-\frac{1}{2a(t)\ell}\sum_{i=1}^{N}\frac{\partial^{2}\Psi}{\partial\phi_{i}^{2}}+\frac{a(t)\ell}{2}\bm{\phi}^{T}\bm{\Omega}^{2}(t)\bm{\phi}\
\Psi\,.$ (42)

The last remaining ingredient is specifying a reasonable initial state. We
choose the field to be in the instantaneous ground state of the $N$ particle
quantum mechanics model at time $t_{0}$: this choice corresponds to the zeroth
order adiabatic vacuum in a field theoretic language.

Now, since our model is quadratic, the state will remain Gaussian as time
evolves. Following [38, 39, 40], we then write

$\Psi(t;\bm{\phi})=\mathcal{N}(t)\ {\rm
exp}\left[\frac{ia(t)\ell}{2}\bm{\phi}^{T}\bm{M}(t)\bm{\phi}\right],$ (43)

where $\bm{M}(t)$ is a $N\times N$ complex symmetric matrix verifying
$M(t_{0})=i\bm{\Omega}^{2}(t_{0})^{1/2}$ and $\mathcal{N}$ is a normalization
factor. Plugging this ansatz in the Schrödinger equation (42) we obtain the
matrix differential equation

$\dot{\bm{M}}(t)+\bm{M}^{2}(t)+\bm{\Omega}^{2}(t)+H(t)\bm{M}(t)=0.$ (44)

With the help of the mode functions $c_{n}(t)$ introduced in Eqs. (6), (7) and
(8) the matrix elements of $\bm{M}(t)$ can then be written as,

$[\bm{M}(t)]_{jk}=\frac{1}{N}\sum_{n=0}^{N-1}c_{n}(t)^{-1}\dot{c}_{n}(t)e^{2i\pi
n(j-k)/N}.$ (45)

Here the fact that $c_{n}=c_{N-n}$ for $1\leq n\leq N-1$ can be used to
convince oneself that $\bm{M}(t)$ is indeed symmetric. Using the Wronskian
constraint

$c_{n}(t)^{*}\dot{c}_{n}(t)-c_{n}(t)\dot{c}_{n}(t)^{*}=\frac{i}{a(t)\ell},$
(46)

which can be inferred from (6), (7) and (8), the probability density
functional $\mathcal{P}(t;\bm{\phi})\equiv|\Psi(t;\bm{\phi})|^{2}$ can finally
be expressed in terms of the covariance matrix (5) as in Eq. (4).

##  Appendix B Numerics

In this appendix we give some details about the numerical equations and
methods used. We also discuss the choice of parameters for our numerical
integration method and its sensitivity thereon.

###  B.1 Equivalent evolution equations

This section focuses on the derivation of the numerical equations we
effectively use to compute the number density of kinks. Analogous equations
can be obtained for domain walls.

The idea is to recast the differential equations (6) in a form that is more
amenable to numerical treatment. Indeed, the exponential growth of the
$c_{n}$s corresponding to unstable modes quickly leads to very large numbers
that exceed machine precision thus limiting the duration for which the
numerical simulations are reasonably accurate. This problem can be
circumvented by factoring out the leading exponential behavior (see [38]).
More precisely, by writing $c_{n}(t)\equiv\rho_{n}(t)e^{i\theta_{n}(t)}$ and
using (46) we find

$\rho_{n}^{2}\dot{\theta}_{n}=\frac{1}{2a(t)\ell},$ (47)

which allows us to rewrite (6) as

$\displaystyle\ddot{\rho}_{n}+H(t)\dot{\rho}_{n}+\left[\frac{4}{(a(t)\ell)^{2}}{\rm
sin}^{2}\left(\frac{\pi n}{N}\right)+m^{2}(t)\right]\rho_{n}$
$\displaystyle\hskip 85.35826pt=\frac{1}{4(a(t)\ell)^{2}\rho_{n}^{3}}.$ (48)

The corresponding initial conditions are

$\displaystyle\rho_{n}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2a_{0}\ell}}\left[\frac{4}{(a_{0}\ell)^{2}}\sin^{2}\left(\frac{\pi
n}{N}\right)+m^{2}(t_{0})\right]^{-1/4}$ (49)
$\displaystyle\dot{\rho}_{n}(t_{0})$ $\displaystyle=$ $\displaystyle 0.$ (50)

Noticing that, for $t>t_{\rm pt}$, the fastest growing unstable mode is the
zero mode and that $\rho_{0}(t)\sim e^{mt}$ for large $t$, we are led to
define $q(t)\equiv\ln(\rho_{0}(t))$ and
$r_{n}(t)\equiv\rho_{n}(t)e^{-q_{n}(t)}$. From (B.1), (49) and (50) we then
obtain

$\ddot{q}+\dot{q}^{2}+H(t)\dot{q}+m^{2}(t)=\frac{e^{-4q}}{4(a(t)\ell)^{2}}\,,$
(51)

with initial conditions

$\displaystyle q(t_{0})$ $\displaystyle=$
$\displaystyle\ln\left[\frac{1}{\sqrt{2a_{0}\ell}}\Big{(}m^{2}(t_{0})\Big{)}^{-1/4}\right],$
(52) $\displaystyle\dot{q}(t_{0})$ $\displaystyle=$ $\displaystyle 0,$ (53)

and

$\displaystyle\ddot{r}_{n}+\left(2\dot{q}+H(t)\right)\dot{r}_{n}+$ (54)
$\displaystyle+\frac{4}{(a(t)\ell)^{2}}\left[\sin^{2}\left(\frac{\pi
n}{N}\right)+\frac{e^{-4q}}{16}\left(1-r_{n}^{-4}\right)\right]r_{n}=0,$ (55)

with initial conditions

$\displaystyle r_{n}(t_{0})$ $\displaystyle=$
$\displaystyle\left[1+\frac{4}{(a_{0}\ell)^{2}m^{2}(t_{0})}\sin^{2}\left(\frac{\pi
n}{N}\right)\right]^{-1/4}$ (56) $\displaystyle\dot{r}_{n}(t_{0})$
$\displaystyle=$ $\displaystyle 0.$ (57)

By construction, Eqs. (51) and (54) do not suffer from the precision issues
that plague the mode function evolution equations in their original form (6)
and allow for long integration times.

###  B.2 Miscellanea

(a) (a)

(b) (c)

(c) (b)

(d) (d)

Figure 6:  Number density of kinks ($n_{K}$, left) and length density of
domain walls ($\mu_{DW}$, right) as a function of time, for different choices
of the numerical evolution parameters $L$, $N$ (for $\alpha=0$ or $0.5$), and
$dt$ (for $\alpha=0.8$ or $0.9$). The _upper panels_ illustrate the UV and IR
cutoff insensitivity of the results, while the _lower panels_ illustrate their
independence with respect to a finer choice of the temporal time step.

We now mention several subtle points that have an impact on how we carry out
our numerical computations.

The fraction $u(t)$ of unstable modes in the system depends monotonically on
$\ell t^{\alpha}$ as can be inferred from (14). To maintain accuracy, one
needs to make sure that $u$ doesn’t reach $\mathcal{O}(1)$ before the end of
the numerical evolution. (In other words we need the unstable modes to not
completely deplete the entire set of modes.) For a given total duration of the
numerical evolution, this requires choosing a sufficiently small $\ell$.
Moreover, the larger $\alpha$ is, the faster the expansion, and the smaller
$\ell$ needs to be. This requirement is compatible with our focus on the
continuum limit of the discretized system. Since we are ultimately interested
in the infinite volume limit we also require $L$ (and consequently $N=L/\ell$)
to be large.

An additional difficulty is related to the time increment $dt$ in our Crank-
Nicholson numerical method. It turns out that the number of modes summed over
in (16) and (17) becomes quite large for $\alpha>0.3$ and the numerical error
in the computation of $c_{n}$ (depending on $dt$ and typically quite small)
adds up to reach unacceptably large values if the time increment is not small
enough. For faster expansion rates we need to substantially decrease $dt$ and
therefore, in order to keep the runtime reasonably short, we are forced to
decrease $N$ as well.999Since $u$ and $\ell$ are fixed by the desired
evolution duration, $L$ needs to be smaller as well.

For kinks in $1+1$ dimensions the computational complexity of the numerical
method is $\mathcal{O}(N)$ and it is easy to satisfy the above conditions. We
choose the lattice spacing to be $\ell=0.0025$, which allows numerical
evolution durations of order $10^{3}m^{-1}$ for even the largest $\alpha$s.
This is achieved by choosing $L=320,N=128000$ for
$\alpha\in\left[0.0,0.3\right]$ and $L=32,N=12800$ for
$\alpha\in\left(0.3,0.9\right]$. We choose $dt=0.5$ for
$\alpha\in\left[0.0,0.3\right]$, $dt=0.25$ for
$\alpha\in\left(0.3,0.7\right]$, and $dt=0.0625$ for
$\alpha\in\left(0.7,0.9\right]$. We see that as $\alpha$ increases we decrease
$dt$, as well as $N$ and $L$, to keep numerical error under control.

The above choices allow an accurate determination of the late time power law
followed by the kink number density as a function of $\alpha$, as depicted in
Fig. 3. However it turns out they are not sufficient for accurately capturing
the immediate aftermath of the phase transition (where boundary effects are
non-negligible for the previous choices of $L$). To achieve the desired
precision (without having to increase $N$ as well) for very early times we
perform shorter, lower-resolution runs (up to $t\sim 50m^{-1}$) with
$L=3200,N=12800,$ and $dt=0.25$ (for all values of $\alpha$). We then combine
the shorter runs with the main runs to obtain the smooth curves shown in Figs.
1 and 2.

For domain walls in $2+1$ dimensions, numerical limitations place even greater
constraints on our choices of parameters for the numerical evolution. Indeed
the computational complexity is $\mathcal{O}(N^{2})$ and we are effectively
limited $N\lesssim\mathcal{O}(10^{3})$. We choose $L=1000$ and $N=2000$ for
the static background ($\alpha=0$). In this case we obtain $\kappa\approx
0.49$, a slightly different value from the exact one, $0.5$. This means our
choice of $L$ and $N$ leads to higher uncertainty results for the values of
$\kappa$, but it also provides a reasonable estimate for the expected error,
$\sim 2\%$. This is encoded in the error bars around the points of Fig. 4. We
also need to contend with the same issues encountered in the case of kinks
(and explained in the previous paragraphs). To ensure a numerical evolution
duration of order $10^{3}m^{-1}$, we fix $N=2000$ and decrease $L$ with
increasing $\alpha$ as follows: $L=1000$ for $\alpha=0.1$, $L=800$ for
$\alpha=0.2$, $L=300$ for $\alpha=0.3$, $L=100$ for $\alpha=0.4$, $L=60$ for
$\alpha=0.5$, $L=30$ for $\alpha=0.6$, $L=10$ for $\alpha=0.7$, $L=6$ for
$\alpha=0.8$, and $L=2$ for $\alpha=0.9$. The choice of such small $L$s is
justified because the expansion increases the physical size of the box and we
expect boundary effects to be negligible in the late time limit where the
domain wall network is sufficiently sparse. To maintain nunmerical precision
we also need to decrease the time increment with increasing $\alpha$ as
follows: $dt=0.5$ for $\alpha\in\left[0.0,0.7\right)$, $dt=0.25$ for
$\alpha\in\left[0.7,0.8\right)$, and $dt=0.125$ for
$\alpha\in\left[0.8,0.9\right]$.

Finally, as a sanity check of the above parameter choices, in Fig. 6, we
demonstrate the UV- and IR-independence of our results for two different
values of $\alpha$ ($0$ and $0.5$). We also show that our choice of time
evolution increment is sufficiently fine for the values of $\alpha$ most prone
to numerical instabilities ($0.8$ and $0.9$). This proves the robustness of
our results despite the computational limitations.

## References

  * Agazie  _et al._ [2023] G. Agazie _et al._ (NANOGrav), The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background, Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023a] J. Antoniadis _et al._ (EPTA), The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals, Astron. Astrophys. 678, A50 (2023a), arXiv:2306.16214 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023b] J. Antoniadis _et al._ (EPTA), The second data release from the European Pulsar Timing Array I. The dataset and timing analysis, Astron. Astrophys. 678, A48 (2023b), arXiv:2306.16224 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023c] J. Antoniadis _et al._ (EPTA), The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe, (2023c), arXiv:2306.16227 [astro-ph.CO] .
  * Reardon  _et al._ [2023] D. J. Reardon _et al._ , Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array, Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Zic  _et al._ [2023] A. Zic _et al._ , The Parkes Pulsar Timing Array Third Data Release, (2023), arXiv:2306.16230 [astro-ph.HE] .
  * Xu  _et al._ [2023] H. Xu _et al._ , Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I, Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Bian  _et al._ [2021] L. Bian, R.-G. Cai, J. Liu, X.-Y. Yang, and R. Zhou, Evidence for different gravitational-wave sources in the NANOGrav dataset, Phys. Rev. D 103, L081301 (2021), arXiv:2009.13893 [astro-ph.CO] .
  * Sakharov  _et al._ [2021] A. S. Sakharov, Y. N. Eroshenko, and S. G. Rubin, Looking at the NANOGrav signal through the anthropic window of axionlike particles, Phys. Rev. D 104, 043005 (2021), arXiv:2104.08750 [hep-ph] .
  * Bian  _et al._ [2022] L. Bian, J. Shu, B. Wang, Q. Yuan, and J. Zong, Searching for cosmic string induced stochastic gravitational wave background with the Parkes Pulsar Timing Array, Phys. Rev. D 106, L101301 (2022), arXiv:2205.07293 [hep-ph] .
  * Ferreira  _et al._ [2023] R. Z. Ferreira, A. Notari, O. Pujolas, and F. Rompineve, Gravitational waves from domain walls in Pulsar Timing Array datasets, JCAP 02, 001, arXiv:2204.04228 [astro-ph.CO] .
  * Afzal  _et al._ [2023] A. Afzal _et al._ (NANOGrav), The NANOGrav 15 yr Data Set: Search for Signals from New Physics, Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Blasi  _et al._ [2023] S. Blasi, A. Mariotti, A. Rase, and A. Sevrin, Axionic domain walls at Pulsar Timing Arrays: QCD bias and particle friction, JHEP 11, 169, arXiv:2306.17830 [hep-ph] .
  * Ellis  _et al._ [2024] J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, What is the source of the PTA GW signal?, Phys. Rev. D 109, 023522 (2024), arXiv:2308.08546 [astro-ph.CO] .
  * Ferreira  _et al._ [2024] R. Z. Ferreira, A. Notari, O. Pujolàs, and F. Rompineve, Collapsing Domain Wall Networks: Impact on Pulsar Timing Arrays and Primordial Black Holes, (2024), arXiv:2401.14331 [astro-ph.CO] .
  * Kibble [1976] T. W. B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9, 1387 (1976).
  * Vilenkin and Shellard [2000] A. Vilenkin and E. P. S. Shellard,  _Cosmic Strings and Other Topological Defects_ (Cambridge University Press, 2000).
  * Vachaspati [2007] T. Vachaspati,  _Kinks and Domain Walls : An Introduction to Classical and Quantum Solitons_ (Oxford University Press, 2007).
  * Saikawa [2017] K. Saikawa, A review of gravitational waves from cosmic domain walls, Universe 3, 40 (2017), arXiv:1703.02576 [hep-ph] .
  * Sikivie [1982] P. Sikivie, Of Axions, Domain Walls and the Early Universe, Phys. Rev. Lett. 48, 1156 (1982).
  * Vilenkin and Everett [1982] A. Vilenkin and A. E. Everett, Cosmic Strings and Domain Walls in Models with Goldstone and PseudoGoldstone Bosons, Phys. Rev. Lett. 48, 1867 (1982).
  * Zeldovich  _et al._ [1974] Y. B. Zeldovich, I. Y. Kobzarev, and L. B. Okun, Cosmological Consequences of the Spontaneous Breakdown of Discrete Symmetry, Zh. Eksp. Teor. Fiz. 67, 3 (1974).
  * Press  _et al._ [1989] W. H. Press, B. S. Ryden, and D. N. Spergel, Dynamical Evolution of Domain Walls in an Expanding Universe, Astrophys. J. 347, 590 (1989).
  * Larsson  _et al._ [1997] S. E. Larsson, S. Sarkar, and P. L. White, Evading the cosmological domain wall problem, Phys. Rev. D 55, 5129 (1997), arXiv:hep-ph/9608319 .
  * Garagounis and Hindmarsh [2003] T. Garagounis and M. Hindmarsh, Scaling in numerical simulations of domain walls, Phys. Rev. D 68, 103506 (2003), arXiv:hep-ph/0212359 .
  * Oliveira  _et al._ [2005] J. C. R. E. Oliveira, C. J. A. P. Martins, and P. P. Avelino, The Cosmological evolution of domain wall networks, Phys. Rev. D 71, 083509 (2005), arXiv:hep-ph/0410356 .
  * Hiramatsu  _et al._ [2010] T. Hiramatsu, M. Kawasaki, and K. Saikawa, Gravitational Waves from Collapsing Domain Walls, JCAP 05, 032, arXiv:1002.1555 [astro-ph.CO] .
  * Leite and Martins [2011] A. M. M. Leite and C. J. A. P. Martins, Scaling Properties of Domain Wall Networks, Phys. Rev. D 84, 103523 (2011), arXiv:1110.3486 [hep-ph] .
  * Kawasaki and Saikawa [2011] M. Kawasaki and K. Saikawa, Study of gravitational radiation from cosmic domain walls, JCAP 09, 008, arXiv:1102.5628 [astro-ph.CO] .
  * Hiramatsu  _et al._ [2014] T. Hiramatsu, M. Kawasaki, and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, JCAP 02, 031, arXiv:1309.5001 [astro-ph.CO] .
  * Correia  _et al._ [2014] J. R. C. C. C. Correia, I. S. C. R. Leite, and C. J. A. P. Martins, Effects of Biases in Domain Wall Network Evolution, Phys. Rev. D 90, 023521 (2014), arXiv:1407.3905 [hep-ph] .
  * Correia  _et al._ [2018] J. R. C. C. C. Correia, I. S. C. R. Leite, and C. J. A. P. Martins, Effects of biases in domain wall network evolution. II. Quantitative analysis, Phys. Rev. D 97, 083521 (2018), arXiv:1804.10761 [astro-ph.CO] .
  * Zambujal Ferreira  _et al._ [2022] R. Zambujal Ferreira, A. Notari, O. Pujolàs, and F. Rompineve, High Quality QCD Axion at Gravitational Wave Observatories, Phys. Rev. Lett. 128, 141101 (2022), arXiv:2107.07542 [hep-ph] .
  * Kitajima  _et al._ [2024] N. Kitajima, J. Lee, K. Murai, F. Takahashi, and W. Yin, Gravitational waves from domain wall collapse, and application to nanohertz signals with QCD-coupled axions, Phys. Lett. B 851, 138586 (2024), arXiv:2306.17146 [hep-ph] .
  * Chang and Cui [2023] C.-F. Chang and Y. Cui, Dynamics of Long-lived Axion Domain Walls and Its Cosmological Implications, (2023), arXiv:2309.15920 [hep-ph] .
  * Hindmarsh [1996] M. Hindmarsh, Analytic scaling solutions for cosmic domain walls, Phys. Rev. Lett. 77, 4495 (1996), arXiv:hep-ph/9605332 .
  * Hindmarsh [2003] M. Hindmarsh, Level set method for the evolution of defect and brane networks, Phys. Rev. D 68, 043510 (2003), arXiv:hep-ph/0207267 .
  * Mukhopadhyay  _et al._ [2020a] M. Mukhopadhyay, T. Vachaspati, and G. Zahariade, Quantum Formation of Topological Defects, Phys. Rev. D 102, 116002 (2020a), arXiv:2009.11480 [hep-th] .
  * Mukhopadhyay  _et al._ [2020b] M. Mukhopadhyay, T. Vachaspati, and G. Zahariade, Emergence of classical structures from the quantum vacuum, Phys. Rev. D 102, 056021 (2020b), arXiv:2004.07249 [hep-th] .
  * Pujolas and Zahariade [2023] O. Pujolas and G. Zahariade, Domain wall annihilation: A QFT perspective, Phys. Rev. D 107, 123527 (2023), arXiv:2212.11204 [hep-th] .
  * Calzetta [1989] E. Calzetta, Spinodal decomposition in quantum field theory, Annals of Physics 190, 32 (1989).
  * Boyanovsky [1993] D. Boyanovsky, Quantum spinodal decomposition, Phys. Rev. E 48, 767 (1993), arXiv:hep-th/9301095 .
  * Ibaceta and Calzetta [1999] D. Ibaceta and E. Calzetta, Counting defects in an instantaneous quench, Phys. Rev. E 60, 2999 (1999), arXiv:hep-ph/9810301 .
  * Boyanovsky  _et al._ [2000] D. Boyanovsky, H. J. de Vega, and R. Holman, Nonequilibrium phase transitions in condensed matter and cosmology: Spinodal decomposition, condensates and defects, NATO Sci. Ser. C 549, 139 (2000), arXiv:hep-ph/9903534 .
  * Watson [1944] G. Watson,  _A Treatise on the Theory of Bessel Functions_ (Cambridge University Press, 1944).
  * Olver  _et al._ [2010] F. Olver, D. Lozier, R. Boisvert, and C. Clark,  _The NIST Handbook of Mathematical Functions_ (Cambridge University Press, New York, NY, 2010).
  * Avelino  _et al._ [2005] P. P. Avelino, C. J. A. P. Martins, and J. C. R. E. Oliveira, One-scale model for domain wall network evolution, Phys. Rev. D 72, 083506 (2005), arXiv:hep-ph/0507272 .
  * Avelino  _et al._ [2022] P. Avelino, D. Grüber, and L. Sousa, Analytical Scaling Solutions for the Evolution of Cosmic Domain Walls in a Parameter-Free Velocity-Dependent One-Scale Model, Symmetry 14, 1799 (2022), arXiv:2203.16173 [astro-ph.CO] .