# Domain wall annihilation – a QFT perspective

Oriol Pujolàs1 George Zahariade1 1Institut de Física d’Altes Energies (IFAE)
and The Barcelona Institute of Science and Technology (BIST)
Campus UAB 08193 Bellaterra (Barcelona), Spain.

###### Abstract

Domain wall networks in the early universe, formed upon spontaneous breaking
of a discrete symmetry, have a rich impact on cosmology. Yet, they remain
somewhat unexplored. We introduce a new analytic strategy to understand better
the domain wall epoch, from formation to annihilation. Our method includes a
quantum field theoretical treatment of the initial state at domain wall
formation, as well as of the time evolution. We find that the domain wall area
density for a network with biased initial condition in $d+1$ dimensional flat
spacetime evolves as $t^{-1/2}\,\exp\big{(}-(t/t_{ann})^{d/2}\big{)}$. We
comment on the relation between this and previous results obtained in
condensed matter and in cosmology. The extrapolation of this law to an
expanding universe applies to networks that are close to the domain wall ‘gas’
limit.

##  I Introduction

Domain wall (DW) networks are interesting probes of physics beyond the
Standard Model (SM). DWs form whenever a discrete symmetry is spontaneously
broken. A variety of SM extensions feature DWs, among others, axion models. In
the early universe, any phase transition with spontaneous breaking of a
discrete symmetry results in a cosmological DW network, which has a strong
impact on cosmology, see [1, 2] for reviews. DW networks typically attain a
scaling or self-similar regime, where at any time $t$ the network is composed
of $O(1)$ Hubble-sized walls, and so the average network density is
$\rho_{{}_{\rm DW}}\sim\sigma_{{}_{\rm DW}}/t$, with $\sigma_{{}_{\rm DW}}$
the DW tension. Thus DWs dilute more slowly than matter, bringing about the
so-called ‘DW problem’. In practice, the DW problem is avoided if i)
$\sigma_{{}_{\rm DW}}$ is tiny ($\sigma_{{}_{\rm DW}}\lesssim MeV^{3}$, [3,
4]), or ii) the network annihilates at some time, $t_{ann}$.

The interest in DW networks is strengthened because it is possible nowadays to
search for the gravitational wave (GW) signal that they generate in GW
observatories. The signal amplitude scales like $\sigma_{{}_{\rm DW}}^{2}$ so
the loudest signals necessarily correspond to networks that annihilate at some
point. Moreover, DWs also arise in (post-inflationary) axion models, where a
hybrid DW-string network is formed below the QCD epoch and its annihilation
affects the final axion dark matter abundance. This renders the annihilating
networks of particular interest, and considerable work has been dedicated to
them recently [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23].

Various mechanisms are known to result in the annihilation of DWs. To
illustrate the various options, let us consider the discrete symmetry to be
$\mathbb{Z}_{2}$ acting on a real scalar field as $\phi\to-\phi$. For
concreteness we can assume

$V(\phi)=\frac{\lambda}{4}\phi^{4}+\frac{m^{2}(t)}{2}\phi^{2}~{}.$ (1)

The potential $V(\phi)$ is $\mathbb{Z}_{2}$ symmetric and at some time
$m^{2}(t)$ becomes negative and $V(\phi)$ develops a double-well shape with
degenerate minima. There are three basic types of annihilation mechanisms
(arising naturally in different models):

  * •

symmetry restoration: the potential returns to a single $\mathbb{Z}_{2}$
preserving vacuum and consequently the DWs disappear. This case does not
require any further breaking of $\mathbb{Z}_{2}$. However, having a DW finite
epoch would require that the mass term becomes first tachyonic and later on
positive. We won’t discuss this option in this work.

  * •

pressure bias: consisting in the addition of small explicit
$\mathbb{Z}_{2}$-breaking terms in the potential that uplift the vacuum
degeneracy. The simplest option corresponds to adding to the potential the
lowest dimension $\mathbb{Z}_{2}$-breaking operator, $\phi$. The potential
difference in the two minima, $\Delta V$, results in a pressure pushing the
walls towards the false vacuum regions.

  * •

population bias: more than including explicit symmetry breaking terms in the
Lagrangian, the symmetry is broken by assuming an uneven initial distribution
of the degenerate vacua. This can be realized for instance by having the field
displaced from $\phi=0$ at the time of the DW forming transition.

The main focus of this work is network annihilation for biased networks,
paying particular attention to how the annihilation proceeds as a function of
time. One expects that the annihilation translates into an exponential
suppression of the DW average density compared to the scaling regime. In other
words,

$\rho_{{}_{\rm DW}}=\frac{\sigma_{{}_{\rm DW}}}{t}\,S(t)$ (2)

with $S(t)$ an exponentially decaying function over some typical annihilation
scale $t_{ann}$.

The annihilation of simple $\mathbb{Z}_{2}$ DW networks has been discussed
before [24, 25, 26, 27, 28, 29, 30], both analytically and numerically (using
discretized field theory simulations). Early numerical simulations of models
with population bias in an expanding universe were performed in [25]. They
were fitted to a phenomenological form

$S_{\rm CLO}=\exp\left(-\eta/\eta_{ann}\right)$ (3)

with $\eta=\int dt/a(t)$ the conformal time. The simulation results fitted
reasonably well [25].

Exploiting the methods developed previously within solid state physics [24],
Hindmarsh obtained quite a remarkable result: the DW network energy density in
an expanding $d+1$ dimensional universe should be suppressed by [26]111The
suppressions (3) or (4) are reminiscent of percolation theory [31, 32, 25, 1],
where the number of DWs of different sizes formed at the symmetry breaking
transition was found to scale with nontrivial powers of the size in the
exponential factor. Let us just emphasize that Eqs. (3), (4) refer to the
network evolution in time, sufficiently late after the transition.

$S_{\rm
H}=\exp\left[-\;U^{2}\left(\frac{\eta}{\eta_{i}}\right)^{d}\,\right]~{}.$ (4)

Here $\eta_{i}$ is some initial time, when the network is still self-similar,
and $U$ is a dimensionless measure of the bias, as a deviation from $Z_{2}$
symmetry in the (Gaussian) distribution of the field. Note that the quadratic
dependence in $U$ can be anticipated, since annihilation must take place
independently of its sign. Writing the exponent as $-\eta^{d}/\eta_{ann}^{d}$,
one identifies that the annihilation time scales with bias as
$\eta_{ann}=\eta_{i}\,U^{-2/d}$.

As it stands, Eq. (4) is intended to apply for population bias. (The method
does not include any pressure term in the equation of motion for the walls.)
However, the decay (4) has been checked against numerical simulations and it
turns out to work well for pressure bias and not so much for population bias.
Indeed, various numerical simulations of both pressure and population biases
were performed in [27] and [28, 29], confirming that pressure bias follows (4)
but population bias instead decays less rapidly.

In an attempt to better understand the annihilation process, in the present
work we aim at solving for the network evolution and annihilation analytically
‘from scratch’ in the quantum field theory (QFT). For the moment we will
restrict to the $\mathbb{Z}_{2}$ model and to flat space.

Let us open a parenthesis to emphasize that there is a physical process
sensitive to the precise form of the exponential suppression $S(t)$. This is
the formation of primordial black holes (PBHs) from the collapse of the DW
network [33]. (See also [34, 35] for a similar PBH formation mechanism from
DWs). In essence, the network annihilation can be pictured as the
reorganization of the DW shapes so that they form closed structures similar to
false vacuum pockets. Once a closed DW fits into a Hubble volume, it collapses
under the effect of both the tension $\sigma_{{}_{\rm DW}}$ and the vacuum
energy difference $\Delta V$. If the closed DWs shrink to small enough size,
they can form black holes (BHs). The criterion for formation is that the
collapsing DW fits within the Schwarzschild radius associated to it. This can
be estimated from the total mass stored by the wall initially, when it is
Hubble-sized and roughly at rest, i.e., $R_{S}\sim(\sigma_{{}_{\rm
DW}}t^{2}+\Delta V\,t^{3})/M_{P}^{2}$. The closer this is to the Hubble scale,
$t$, the more likely it will be to form PBHs. BH formation then strongly
depends on the ‘figure of merit’ [33] defined as the ratio of the
Schwarzschild to Hubble scales, $p=R_{S}/t$. It’s easy to see that the
condition $p\gtrsim 1$ actually coincides with DW domination. This must be
avoided for a viable cosmology. The largest value $p$ reaches is at $t_{ann}$,
so we must have $p_{ann}\ll 1$. At $t_{ann}=\sigma_{{}_{\rm DW}}/\Delta V$
both bulk and surface terms contribute equally to the total DW mass so one
has, in radiation domination, $p_{ann}\sim\Delta V/T^{4}_{ann}$. Avoidance of
DW domination, $p_{ann}\ll 1$, is certainly feasible but seems an
insurmountable obstacle to produce PBHs. However, this is not the end of the
story. The reason is that the DW network contains closed walls of various
sizes. The bigger DWs in the network from the tail of the distribution (called
‘late birds’ in [33]), only fit into a Hubble volume later on, at temperatures
$T_{*}<T_{ann}$. The figure of merit for them becomes very quickly enhanced,
$p_{*}=p_{ann}(T_{ann}/T_{*})^{4}$, possibly above the collapse threshold
$p_{*}\gtrsim 1$ (without running into DW domination because this only happens
in a small fraction of the Hubble volumes). So, sufficiently late birds can
form PBHs, and their fraction is encoded in the exponential suppression in
$\rho_{{}_{\rm DW}}$. The details in the exponent are, then, important.

Returning to our computation, our method to solve for the network annihilation
is simply the extension of [36, 37] to include the bias. Ref. [36, 37] discuss
the formation of walls in $1+1$ dimensions at the $\mathbb{Z}_{2}$ symmetry
breaking transition (when $m^{2}(t)<0$) and their subsequent behaviour. At the
transition, a collection of kinks and antikinks are produced, which later on
collide and annihilate. The result, in the absence of bias and in flat space,
is that the kink-antikink network (or ‘plasma’), quickly sets into a self-
similar evolution, in which the kink/antikink number density scales as [36,
37]

$n_{\rm K}\sim\frac{1}{\sqrt{t}}~{}.$ (5)

Notice that this is a diffusive behaviour: the typical kink-antikink
separation $L\equiv 1/n_{\rm K}$ scales like $t^{1/2}$, we return to this
below.

This result is obtained from a full QFT treatment, which relies mostly on the
following observation [36, 37]: across the $\mathbb{Z}_{2}$ breaking
transition, the kinks’ positions and their evolution are already well defined
at the free theory, $\lambda=0$, level. Walls can be defined simply as places
where $\phi$ changes sign. It is then possible to track the motion of walls by
keeping only the quadratic terms in the Lagrangian throughout the
$\mathbb{Z}_{2}$ breaking transition where the mass term changes sign.

To make clear that we will be working in this free field approximation, we
introduce the terminology precursor kinks or precursor domain walls understood
as the positions/surfaces of the zeros of the free fluctuating (tachyonic)
field. Precursor walls are related to the standard DWs, but differ in
important aspects too. In the discrete symmetry breaking transition, first
precursor walls are formed. Then, as the field grows to nonlinear values they
eventually become standard domain walls. Once formed, standard DWs obey the
Nambu Goto (NG) equation to a good approximation. Precursor walls, instead, do
not obey the NG equation.

Our results are relevant for cosmological models accommodating a long
precursor wall network epoch. (A wall precursor regime can be realized
physically, for instance, if the scalar mass scale near the symmetry breaking
transition is much lighter than the Hubble rate $H$ at the symmetry breaking
transition.) The full computation in the case of an expanding universe is
beyond the scope of the present work, but we shall make some comments.

Working at the level of precursor walls will serve as an illustrative exercise
to understand the details of the network annihilation. As we will see, already
in flat space precursor wall networks follow a self-similar scaling regime –
of diffusive type similar to 5. An advantage of this QFT formalism (that can
be exported to the cosmological case) is that the initial state for the
network is properly described as the quantum vacuum of the field.

In the present work we shall do a first step in this analysis by analyzing the
statistics of annihilating precursor wall networks in flat space in presence
of a population bias. In other words, our goal is to extend the result of [36,
37] to a model with bias, and to $d+1$ dimensions. We add the bias as a
$\mathbb{Z}_{2}$ breaking term in the Lagrangian. The simplest (and most
relevant) operator is linear in the field $\phi$ with a possibly time
dependent coefficient. With no loss of generality, the resulting quadratic
potential with linear bias is

$V(\phi)\simeq\frac{1}{2}m^{2}(t)(\phi-\delta\phi(t))^{2}$ (6)

and the c-number $\delta\phi(t)$ encodes the bias directly as the position of
the maximum (for $m^{2}<0$). The formalism allows to consider a separate time
dependence for $\delta\phi(t)$ and $m^{2}(t)$, so we will keep it unspecified
as much as possible.

The simplest example of a constant tilt in the potential, $\mu^{3}\phi$,
translates to $\delta\phi(t)=-\mu^{3}/m^{2}(t)$. Thus, the instantaneous
minimum/maximum of the potential moves to $-\infty$ and reemerges from
$+\infty$, for a transient non adiabatic regime. For $m^{2}(t)$ approaching a
negative constant value, the evolution goes back to adiabatic at some point
and to a good approximation (if the quench is fast enough) the initial
wavefunctional for the field is off-center with respect to the maximum in the
potential. In other words, the sudden quench approximation in the presence of
a bias reduces to initiating the evolution as illustrated in Fig. 1.

We postpone a more general study for future work, but here we limit ourselves
to this case. As the picture suggests, we are going to capture a bias of the
population type. Indeed, since in our description the stabilization terms do
not appear, the pressure difference $\Delta V$ will not enter at any stage in
the computation. In other words, in this paper we will only be interested in
the so-called spinodal decomposition or spinodal instability phase, where the
homogeneous $\phi=0$ vacuum phase splits into different domains with $\phi<0$
and $\phi>0$ separated by domain walls (kinks and antikinks) with arbitrary
positions and relative velocities. During this phase, the exact vacuum
structure is irrelevant as the field fluctuations aren’t yet large enough to
‘feel’ the presence of the stabilizing $\lambda\phi^{4}$ interaction. This
corresponds to the weak-coupling, fast phase transition limit [37].

Figure 1:  Sketch of the wavefunctional for the field $\Psi(\phi)$ (dashed
blue) and the potential $V(\phi)$ (solid black) as a function of the field
$\phi$ after the $\mathbb{Z}_{2}$ symmetry breaking transition. Generically,
after $m^{2}(t)$ changes sign, $\Psi(\phi)$ is off-center with respect to the
maximum of the potential. (This is indeed the case when the bias arises from a
constant tilt, see the main text.) The wavefunction then evolves in time by
spreading and rolling, as indicated by the arrows. The DWs are field
configurations that interpolate between the two sides of the maximum. In
practice, we find ourselves in the standard notion of a population bias, with
an off-center initial probability distribution.

The sketch of how the computation proceeds is as follows. The starting point
is the quadratic theory (6) defined by two functions of time $m^{2}(t)$,
$\delta\phi(t)$. We assume periodic boundary conditions and discretize the
spatial coordinate. Effectively, this truncates to a finite number of harmonic
oscillators, which allows for a treatment in the Schrödinger picture. Since we
work at quadratic level, time evolution only amounts to keeping track of how
the Gaussian wavefunctional evolves. Following [36, 37], we identify a kink
number density operator, we extend it to area density for DWs, and simply
evaluate its expectation value.

The rest of this article is organized as follows. We review in Sec. II a close
analogue to DW networks, an ensemble of kinks and antikinks in $1+1$
dimensions, which even in the classical point-particle approximation displays
diffusive scaling. We discuss the quantum version of the problem in $1+1$
dimensions for kinks in Sec. III, and for (precursor) DWs in $d+1$ dimensions
in Sec. IV. We present the effective description of the precursor wall network
dynamics in the VOS ‘one-scale’ model language in V. We briefly discuss our
results in Sec. VI.

##  II Kink diffusion in $1+1$ dimensions

Kinks in $1+1$ dimensions differ qualitatively from DWs because they aren’t
extended objects and so they can only behave as massive point particles. As
such, a kink/antikink ensemble (or ‘plasma’) must not follow the self-similar
scaling common to true DWs. Instead, it falls under a scaling regime of
diffusive type. In the following we overview how this happens. Let us
emphasize that in this Section we consider stabilized kinks (as opposed to
kink precursors). Therefore we are having in mind some nonlinear potential for
the scalar, the details of which are unimportant beyond providing the kinks as
semiclassical solutions that interpolate between the two degenerate minima.
For simplicity we stick to a $\mathbb{Z}_{2}$ invariant model.

Both kinks and antikinks reduce to point particles of the same mass (by
$\mathbb{Z}_{2}$ symmetry). At large distances and late times, (anti-)kinks
are only characterized by the mass and the point particle approximation must
capture well the dynamics. This is the equivalent of the Nambu-Goto
description in $1+1$. Moreover, they have two more important properties: i)
they annihilate when they meet each other (also known as ballistic
annihilation), and ii) the $\mathbb{Z}_{2}$ symmetry breaking transition
creates an alternating array of kinks and antikinks. To a very good
approximation, the kink/antikink interaction is short range, so they only
annihilate when they step on each other. (This holds in the absence of a
pressure bias, which would introduce a long range constant force.)

The system then behaves similarly to an ‘ionised’ plasma, of kinks and
antikinks. The time evolution is almost trivial since the kink motions are
ballistic. One only needs to keep track of when kinks-antikinks collide, and
disappear thereafter. See e.g. Fig. 2 for an illustrative example of the
evolution of the kink-antikink plasma in this approximation.

This kind of problem appeared in the chemistry literature [38, 39], with a
focus on whether self-similar evolution and ‘universal’ exponents emerge. The
outcome is that the late time behaviour is sensitive to the initial state, the
initial distribution of velocities $v$ and inter-particle separations $\Delta
x$. This can be somewhat visualized already in Fig. 2: pairs that survive for
a long time must have small relative velocities. Then, a larger abundance of
small velocities in the initial state must translate into a slower dilution.

Figure 2: Evolution in time (horizontal axis) of the ‘kink plasma’ – a
distribution of kinks and anti-kinks in 1 space direction (vertical axis),
with random initial separations and velocities in a representative
realization. Consecutive lines correspond to alternating kinks/anti-kinks. In
the point-particle limit, kinks and antikinks annihilate when they meet each
other.

We can attempt to describe the statistical properties of the kink-antikink
plasma with a Boltzmann equation. In the simplest description, one keeps track
of the kink number density $n_{K}(t)$ and the (root-mean-squared) typical kink
velocity $v(t)$.

This is very much parallel to the so-called ‘velocity dependent one-scale’
(VOS) model, developed for cosmic strings [1] and later for DW networks [40],
to which we return in Sec.V. Then, the annihilation process is expected to be
captured by a simple ‘collision term’

$\dot{n}_{\rm K}=-\sigma\,v\,n_{\rm K}^{2}$ (7)

where $\sigma$ measures the scattering cross section for annihilation (which
is dimensionless in $1+1$).

If the system admits a self-similar solution one expects both $n$ and $v$ to
scale as power laws. Another way to look at (7) is that the collision rate
$\Gamma$ that keeps the plasma ‘equilibrated’ (rather, in a steady-state) is
$\Gamma\equiv-\dot{n}_{\rm K}\,/\,n_{\rm K}=\sigma\,v\,n_{\rm K}$. In a
scaling regime, it is reasonable to assume that $\dot{v}/v$ is also
proportional to $\Gamma$. This then provides an evolution equation for $v$,

$\dot{v}=-\kappa\,n_{\rm K}\,v^{2}$ (8)

with $\kappa$ another dimensionless constant.

By construction, (7) and (8) must contain scaling solutions. They are
parametized by the two constants: the value of $n_{\rm K}\,v\,t$, which must
be constant by dimensional analysis [38]; and the exponent $\mu$, entering as

$n_{\rm K}\sim t^{-\mu}\,,\qquad v\sim t^{\mu-1}~{}.$

For general $\sigma$, $\kappa$ one indeed finds nontrivial scaling solutions
with

$\mu=\frac{\sigma}{\kappa+\sigma}~{}.$

The non-universality found in annihilation kinetics [38, 39] is the statement
that the values of $\sigma$ and $\kappa$ depend on the initial state. For
$\kappa=\sigma$, Eqns. (7) and (8) give the diffusive scaling where the
separation between pairs, $L\equiv 1/n_{\rm K}$, exhibits Brownian-motion
growth

$L\sim t^{1/2}~{},$

typical of diffusion processes.

Interestingly, this is precisely the behaviour (5) found for precursor kinks
in $1+1$ using QFT methods [36, 37]. Note that precursor kinks do not
necessarily move ballistically. However they also behave like heavy localized
objects moving inertially. For this reason, and taking also into account that
for precursor kinks the initial condition is fixed by the physical vacuum
state of the quantum field [36, 37], it is not so surprising that the
precursor kinks don’t differ that much from stabilized kinks and that they
actually reduce to a particular choice of $\kappa$, $\sigma$. For real DWs
(kinks in $d+1$ dimensions with $d>1$) instead we expect that stabilized walls
and precursor walls behave differently, see Sec. V.

As a final comment, we just remark that at least qualitatively the effect a
(pressure) bias is also easy to visualize for the kink-antink ensemble.
Indeed, this introduces a constant proper acceleration $a$ that pushes kinks
towards antikinks (so that the true vacuum occupies more volume). One expects
that this is more efficient than just by diffusion when $\Gamma\lesssim a$.
Clearly, this has catastrophic consequences. At most, after a time of order
$1/a$ the velocities become relativistic and the encounter time then is only
dictated by the kink-antikink separation. It is clear then that the
(exponentially suppressed) abundance of kinks at late times strongly depends
on the distribution of separations near the annihilation time $t_{ann}$. We
expect that this feature holds for DWs in cosmology too. That is, that the
abundance of DWs after annihilation be set by the distribution of closed DWs
larger than the Hubble length at $t_{ann}$.

##  III Kink precursors and scaling from free QFT

We consider a scalar field $\phi(t,x)$ in $1+1$ dimensions undergoing a
quantum phase transition. Assuming a potential of the form (6), the model is
parameterized by 2 functions: $m^{2}(t)$ and $\delta\phi(t)$. To set the
stage, we assume that initially the scalar potential is
$\mathbb{Z}_{2}$-symmetric and has a unique vacuum at $\phi=0$ (where the mass
of $\phi$ is $m_{0}$). As time goes on, the potential flips (under the effect
of some external time-varying parameter for instance) and it acquires two
local minima. To keep it quite general, we allow the local maximum of this new
potential (where the tachyonic mass of $\phi$ is assumed to be of magnitude
$m_{0}$) to be shifted by some amount $\delta\phi_{0}$. Without loss of
generality we also assume that the potential vanishes at this point. The
$\mathbb{Z}_{2}$-symmetry is thus both explicitly and spontaneously broken
during the phase transition.

This stage of the phase transition can thus described by the quadratic
Lagrangian

$\mathcal{L}=\int
dx\left[\frac{1}{2}\dot{\phi}^{2}-\frac{1}{2}\phi^{\prime}{}^{2}-\frac{1}{2}m^{2}(t)(\phi-\delta\phi(t))^{2}\right]\,,$
(9)

where $m^{2}(t)$ and $\delta\phi(t)$ are functions of time only that verify,
$m^{2}=m_{0}^{2}$ and $\delta\phi=0$ at $t\to-\infty$, while
$m^{2}=-m_{0}^{2}$ and $\delta\phi=\delta\phi_{0}$ at $t\to+\infty$. This can
be achieved for instance by taking
$m^{2}(t)=-m_{0}^{2}\tanh\left(t/\tau_{m}\right)$ and
$\delta\phi(t)=\delta\phi_{0}\left[1+\tanh\left((t-t_{s})/\tau_{s}\right)\right]$
where the $t_{s}$, $\tau_{m}$ and $\tau_{s}$ are time scales that parametrize
the details of the quench (we shall focus on the instantaneous quench limit
$\tau_{m,\,s}\to 0$ with $t_{s}=0$ later on). Dotted and primed quantities
respectively denote time and space derivatives.

In this context, the natural definition of kinks are as the zeros of
$\phi-\delta\phi$. (We work at quadratic level, so keep in mind that these are
kink precursors really.) Our task then is to estimate the evolution of the
average kink number density over the course of the phase transition. Since the
model simply corresponds to a free quantum scalar field evolving in a time
dependent homogeneous background, we will proceed by fully solving for its
quantum dynamics in the Schrödinger picture. We will closely follow the
methods of [36, 37], by solving for the time-dependent wavefunctional
describing the quantum state of the field $\phi$ in the Schrödinger picture
and using this to semi-analytically compute the quantum average of a properly
defined kink number density operator.

###  III.1 Dynamics in the Schrödinger picture

We start by compactifying the model i.e. imposing periodic boundary conditions
at $x=0$ and $x=L$. After an integration by parts, (9) can thus be recast as

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

We then discretize it on a finite lattice made up of $N$ evenly spaced points
labelled by an index $\mathscr{i}$ running from $1$ to $N$. The lattice
spacing $L/N$ is denoted by $a$. If we define the discretized field values
$\phi_{\mathscr{i}}\equiv\phi(t,\mathscr{i}a)$ and replace the second
derivative by its lowest order central finite difference approximation,

$\phi_{\mathscr{i}}^{\prime\prime}\to\frac{\phi_{\mathscr{i}-1}-2\phi_{\mathscr{i}}+\phi_{\mathscr{i}+1}}{a^{2}}\,,$
(11)

we can rewrite (10) as

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\sum_{\mathscr{i}=1}^{N}a\left[\frac{1}{2}\dot{\phi}_{\mathscr{i}}^{2}+\frac{1}{2a^{2}}\phi_{\mathscr{i}}(\phi_{\mathscr{i}-1}-2\phi_{\mathscr{i}}+\phi_{\mathscr{i}+1})\right.$
(12) $\displaystyle\hskip
85.35826pt\left.-\frac{1}{2}m^{2}(t)(\phi_{\mathscr{i}}-\delta\phi(t))^{2}\right]\,.$

Notice that $\phi_{N+1}=\phi_{1}$ and $\phi_{0}=\phi_{N}$ by virtue of the
periodicity of the lattice.

We can give a more compact expression of this Lagrangian by assembling the
discretized field values in a column vector
$\bm{\phi}\equiv(\phi_{1},\phi_{2},\dots,\phi_{N})^{T}$ and defining the
matrix $\bm{\Omega}^{2}(t)$ by

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

(Here and henceforth, bold quantities denote vectors and matrices while T
denotes matrix transposition.) With these definitions and conventions Eq. (12)
reads

$\mathcal{L}=\frac{a}{2}\dot{\bm{\phi}}^{T}.\dot{\bm{\phi}}-\frac{a}{2}(\bm{\phi}-\delta\phi(t)\bm{1})^{T}.\bm{\Omega_{2}}(t).(\bm{\phi}-\delta\phi(t)\bm{1})\,,$
(14)

where we have introduced the ‘vector of ones’ $\bm{1}\equiv(1,1,\dots,1)^{T}$
and used the property $\bm{\Omega}^{2}(t).\bm{1}=m^{2}(t)\bm{1}$ (i.e.,
$\bm{1}$ is an eigenvector associated to the eigenvalue $m^{2}(t)$).

Since the Lagrangian is quadratic, and assuming that at $t=t_{0}\ll-
m_{0}^{-1}$ the field starts in its (gaussian) quantum vacuum, we know that
the state will remain gaussian during time evolution. A good ansatz for the
Schrödinger wavefunctional describing the state of the system at time $t$ is
therefore

$\Psi(\bm{\phi},t)=\mathcal{N}(t)\exp\left[ia\bm{D}(t)^{T}.\bm{\phi}+\frac{ia}{2}\bm{\phi}^{T}.\bm{M}(t).\bm{\phi}\right]\,,$
(15)

where $\bm{M}(t)$ is an $N\times N$ complex symmetric matrix, $\bm{D}(t)$ a
complex $N$-component column vector and $\mathcal{N}(t)$ a normalization
factor. The requirement that the field $\phi$ be in its vacuum long before the
phase transition imposes $\bm{M}(t_{0})=i\bm{\Omega}^{2}(t_{0})^{1/2}$,
$\bm{D}(t_{0})=0$, and
$\mathcal{N}(t_{0})=(a/\pi)^{N/4}\det(\bm{\Omega}^{2}(t_{0}))^{1/8}$. Here the
matrix powers are unambiguously defined since the matrix
$\bm{\Omega}^{2}(t_{0})$ is symmetric positive definite.

Now, $\Psi$ verifies the Schrödinger equation

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

where we have set $\hbar=1$. Plugging (15) into (16) and identifying the
quadratic terms in $\phi_{\mathscr{i}}\phi_{\mathscr{j}}$ yields the matrix
dfiferential equation

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

It is easy to convince oneself that the unique solution to this equation
subject to the specified initial conditions is symmetric for all times.
Analogously, identifying the linear terms in $\phi_{\mathscr{i}}$ (and using
the fact that $\bm{1}$ is an eigenvector of $\bm{\Omega}^{2}$) yields

$\dot{\bm{D}}+\bm{M}.\bm{D}-m^{2}(t)\delta\phi(t)\bm{1}=0\,.$ (18)

The normalization factor is also easy to compute as a function of $\bm{M}$ but
we will not need it explicitly. It turns out that one can solve for $\bm{M}$
and $\bm{D}$ in terms of $N$ complex functions of time.

We now introduce the complex mode functions $c_{n}(t)$ (with $n$ running from
$0$ to $N-1$) verifying

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

and with initial conditions

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

Then it is straightforward to verify that

$\displaystyle[\bm{M}(t)]_{\mathscr{i}\mathscr{j}}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=0}^{N-1}c_{n}(t)^{-1}\dot{c}_{n}(t)e^{2i\pi
n(\mathscr{i}-\mathscr{j})/N}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=0}^{N-1}c_{n}(t)^{-1}\dot{c}_{n}(t)\cos(2\pi
n(\mathscr{i}-\mathscr{j})/N)\,,$

where we have used the fact that $c_{n}(t)=c_{N-n}(t)$ for $1\leq n\leq N-1$
to write $\bm{M}$ in manifestly symmetric form. Similarly we can see that

$\bm{D}(t)=c_{0}(t)^{-1}\int_{t_{0}}^{t}ds\,m^{2}(s)\,\delta\phi(s)c_{0}(s)\bm{1}\,.$
(23)

We are finally in a position to write the probability density functional of
the field configuration $\bm{\phi}$,
$\mathcal{P}(\bm{\phi},t)\equiv|\Psi(\bm{\phi},t)|^{2}$, in terms of the
$c_{n}(t)$. We have

$\mathcal{P}(\bm{\phi},t)=\frac{1}{\sqrt{\det(2\pi\bm{K})}}e^{-(\bm{\phi}-\Delta(t)\bm{1})^{T}.\bm{K}(t)^{-1}.(\bm{\phi}-\Delta(t)\bm{1})/2}\,,$
(24)

where

$[\bm{K}(t)]_{\mathscr{i}\mathscr{j}}=\frac{1}{N}\sum_{n=0}^{N-1}|c_{n}(t)|^{2}\cos(2\pi
n(\mathscr{i}-\mathscr{j})/N)$ (25)

is the spatial correlation function between two lattice points, and

$\displaystyle\Delta(t)$ $\displaystyle=$ $\displaystyle
ia\int_{t_{0}}^{t}ds\,m^{2}(s)\,\delta\phi(s)$ (26) $\displaystyle\hskip
34.14322pt\times\left(c_{0}(t)^{*}c_{0}(s)-c_{0}(t)c_{0}(s)^{*}\right)$

represents the (homogeneous) shift of the vacuum expectation value of $\phi$
triggered by the potential bias $\delta\phi$. To derive these expressions we
have used (III.1) and (23) as well as the fact that
$c_{n}(t)^{*}\dot{c}_{n}(t)-c_{n}(t)\dot{c}_{n}(t)^{*}$ is a conserved
quantity equal to $i/a$ for all $n$. As a sanity check, it is easy to verify
that for constant positive $m^{2}(t)=m_{0}^{2}$ and slowly varying
$\delta\phi(t)$, $\Delta(t)\approx\delta\phi(t)$.

We now have all the necessary tools at our disposal to evaluate the kink
number density.

###  III.2 Kink number density

We will follow the method in Refs. [36, 37] and look for kinks and antikinks
among zeros of $\phi-\delta\phi(t)$. Of course some of those will be quantum
fluctuations and there will be a vast overcounting but we will turn to that
problem later on. For the time being we define the quantum operator

$\displaystyle\hat{n}_{\rm Z}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{L}\sum_{\mathscr{i}=1}^{N}\frac{1}{4}\left[\rm{sgn}(\hat{\phi}_{\mathscr{i}}-\delta\phi(t))-\rm{sgn}(\hat{\phi}_{\mathscr{i}+1}-\delta\phi(t))\right]^{2}$
(27) $\displaystyle=$
$\displaystyle\frac{N}{2L}-\frac{1}{2L}\sum_{\mathscr{i}=1}^{N}\rm{sgn}\left((\hat{\phi}_{\mathscr{i}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1}-\delta\phi)\right)\,,$

where $\rm{sgn}$ denotes the signum function. This operator counts the number
density of $\phi-\delta\phi(t)$ sign changes between two consecutive lattice
points thus providing a lower estimate for the number density of zeros of this
quantity. Of course this estimate will change as the lattice gets finer (or in
other words as $a$ becomes smaller) and, given the fundamentally quantum
nature of the problem, may potentially diverge in the continuum limit. We will
come back to this problem and for the time being will simply disregard any
subtleties related to the coarseness of the lattice. In fact we want to
calculate the quantum average of this operator in the state whose
wavefunctional we computed previously. This reads

$\langle\hat{n}_{\rm
Z}\rangle=\frac{N}{2L}-\frac{1}{2L}\sum_{\mathscr{i}=1}^{N}\left\langle\rm{sgn}\left((\hat{\phi}_{\mathscr{i}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1}-\delta\phi)\right)\right\rangle\,,$
(28)

where, after a shift in the integration variables,

$\displaystyle\left\langle\rm{sgn}\left((\hat{\phi}_{\mathscr{i}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1}-\delta\phi)\right)\right\rangle=$
$\displaystyle\hskip 14.22636pt\int d\phi_{1}\dots
d\phi_{N}\rm{sgn}\left(\phi_{\mathscr{i}}\phi_{\mathscr{i}+1}\right)\mathcal{P}(\bm{\phi}+\delta\phi(t)\mathbf{1},t)\,.$
(29)

Now, since the matrix $\bm{K}^{-1}$ is circulant i.e. its coefficients
$[\bm{K}^{-1}]_{\mathscr{i}\mathscr{j}}$ only depend on
$\mathscr{i}-\mathscr{j}$, it is easy to see that
$\left\langle\rm{sgn}\left((\hat{\phi}_{\mathscr{i}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1}-\delta\phi)\right)\right\rangle=\left\langle\rm{sgn}\left((\hat{\phi}_{1}-\delta\phi)(\hat{\phi}_{2}-\delta\phi)\right)\right\rangle$.
(This is the algebraic signature of the fact that the background is
translationally invariant.) If we define the four quadrants of the
$(\phi_{1},\phi_{2})$ plane in the conventional way, Eq. 28 reduces to

$\langle\hat{n}_{\rm
Z}\rangle=\frac{N}{2L}\left[1+\sum_{Q=1}^{4}(-1)^{Q}I_{Q}\right]\,,$ (30)

where

$I_{Q}\equiv\iint_{Q\
\rm{quadrant}}d\phi_{1}d\phi_{2}\tilde{\mathcal{P}}(\phi_{1}+\delta\phi(t),\phi_{2}+\delta\phi(t))$
(31)

and

$\tilde{\mathcal{P}}(\phi_{1},\phi_{2})\equiv\int d\phi_{3}\dots
d\phi_{N}\mathcal{P}(\bm{\phi},t)$ (32)

is the marginal probability density of $\phi_{1}$ and $\phi_{2}$. We can
easily compute this quantity by using a well-known property of multi-variate
normal distributions such as $\mathcal{P}$ which states that in order to
obtain the marginal distribution over a subset of variables, one simply needs
to drop the variables that are intregrated out from the covariance matrix and
the mean vector (see for instance [41]). In our case the mean vector is zero
and the marginal covariance matrix reduces to the $2\times 2$ upper left block
of $\bm{K}$ or $\bm{K}_{2\times 2}\equiv\alpha\bm{I}_{2}+\beta\bm{\sigma}_{1}$
where

$\displaystyle\alpha(t)$ $\displaystyle\equiv$
$\displaystyle[\bm{K}(t)]_{11}=\frac{1}{N}\sum_{n=0}^{N-1}|c_{n}(t)|^{2}\,,$
(33) $\displaystyle\beta(t)$ $\displaystyle\equiv$
$\displaystyle[\bm{K}(t)]_{12}=\frac{1}{N}\sum_{n=0}^{N-1}|c_{n}(t)|^{2}\cos(2\pi
n/N)\,.$ (34)

With these notations we obtain

$\displaystyle\tilde{\mathcal{P}}(\phi_{1},\phi_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\sqrt{\alpha^{2}-\beta^{2}}}\exp\biggl{[}-\frac{1}{2(\alpha^{2}-\beta^{2})}\Big{(}\alpha(\phi_{1}-\Delta)^{2}$
$\displaystyle\hskip
22.76228pt-2\beta(\phi_{1}-\Delta)(\phi_{2}-\Delta)+\alpha(\phi_{2}-\Delta)^{2}\Big{)}\biggr{]}\,.$

Noticing that $\sum_{Q=1}^{4}I_{Q}=1$ and that $I_{2}=I_{4}$ we can now
rewrite (30) in the more explicit manner,

$\langle\hat{n}_{\rm
Z}\rangle=\frac{2N}{L}I_{2}=\frac{N}{L\sqrt{\pi}}\int_{0}^{\infty}dx\,e^{-x^{2}}\left[{\rm
erf}\left(x\sqrt{\frac{\alpha-\beta}{\alpha+\beta}}-\frac{\delta\phi-\Delta}{\sqrt{\alpha+\beta}}\right)+{\rm
erf}\left(x\sqrt{\frac{\alpha-\beta}{\alpha+\beta}}+\frac{\delta\phi-\Delta}{\sqrt{\alpha+\beta}}\right)\right]\,,$
(35)

where ${\rm erf}(z)=2/\sqrt{\pi}\int_{0}^{z}dt\,e^{-t^{2}}$ is the standard
error function. As stated at the beginning of this section, if we simply use
this formula as is, we will vastly overcount the number density of kinks on
the lattice. This is due to small fluctuations of the quantum field that are
always present (whether there is a phase transition or not) and that should be
disregarded. Following Refs. [36, 37] we notice that such spurious zeros of
the field that do not correspond to a kink or antikink are due to the presence
of oscillatory modes in the expressions for $\alpha$ and $\beta$. We are thus
led to define $\bar{\alpha}$ and $\bar{\beta}$ by restricting the sums to
those modes $c_{n}$ that are unstable i.e. with $n$ such that $4\sin^{2}(\pi
n/N)/a^{2}+m^{2}(t)\leq 0$,

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

where $n_{c}(t)\equiv\lfloor N\sin^{-1}(a\sqrt{-m^{2}(t)}/2)/\pi\rfloor$ and
$\lfloor\rfloor$ denotes the integer part function. (When there are no
unstable modes, $\bar{\alpha}$ and $\bar{\beta}$ are understood to vanish by
convention.) Then the average number density of kinks $n_{\rm K}$ can be
computed from $\langle\hat{n}_{\rm Z}\rangle$ by replacing $\alpha$ and
$\beta$ by $\bar{\alpha}$ and $\bar{\beta}$ in (35),

$\displaystyle n_{\rm K}$ $\displaystyle=$
$\displaystyle\frac{N}{L\sqrt{\pi}}\int_{0}^{\infty}dx\,e^{-x^{2}}\left[{\rm
erf}\left(x\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}-\frac{\delta\phi-\Delta}{\sqrt{\bar{\alpha}+\bar{\beta}}}\right)\right.$
(38) $\displaystyle\hskip 42.67912pt\left.+{\rm
erf}\left(x\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}+\frac{\delta\phi-\Delta}{\sqrt{\bar{\alpha}+\bar{\beta}}}\right)\right]\,.$

We can already notice that for zero bias, $\delta\phi(t)=\Delta(t)=0$ and we
recover the result of Refs. [36, 37],

$\left.n_{\rm K}\right|_{\rm no\,bias}=\frac{2N}{\pi
L}\tan^{-1}\left(\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}\right)=\frac{N}{\pi
L}\cos^{-1}\left(\frac{\bar{\beta}}{\bar{\alpha}}\right)\,.$ (39)

Before going any further it will be interesting to see how (38) and (39)
simplify in the continuum limit i.e. when $N\to\infty$ at fixed $L$. For this
we notice that the quantity

$\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}=\tan\left(\frac{\pi
L}{2N}\left.n_{\rm K}\right|_{\rm no\,bias}\right)\approx\frac{\pi
L}{2N}\left.n_{\rm K}\right|_{\rm no\,bias}$ (40)

vanishes in this limit since $\left.n_{\rm K}\right|_{\rm no\,bias}$ is a
physical quantity that has a well-defined, $N$ independent, finite limit for
all times [36, 37]. Moreover, the integral in (38) is dominated by values of
the integrand close to $x=0$ (because of the $e^{-x^{2}}$ factor). We can
therefore Taylor expand the error functions in powers of the vanishingly small
quantity $x\sqrt{(\bar{\alpha}-\bar{\beta})/(\bar{\alpha}+\bar{\beta})}$ to
obtain

$n_{\rm K}=\left.n_{\rm K}\right|_{\rm
no\,bias}\exp\left[-\frac{(\delta\phi-\Delta)^{2}}{\bar{\alpha}+\bar{\beta}}\right]\,.$
(41)

It is worth noticing that this expression only depends on $L$ through
$\bar{\alpha}$ and $\bar{\beta}$ which are sums of a finite number
($2n_{c}(t)+1\leq L\sqrt{-m^{2}(t)}/4$) of terms. As $L\to\infty$ these sums
become integrals.

We now turn to the evaluation of $n_{\rm K}$ which can be done either semi-
analytically in the particular case of a sudden phase transition, where
$m^{2}(t)=-m_{0}^{2}(2\Theta(t)-1)$ and
$\delta\phi(t)=\delta\phi_{0}\Theta(t)$ ($\Theta$ being the standard Heaviside
step function). In this case the mode functions $c_{n}(t)$ are exactly
calculable and, taking first the continuum limit $N\to\infty$ and then the
infinite volume limit $L\to\infty$,we obtain

$\displaystyle\delta\phi(t)-\Delta(t)=\delta\phi_{0}\cosh\left[m_{0}(t-t_{s})\right]\Theta(t-t_{s})\,,$
(42)
$\displaystyle\bar{\alpha}+\bar{\beta}=\frac{1}{\pi}\bigintsss_{0}^{m_{0}}dk\,\left[\frac{m_{0}^{2}\cosh\left[2t\sqrt{m_{0}^{2}-k^{2}}\right]-k^{2}}{(m_{0}^{2}-k^{2})\sqrt{m_{0}^{2}+k^{2}}}\right]$
(43)

and $\left.n_{\rm K}\right|_{\rm no\,bias}$ is known from Eq. (93) of Ref.
[37] to be

$\displaystyle\frac{1}{\pi}\left\\{\bigintsss_{0}^{m_{0}}dk\,k^{2}\left[\frac{m_{0}^{2}\cosh\left(2t\sqrt{m_{0}^{2}-k^{2}}\right)-k^{2}}{(m_{0}^{2}-k^{2})\sqrt{k^{2}+m_{0}^{2}}}\right]\right\\}^{1/2}$
$\displaystyle\hskip
14.22636pt\times\left\\{\bigintsss_{0}^{m_{0}}dk\,\left[\frac{m_{0}^{2}\cosh\left(2t\sqrt{m_{0}^{2}-k^{2}}\right)-k^{2}}{(m_{0}^{2}-k^{2})\sqrt{k^{2}+m_{0}^{2}}}\right]\right\\}^{-1/2}\,.$
(44)

Now, plugging these expression into (41) yields an explicit analytic
expression that we plot in Fig. 3. As expected we notice that the initial kink
number density is orders of magnitude lower than in the zero bias case. This
is due to the presence of a potential barrier that suppresses kink formation.
A surprising feature of this plot is the small bump appearing around $t\sim
m_{0}^{-1}$ and that decays rapidly afterwards. Technically, this is due to
the vacuum expectation value of the field $\Delta(t)$ growing more slowly than
the “width” $\sim\bar{\alpha}+\bar{\beta}$ of the probability density
functional. Physically, it seems that we can interpret this are bubble
nucleation in the first stages of evolution, when the potential barrier is
still not too high and tunnelling shouldbe not suppressed. The late time
behavior is indicative of a faster than power law suppression. In fact, in
this limit, $n_{\rm K}/\left.n_{\rm K}\right|_{\rm no\,bias}$ is well fit by
$\exp[-1.72\,t^{0.51}]$ in units where $m_{0}=1$, as can be seen in Fig. 4.

Figure 3:  Kink number density as a function of time in $m_{0}=1$ units for
different values of the bias $\delta\phi_{0}$, in $1+1$ dimensions. The
topmost curve corresponds to the zero bias case, and $\delta\phi_{0}$
increases in increments of $0.2$ as we move down to the lower curves.  Figure
4:  Plot of the natural logarithm of the kink number density as a function of
the square root of time in $m_{0}=1$ units and for $\delta\phi_{0}=1$.

We can even recover explicitly the late time behavior of $n_{\rm K}(t)$.
Realizing that in the limit $t\gg m_{0}$ all the $k$ integrals in (43) and
(44) are dominated by values of $k\ll m_{0}$ we can replace the integrands
with their lowest order $k/m_{0}$ expansions to obtain

$n_{\rm K}(t)\sim\sqrt{\frac{m_{0}}{t}}e^{-\delta\phi_{0}^{2}\sqrt{\pi
m_{0}t}}\,.$ (45)

This is in good agreement with the above numerically determined coefficients.

##  IV Domain wall networks

The same method outlined above for kinks in $1+1$ dimensions can be extended
rather directly to $d+1$ dimensions. For concreteness we will show explicitly
how the extension works in $2+1$ dimensions and then give the general result.

We start with the Lagrangian for the $2+1$ dimensional real scalar field
$\phi(t,x,y)$

$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle\int
dxdy\left[\frac{1}{2}(\partial_{t}\phi)^{2}-\frac{1}{2}(\partial_{x}\phi)^{2}-\frac{1}{2}(\partial_{y}\phi)^{2}\right.$
(46) $\displaystyle\hskip
56.9055pt\left.-\frac{1}{2}m^{2}(t)(\phi-\delta\phi(t))^{2}\right]\,,$

where the functions of time $m^{2}(t)$ and $\delta\phi(t)$ obey the same
properties as in the $1+1$ dimensional case. (The field $\phi$ is now
dimnensionful however.) This is seen to be a model with a broken
$\mathbb{Z}_{2}$ symmetry that typically would feature domain walls. Next, we
compactify space on a 2-torus of area $L^{2}$ by assuming periodic boundary
conditions ($\phi(x+L,y)=\phi(x,y+L)=\phi(x,y)$) and discretize it on a
regular square lattice made up of $N^{2}$ (with lattice spacing $a=L/N$). At
each lattice point
$(x_{\mathscr{i}},y_{\mathscr{j}})\equiv(\mathscr{i}a,\mathscr{j}a)$ we define
the discretized field values
$\phi_{\mathscr{i}\mathscr{j}}\equiv\phi(x_{\mathscr{i}},y_{\mathscr{j}})$.

The discretized version of (46) has a form analogous to (14),

$\mathcal{L}=\frac{a^{2}}{2}\dot{\bm{\phi}}^{T}.\dot{\bm{\phi}}-\frac{a^{2}}{2}(\bm{\phi}-\delta\phi(t)\bm{1})^{T}.\bm{\Omega_{2}}(t).(\bm{\phi}-\delta\phi(t)\bm{1})\,,$
(47)

as long as it is understood that any vectors and matrices are now $N^{2}$ and
$N^{2}\times N^{2}$ dimensional respectively. For instance,

${\bm{\phi}}\equiv(\phi_{11},\phi_{12},...,\phi_{1N},\phi_{21},...,\phi_{2N},...,\phi_{NN-1},\phi_{NN})^{T}.$
(48)

and $\bm{\Omega}^{2}(t)$ is given by

$[\bm{\Omega}^{2}]_{\mathscr{i}\mathscr{j},\mathscr{k}\mathscr{l}}=\begin{cases}+{2}/{a^{2}}+m^{2}(t)\,,&\mathscr{i}=\mathscr{k},\mathscr{j}=\mathscr{l}\\\
-{1}/{a^{2}}\,,&\mathscr{i}=\mathscr{k}\pm 1,\mathscr{j}=\mathscr{l}\pm 1\\\
0\,,&\text{otherwise}\,,\end{cases}$ (49)

where equality relations are understood to modulo $N$. More generally, any
$N^{2}\times N^{2}$ matrix $\bm{A}$ would be represented by a two-dimensional
array of matrix elements $A_{ij,kl}$ arranged in the following way:

$A=\pNiceMatrix[xdots/shorten=0.6em]A_{11,11}&A_{11,12}\Cdots
A_{11,1N}A_{11,21}A_{11,22}\Cdots\\\ A_{12,11}A_{12,12}\Cdots
A_{12,1N}A_{12,21}A_{12,22}\Cdots\\\ \Vdots\Vdots\Vdots\Vdots\Vdots\\\
A_{1N,11}A_{1N,12}\Cdots A_{1N,1N}A_{1N,21}A_{1N,22}\Cdots\\\ A_{21,11}\\\
A_{22,11}\\\ \Vdots$ (50)

The generalization of the results of the previous sections is now
straightforward. The functional Schrödinger equation again has a Gaussian
solution $\Psi(t,\bm{\phi})$ which can be expressed in terms of the two-
deimansional mode functions $c_{n,m}(t)$ (with both $n$ and $m$ running from
$0$ to $N-1$). These verify

${\ddot{c}}_{n,m}+\left[k_{n,m}^{2}+m^{2}(t)\right]c_{n,m}=0\,,$ (51)

and with initial conditions

$\displaystyle c_{\mathscr{n},m}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{a\sqrt{2}}\left[k_{n,m}^{2}+m^{2}(t_{0})\right]^{-1/4}\,,$
(52) $\displaystyle{\dot{c}}_{n,m}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{i}{a\sqrt{2}}\left[k_{n,m}^{2}+m^{2}(t_{0})\right]^{1/4}\,.$
(53)

Here we have introduced the discretized momentum

$k_{n,m}=\frac{2}{a}\left\\{\sin^{2}\left(\frac{\pi
n}{N}\right)+\sin^{2}\left(\frac{\pi m}{N}\right)\right\\}^{1/2}\,,$ (54)

for notational simplicity. (Notice also that the normalization of the mode
functions has been modified with respect to (20) and (21) in such a way that
they have the same dimension as the field $\phi$ i.e. 1/2.) The probability
density functional $\mathcal{P}(t,\bm{\phi})=|\Psi(t,\bm{\phi})|^{2}$ is of
the same form as (24) but with

$\displaystyle[\bm{K}(t)]_{\mathscr{i}\mathscr{j},\mathscr{k}\mathscr{l}}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\sum_{n,m=0}^{N-1}|c_{n,m}(t)|^{2}$ (55)
$\displaystyle\times\cos\left[\frac{2\pi}{N}\left(n(\mathscr{i}-\mathscr{k})+m(\mathscr{j}-\mathscr{l})\right)\right]\,,$

and

$\displaystyle\Delta(t)$ $\displaystyle=$ $\displaystyle
ia^{2}\int_{t_{0}}^{t}ds\,m^{2}(s)\,\delta\phi(s)$ (56) $\displaystyle\hskip
22.76228pt\times\left(c_{0,0}(t)^{*}c_{0,0}(s)-c_{0,0}(t)c_{0,0}(s)^{*}\right)\,.$

This contains all the information about the quantum dynamics of the field
theory we are considering. In the case of domain walls in two or higher
dimensions, the relevant quantity to compute is the average area density
(strictly speaking length density in two dimensions) of such extended objects.
It remains then to write a quantum operator generalizing $\hat{n}_{\rm Z}$ and
which, in some limit, will describe the desired observable. A possible choice
is

$\displaystyle\hat{\mathcal{A}}_{Z}$ $\displaystyle\equiv$
$\displaystyle\frac{a}{L^{2}}\sum_{\mathscr{i,\mathscr{j}}=1}^{N}\frac{1}{4}\left[\rm{sgn}(\hat{\phi}_{\mathscr{i}\mathscr{j}}-\delta\phi(t))-\rm{sgn}(\hat{\phi}_{\mathscr{i}+1,\mathscr{j}}-\delta\phi(t))\right]^{2}+\frac{a}{L^{2}}\sum_{\mathscr{i,\mathscr{j}}=1}^{N}\frac{1}{4}\left[\rm{sgn}(\hat{\phi}_{\mathscr{i}\mathscr{j}}-\delta\phi(t))-\rm{sgn}(\hat{\phi}_{\mathscr{i},\mathscr{j}+1}-\delta\phi(t))\right]^{2}$
(57) $\displaystyle=$
$\displaystyle\frac{N}{L}-\frac{a}{2L^{2}}\sum_{\mathscr{i},\mathscr{j}=1}^{N}\rm{sgn}\left((\hat{\phi}_{\mathscr{i},\mathscr{j}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1,\mathscr{j}}-\delta\phi)\right)-\frac{a}{2L^{2}}\sum_{\mathscr{i},\mathscr{j}=1}^{N}\rm{sgn}\left((\hat{\phi}_{\mathscr{i},\mathscr{j}}-\delta\phi)(\hat{\phi}_{\mathscr{i},\mathscr{j}+1}-\delta\phi)\right)\,.$

Figure 5:  Sketch of a generic domain wall crossing two adjacent sides of a
lattice cell. The field vanishes along the lower and right sides, and the
length of the corresponding portion of domain wall is estimated by the
operator (57) to be equal to $2a$. Assuming the domain wall to be smooth on
scales of the size of the cell, this estimate may differ from the exact value
by a factor of at most $\sqrt{2}$.

This operator counts the number of sides of cells of our regular square
lattice that are traversed by a domain wall i.e. that are such that the
quantity $\phi_{\mathscr{i}\mathscr{j}}-\delta\phi$ changes sign along them,
multiplies the result by $a$ (the average length of domain wall traversing a
side) and divides it by the total area $L^{2}$ (see Fig. refplaquette). Of
course it suffers from the same overcounting and undercounting disadvantages
as the operator $\hat{n}_{\rm Z}$ but, on top of that, its value can only be
trusted up to factors of order 1 since a domain wall that traverses cell
diagonally will contribute a length $2a$ instead of $a\sqrt{2}$. The
translational and rotational symmetries of the theory (more precisely, the
residual symmetries of the discretized theory: discrete translations and
rotations by multiples of $\pi/4$) imply that

$\displaystyle\left\langle\rm{sgn}\left((\hat{\phi}_{\mathscr{i},\mathscr{j}}-\delta\phi)(\hat{\phi}_{\mathscr{i}+1,\mathscr{j}}-\delta\phi)\right)\right\rangle$
(58)
$\displaystyle=\left\langle\rm{sgn}\left((\hat{\phi}_{\mathscr{i},\mathscr{j}}-\delta\phi)(\hat{\phi}_{\mathscr{i},\mathscr{j}+1}-\delta\phi)\right)\right\rangle$
$\displaystyle=\left\langle\rm{sgn}\left((\hat{\phi}_{11}-\delta\phi)(\hat{\phi}_{12}-\delta\phi)\right)\right\rangle\,,$

and thus the average value of the operator $\hat{\mathcal{A}}_{Z}$ simplifies
considerably:

$\left\langle\hat{\mathcal{A}}_{Z}\right\rangle=\frac{N}{L}\left[1-\left\langle\rm{sgn}\left((\hat{\phi}_{11}-\delta\phi)(\hat{\phi}_{12}-\delta\phi)\right)\right\rangle\right]\,.$
(59)

It is clear from this point onward that the computation will proceed along the
same lines as in Sec. III. With the obvious replacements of
$\phi_{1}\rightarrow\phi_{11}$, $\phi_{2}\rightarrow\phi_{12}$ and the
redefinition of

$\displaystyle\alpha(t)$ $\displaystyle\equiv$
$\displaystyle[\bm{K}(t)]_{11,11}=\frac{1}{N^{2}}\sum_{n,m=0}^{N-1}|c_{n,m}(t)|^{2}\,,$
(60) $\displaystyle\beta(t)$ $\displaystyle\equiv$
$\displaystyle[\bm{K}(t)]_{11,12}=\frac{1}{N^{2}}\sum_{n,m=0}^{N-1}|c_{n,m}(t)|^{2}$
(61) $\displaystyle\hskip 99.58464pt\times\cos(2\pi m/N)\,,$

we can obtain analogous equations to those in (30) and (35). Introducing the
cutoff versions of $\alpha$ and $\beta$ i.e. restricting the sums to those
modes with negative frequency yields $\bar{\alpha}$ and $\bar{\beta}$ and
allows us to write the average domain wall area density as

$\displaystyle\mathcal{A}_{DW}$ $\displaystyle=$
$\displaystyle\frac{N}{2L\sqrt{\pi}}\int_{0}^{\infty}dx\,e^{-x^{2}}\left[{\rm
erf}\left(x\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}-\frac{\delta\phi-\Delta}{\sqrt{\bar{\alpha}+\bar{\beta}}}\right)\right.$
(62) $\displaystyle\hskip 42.67912pt\left.+{\rm
erf}\left(x\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}+\frac{\delta\phi-\Delta}{\sqrt{\bar{\alpha}+\bar{\beta}}}\right)\right]\,.$

The no-bias case is analytically integrable and gives

$\left.\mathcal{A}_{DW}\right|_{\rm no\,bias}=\frac{N}{\pi
L}\tan^{-1}\left(\sqrt{\frac{\bar{\alpha}-\bar{\beta}}{\bar{\alpha}+\bar{\beta}}}\right)\,,$
(63)

while, in the $N\rightarrow\infty$ limit, the general result reads

$\mathcal{A}_{DW}=\left.\mathcal{A}_{DW}\right|_{\rm
no\,bias}\exp\left[-\frac{(\delta\phi-\Delta)^{2}}{\bar{\alpha}+\bar{\beta}}\right]\,.$
(64)

We can also give an analytical estimate of the late time behavior of the
average area density of domain walls for the case of a sudden phase transition
($m^{2}(t)=-m_{0}^{2}(2\Theta(t)-1)$ and
$\delta\phi(t)=\delta\phi_{0}\Theta(t)$) in the limit of infinite volume
$L\rightarrow\infty$, where

$\displaystyle\delta\phi(t)-\Delta(t)=\delta\phi_{0}\cosh\left(m_{0}t\right)\Theta(t)\,,$
(65)
$\displaystyle\bar{\alpha}+\bar{\beta}=\frac{1}{2\pi}\bigintsss_{0}^{m_{0}}kdk\left[\frac{m_{0}^{2}\cosh\left[2t\sqrt{m_{0}^{2}-k^{2}}\right]-k^{2}}{(m_{0}^{2}-k^{2})\sqrt{m_{0}^{2}+k^{2}}}\right]\,,$
(66)
$\displaystyle\bar{\alpha}-\bar{\beta}=\frac{a^{2}}{16\pi}\bigintsss_{0}^{m_{0}}k^{3}dk\left[\frac{m_{0}^{2}\cosh\left[2t\sqrt{m_{0}^{2}-k^{2}}\right]-k^{2}}{(m_{0}^{2}-k^{2})\sqrt{m_{0}^{2}+k^{2}}}\right]\,.$
(67)

Indeed, using the above expressions and taking the late time limit we find
that

$\mathcal{A}_{DW}(t)\sim\sqrt{\frac{m_{0}}{t}}\;\exp\left[-2\pi
t\,\delta\phi_{0}^{2}\right]\,.$ (68)

In $d+1$ dimensions (where the above results can be extended with minimal
modifications) one would get

$\mathcal{A}_{DW}^{d+1}(t)\sim\sqrt{\frac{m_{0}}{t}}\;\exp\left[-{\cal
C}_{d}\,\delta\phi_{0}^{2}m_{0}\left(\frac{t}{m_{0}}\right)^{d/2}\right]\,.$
(69)

with

${\cal C}_{d}=2^{d-1}\,\pi^{d/2}\,.$ (70)

Interestingly, for $d=3$ this expression coincides with the result obtained
for ‘phase ordering’ in condensed matter [24]. We comment on why in Sec. VI.

##  V VOS model

The precursor domain walls introduced above are related to the standard DWs
created during discrete symmetry breaking transitions, but they differ in an
important aspect. Precursor walls are simply the zero iso-surfaces of free
tachyonic fluctuating fields and so they don’t obey the Nambu Goto (NG)
equation. The results from Sec. IV imply that these precursor wall network can
also enter a self-similar regime, yet with different properties compared to
the standard DW network scaling. Let us try to understand the difference in
the language of the effective velocity-dependent one-scale (VOS) models.

First, recall the VOS model for standard DW networks [40]. It consists in a
simplified description of the DW network in terms of two quantities: $L(t)$,
the correlation length or the average separation between walls; and $v(t)$ the
(root-mean-squared) velocity of the walls. Using energy conservation arguments
together with general properties of the Nambu-Goto equation of motion [40, 42]
for DWs in $d+1$ dimensions, one arrives at

$\displaystyle\dot{L}$ $\displaystyle=$
$\displaystyle\left(1+d\,v^{2}\right)H\,L+c\,v\,,$ (71) $\displaystyle\dot{v}$
$\displaystyle=$
$\displaystyle\left(1-v^{2}\right)\left(\frac{k}{L}-d\,Hv\right)\quad{\rm(standard~{}walls)}$
(72)

where $H=\dot{a}/a$ is the expansion rate and $c$, $k$ are constants – the so-
called energy-loss and momentum parameters respectively. Comparing (71) with
(7) (with $n_{\rm K}=1/L$) it is clear that the energy loss parameter $c$ can
be interpreted as an effective cross section, for wall-wall interactions. For
power-law cosmologies, $a(t)\sim t^{\gamma}$, these equations lead to an
attractor scaling solution where $\dot{L}$ and $v$ equal a constant, in
agreement with field theory numerical simulations of the networks.

We can compare this to the equations that control the precursor wall networks.
We can obtain these equations by arguing as in Sec. II. First of all, note
that for kinks in $1+1$ dimensions Eqs. (7) and (8) take the following
suggestive form when written in terms of the correlation length ($L=1/n_{\rm
K}$ in $1+1$ dimensions),

$\displaystyle\dot{L}$ $\displaystyle=$ $\displaystyle\sigma v\,,$ (73)
$\displaystyle\dot{v}$ $\displaystyle=$
$\displaystyle-\kappa\,\frac{v^{2}}{L}\,.$ (74)

In $d+1$ dimensions, the correlation length in the DW network is identified as
$L\equiv 1/{\cal A}$ with ${\cal A}$ the physical area density. In terms of
this, the equations would be basically unchanged, except that the values of
$\sigma,\kappa$ might depend on $d$.

The extension of the results in Sections II and IV to an expanding universe is
beyond the scope of this work. Yet, in the language of the VOS model (73)-(74)
the extension seems to suggest itself. Neglecting relativistic corrections (as
we are interested in diffusive solutions, approaching $v\to 0$), the natural
expectation is

$\displaystyle\dot{L}$ $\displaystyle=$ $\displaystyle H\,L+\sigma v\,,$ (75)
$\displaystyle\dot{v}$ $\displaystyle=$
$\displaystyle-\kappa\frac{v^{2}}{L}-dHv\qquad{\rm(precursor~{}walls)}$ (76)

where $d$ is the number of space dimensions.

The main difference between (75)-(76) and (71)-(72) is in the tension term,
$k/L$, present in (72). This tends to increase $v$ in proportion to the DW
curvature, and results from the NG equation of motion [40, 42]. The term is
absent in the VOS for the precursor walls, which are subject to some
frictional force but no accelerating force from the tension.

Equations (75)-(76) have 2 solutions. The trivial one,

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

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)=t^{\gamma}$, is

$\displaystyle L(t)$ $\displaystyle\propto$ $\displaystyle t^{\mu}\,,\qquad
v(t)=v_{0}\,t^{\mu-1}\quad{\rm with}$ (78) $\displaystyle\mu$ $\displaystyle=$
$\displaystyle\frac{\sigma-(d\,\sigma-\kappa)\,\gamma}{\kappa+\sigma}\,,$ (79)
$\displaystyle v_{0}$ $\displaystyle=$
$\displaystyle\frac{1-(d+1)\,\gamma}{\kappa+\sigma}~{}.$ (80)

This is the generalization to an expanding universe (and to walls) of the
diffusive plasma of kinks of Sec. II. Recall that $v$ is a measure of the
magnitude of the typical velocity. The positivity of $v_{0}$ then implies that
(assuming $\kappa+\sigma>0$) this solution exists only for sufficiently slow
expansion $\gamma<\gamma_{max}\equiv 1/(d+1)$. At $\gamma_{max}$,
$\mu_{max}=\gamma_{max}$ and so this solution merges with the ‘DW gas’
$L(t)\propto a(t)$ solution. Radiation-domination and matter-domination
correspond to $\gamma=2/(d+1)$ and $\gamma=2/d$, too fast an expansion to
allow for the nontrivial diffusive scaling.

One must keep in mind, however, that the phenomenological parameters $\kappa$,
$\sigma$ can actually depend on the expansion rate, i.e. on $\gamma$ [43, 44].
As we argue in the next section, there is a reason to expect that the correct
scaling for precursor walls in an expanding universe retains the form

$L\sim t^{1/2}~{},$

(that is, $\mu=1/2$), in terms of the proper cosmic time $t$. We can translate
this into the following condition for $\kappa$,

$\kappa=\sigma\frac{1-2\,d\,p}{1-2\,p}~{},$ (81)

leading to $v_{0}=(1-2\,p)/\sigma$. This scaling exists and is distinct from
the DW gas only for $p<1/2$. This makes sense: the threshold separation
between diffusive and DW gas behaviours must be when the expansion $a$ is
slower/faster than the diffusion itself. Note that $\kappa$ is negative for
$1/2\,d<p<1/2$. In order to keep up with diffusion $\kappa$ needs to become an
accelerating term.

##  VI Discussion

We have presented a new analytic method to compute the annihilation of DW
networks driven by population bias (an asymmetric distribution of the nearly-
degenerate vacua). We have used standard QFT methods in flat spacetime to i)
identify the appropriate initial condition for the network from the quantum
field ground state at the transition, and ii) we have computed the time
evolution of the most important quantity that characterizes the DW network,
namely, the area density.

Our main result is the computation of the DW area per unit volume ${\cal
A}^{d+1}=({\rm DW\,area})/({\rm volume})$ in $d+1$ dimensions, depicted in
Figs. 3 and 4 for the $d=1$ case. (Qualitatively a similar behaviour occurs in
higher dimensions.) Asymptotically, ${\cal A}^{d+1}$ obeys the decay law

${\cal
A}^{d+1}\sim\frac{1}{\sqrt{t}}\,\exp\left[-(t/t_{ann})^{d/2}\right]~{}.$ (82)

Let us also note that the computation captures an intriguing transient
increase which apparently can be understood as nucleation.

Let us comment on a few points. First, note that our flat space result (85)
agrees with the condensed matter result [24]. The method used in [24] (which
later inspired [26]) appears to be radically different from ours but it shares
important ingredients. In [24] the walls are defined by an auxiliary field $u$
with a statistical averaging with Gaussian statistics, and they move according
to a given equation of motion. The equation of motion, however, is not the
Nambu Goto (NG) equation but the so-called Allen-Cahn equation [45], which is
relevant in friction dominated finite density systems, see [46] for a review.
It is a nonlinear dissipative equation which, ignoring nonlinearities, has the
structural form $\dot{u}=D\,\partial_{i}^{2}u$ with $D$ a constant and
$\partial_{i}^{2}$ the laplacian.

Our method builds upon a Gaussian relativistic quantum field, not from an
equation of motion for the walls. Instead, we track the precursor walls (the
‘zeros’ of the field $\phi-\delta\phi(t)$). Such walls are not expected to
obey a standard NG equation. Still, at late times the overall motion is
expected to obey non-relativistic scaling, simply because the motion of walls
that survive for a long time is encoded in long-wavelength modes of a massive
field. Indeed, the non-relativistic limit of a tachyonic massive field is
easily obtained by introducing the decomposition
$\phi(t,x)=e^{|m|t}\psi(t,x)$, leading to a diffusion-like equation
$|m|\dot{\psi}=\partial_{i}^{2}\psi$ at late times (when
$\ddot{\psi}\ll|m|\dot{\psi}$).

Similary, it is possible to obtain the equation of motion that precursor DWs
obey by treating them as semiclassical objects, that is, as DW-like solutions
to the classical tachyonic Klein Gordon equation. The analogue to the kink
solution is

$\phi=\exp{\big{(}\sqrt{k^{2}+|m^{2}|}\,t\big{)}}\sinh{\big{(}kz\big{)}}~{},$
(83)

with arbitrary $k$ and overall amplitude. Indeed, there is a zero in the $z=0$
plane and the exponential time dependence results from the field not being
stabilized yet. It is easy to see that the ansatz
$\exp{\big{(}\sqrt{k^{2}+|m^{2}|}\,t\big{)}}\sinh{\big{(}k[z-z_{0}(t,x,y)]\big{)}}$
is also a solution at linear order in the bending $z_{0}(t,x,y)$ provided it
satisfies222An improved ansatz exists giving an exact bent wall solution that
holds everywhere, and which leads to a nonlinear equation for $z_{0}$.
However, it doesn’t change much the qualitative properties of (84) so we do
not show it here.

$2\sqrt{k^{2}+|m^{2}|}\;\dot{z}_{0}=\partial_{i}^{2}\,z_{0}~{}.$ (84)

Again, the $\ddot{z}_{0}$ term has been neglected since we consider late
times/long wavelengths. Thus, the nonrelativistic limit makes these walls obey
an Allen-Cahn-like equation, and this explains the agreement with the
condensed matter result [24] and the origin of the ‘diffusive’ scaling $L\sim
t^{1/2}$ of precursor walls. (Notice that for scaling solutions the value of
the diffusion constant factors out.)

The extension of our results to an expanding FRW cosmology is beyond the scope
of this work, however, the reasoning above suggests some expectations.
Clearly, there are two regimes according to whether the scalar mass is
bigger/smaller than the Hubble parameter $H$. For $H\gg|m|$, the field motion
is frozen (even for $m^{2}<0$), and for $H\ll|m|$ we expect precursor walls
with non-relativistic limit of the form $\dot{z}_{0}\propto
a(t)^{-2}\partial_{i}^{2}\,z_{0}$ with $t$ the cosmic time and $\partial_{i}$
the comoving spatial gradient. Thus, we expect that without bias there might
be a scaling regime where the physical correlation length and velocity keep
the same behaviour, $L\sim t^{1/2}$ and $v\sim t^{-1/2}$ in terms of cosmic
time. This scaling differs from the ‘DW gas’ ($L\sim a(t)$, $v=0$) and it
might be realized if the expansion rate is slower than $a(t)\sim t^{1/2}$. The
especially interesting case of radiation domination in $3+1$ seems marginal.
In matter domination, $t^{2/3}$, the diffusive scaling $L\sim t^{1/2}$ is not
expected to be realized.

We can turn now to the network annihilation – to the exponential decay induced
by the bias. Both (82) and Hindmarsh’s result (4) suggest that the exponent is
proportional to the correlation volume, $L^{d}$ (see also [27]) – the comoving
volume in the cosmological case [26]. This suggests that for a precursor wall
network the bias would lead to a suppression of the form

$\exp\left[-{\rm const}\,\left(\frac{\eta}{a(\eta)}\right)^{d/2}\right]\,.$
(85)

This suppression is much milder than (4), more so for faster expansion rates.
In fact, for $a(t)$ faster than (or equal to) $t^{1/2}$ we would expect DW gas
behaviour, $L\sim a(t)$, and the exponent becomes time independent.

Let us emphasize that (85) (if confirmed) is not in contradiction with
Hindmarsh’s result (4). They simply refer to different regimes: (4) holds for
stabilized DWs that obey the NG equation whereas (85) would apply for
precursor walls, that are similar to the DW gas limit. As mentioned above, a
physical precursor wall regime is feasible in cosmology. This should happen
for instance if the Hubble rate $H$ at the symmetry breaking transition is
significantly larger than the typical mass scale near the symmetry breaking
transition. (Another possibility is that the symmetry breaking field is light
during inflation.) This can result in the formation of an underdense DW
network (less than one DW per Hubble patch). Even for stabilized DWs the
initial evolution should be close to the DW gas regime $L\sim a(t)$. As in the
DW gas limit of precursor walls, a population bias is then also expected not
to lead to exponential suppression in time.

It seems, then, that different decay laws, (3), (4), (85) (or even different
ones), can apply for different network models/realizations. It is clear from
the previous paragraph that a crucial ingredient that specifies the DW network
is its initial condition, say, at the symmetry breaking transition.

After all, the network annihilation process consists in the collapse of (fewer
and fewer) closed DW structures larger in size than the horizon at the
annihilation time $t_{ann}$. By definition, then, these are encoded in super-
horizon modes, and they remain basically frozen during the scaling period.
When they evolve, they do so once the scaling period is over. So, it seems
rather plausible that the statistics of how many DWs are present after the
typical annihilation time is also considerably affected by the initial
condition, for deeply super-horizon modes.

It is relevant to compare the assumptions on the initial state in some
previous works. In Ref. [26] a scale invariant (white noise) initial power is
assumed. The numerical simulations in [28, 29] set the field at the two vacua
with (biased) random probability at each lattice sites It isn’t entirely
obvious how to map from one to the other. It seems possible that this is why
(4) fails for the population bias simulations in $2+1$ of [28, 29]. In our
treatment, the initial condition is encoded in the QFT vacuum of the free
massive field. This assumption accounts for the formation of the DWs at the
spontaneous symmetry breaking transition.

The present work can be extended in several directions. Whether the guess for
the decay (85) really applies in an expanding space, even for precursor walls,
requires confirmation. It is possible in principle to include different
initial states (for instance a thermal state at formation), as well as a
different time dependence for both the explicit and spontaneous breaking of
the discrete symmetry. We have considered the simplest model that gives rise
to a DW network, with a $\lambda\phi^{4}$ double-well potential. In principle
the method described in this work to count the DW area density can be extended
to other models, an interesting target being axionic models. We leave these
questions for future work.

###### Acknowledgements.

We thank R. Ferreira, F. Rompineve, T. Vachaspati and M. Mukhopadhyay for
useful discussions and comments. This work is supported by projects
PID2020-115845GB-I00/AEI/10.13039/501100011033 and 2017-SGR-1069. IFAE is
partially funded by the CERCA program of the Generalitat de Catalunya. The
work of GZ was supported by a fellowship from “La Caixa” Foundation (ID
100010434) and from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 847648. The
fellowship code is LCF/BQ/PI20/11760021.

## References

  * Vilenkin and Shellard [2000] A. Vilenkin and E. P. S. Shellard,  _Cosmic Strings and Other Topological Defects_ (Cambridge University Press, 2000).
  * Vachaspati [2010] T. Vachaspati,  _Kinks and domain walls: An introduction to classical and quantum solitons_ (Cambridge University Press, 2010).
  * Zeldovich  _et al._ [1974] Y. B. Zeldovich, I. Y. Kobzarev, and L. B. Okun, Zh. Eksp. Teor. Fiz. 67, 3 (1974).
  * Lazanu  _et al._ [2015] A. Lazanu, C. J. A. P. Martins, and E. P. S. Shellard, Phys. Lett. B 747, 426 (2015), arXiv:1505.03673 [astro-ph.CO] .
  * Hiramatsu  _et al._ [2010] T. Hiramatsu, M. Kawasaki, and K. Saikawa, JCAP 05, 032 (2010), arXiv:1002.1555 [astro-ph.CO] .
  * Hiramatsu  _et al._ [2011] T. Hiramatsu, M. Kawasaki, and K. Saikawa, JCAP 08, 030 (2011), arXiv:1012.4558 [astro-ph.CO] .
  * Kawasaki and Saikawa [2011] M. Kawasaki and K. Saikawa, JCAP 09, 008 (2011), arXiv:1102.5628 [astro-ph.CO] .
  * Hiramatsu  _et al._ [2014] T. Hiramatsu, M. Kawasaki, and K. Saikawa, JCAP 02, 031 (2014), arXiv:1309.5001 [astro-ph.CO] .
  * Sousa and Avelino [2015] L. Sousa and P. P. Avelino, Phys. Rev. D 92, 083520 (2015), arXiv:1507.01064 [astro-ph.CO] .
  * Krajewski  _et al._ [2016] T. Krajewski, Z. Lalak, M. Lewicki, and P. Olszewski, JCAP 12, 036 (2016), arXiv:1608.05719 [astro-ph.CO] .
  * Nakayama  _et al._ [2017] K. Nakayama, F. Takahashi, and N. Yokozaki, Phys. Lett. B 770, 500 (2017), arXiv:1612.08327 [hep-ph] .
  * Krajewski  _et al._ [2018] T. Krajewski, Z. Lalak, M. Lewicki, and P. Olszewski, JCAP 05, 007 (2018), arXiv:1709.10100 [hep-ph] .
  * Saikawa [2017] K. Saikawa, Universe 3, 40 (2017), arXiv:1703.02576 [hep-ph] .
  * Gelmini  _et al._ [2021] G. B. Gelmini, S. Pascoli, E. Vitagliano, and Y.-L. Zhou, JCAP 02, 032 (2021), arXiv:2009.01903 [hep-ph] .
  * Takahashi and Yin [2021] F. Takahashi and W. Yin, JCAP 04, 007 (2021), arXiv:2012.11576 [hep-ph] .
  * Babichev  _et al._ [2022] E. Babichev, D. Gorbunov, S. Ramazanov, and A. Vikman, JCAP 04, 028 (2022), arXiv:2112.12608 [hep-ph] .
  * Dunsky  _et al._ [2022] D. I. Dunsky, A. Ghoshal, H. Murayama, Y. Sakakihara, and G. White, Phys. Rev. D 106, 075030 (2022), arXiv:2111.08750 [hep-ph] .
  * Zambujal Ferreira  _et al._ [2022] R. Zambujal Ferreira, A. Notari, O. Pujolàs, and F. Rompineve, Phys. Rev. Lett. 128, 141101 (2022), arXiv:2107.07542 [hep-ph] .
  * Ramberg  _et al._ [2022] N. Ramberg, W. Ratzinger, and P. Schwaller, (2022), arXiv:2209.14313 [hep-ph] .
  * Gelmini  _et al._ [2022] G. B. Gelmini, A. Simpson, and E. Vitagliano, (2022), arXiv:2207.07126 [hep-ph] .
  * Ferreira  _et al._ [2022] R. Z. Ferreira, A. Notari, O. Pujolas, and F. Rompineve, (2022), arXiv:2204.04228 [astro-ph.CO] .
  * Jiang and Huang [2022] Y. Jiang and Q.-G. Huang, (2022), arXiv:2208.00697 [astro-ph.CO] .
  * Blasi and Mariotti [2022] S. Blasi and A. Mariotti, (2022), arXiv:2203.16450 [hep-ph] .
  * Ohta  _et al._ [1982] T. Ohta, D. Jasnow, and K. Kawasaki, Phys. Rev. Lett. 49, 1223 (1982).
  * Coulson  _et al._ [1996] D. Coulson, Z. Lalak, and B. A. Ovrut, Phys. Rev. D 53, 4237 (1996).
  * Hindmarsh [1996] M. Hindmarsh, Phys. Rev. Lett. 77, 4495 (1996), arXiv:hep-ph/9605332 .
  * Larsson  _et al._ [1997] S. E. Larsson, S. Sarkar, and P. L. White, Phys. Rev. D 55, 5129 (1997), arXiv:hep-ph/9608319 .
  * Correia  _et al._ [2014] J. R. C. C. C. Correia, I. S. C. R. Leite, and C. J. A. P. Martins, 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, Phys. Rev. D 97, 083521 (2018), arXiv:1804.10761 [astro-ph.CO] .
  * Krajewski  _et al._ [2021] T. Krajewski, J. H. Kwapisz, Z. Lalak, and M. Lewicki, Phys. Rev. D 104, 123522 (2021), arXiv:2103.03225 [astro-ph.CO] .
  * Stauffer [1979] D. Stauffer, Phys. Rept. 54, 1 (1979).
  * Lalak  _et al._ [1995] Z. Lalak, B. A. Ovrut, and S. Thomas, Phys. Rev. D 51, 5456 (1995).
  * Ferrer  _et al._ [2019] F. Ferrer, E. Masso, G. Panico, O. Pujolas, and F. Rompineve, Phys. Rev. Lett. 122, 101301 (2019), arXiv:1807.01707 [hep-ph] .
  * Garriga  _et al._ [2016] J. Garriga, A. Vilenkin, and J. Zhang, JCAP 02, 064 (2016), arXiv:1512.01819 [hep-th] .
  * Deng  _et al._ [2017] H. Deng, J. Garriga, and A. Vilenkin, JCAP 04, 050 (2017), arXiv:1612.03753 [gr-qc] .
  * Mukhopadhyay  _et al._ [2020a] M. Mukhopadhyay, T. Vachaspati, and G. Zahariade, Phys. Rev. D 102, 056021 (2020a), arXiv:2004.07249 [hep-th] .
  * Mukhopadhyay  _et al._ [2020b] M. Mukhopadhyay, T. Vachaspati, and G. Zahariade, Phys. Rev. D 102, 116002 (2020b), arXiv:2009.11480 [hep-th] .
  * Ben-Naim  _et al._ [1993] E. Ben-Naim, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 70, 1890 (1993).
  * Piasecki [1995] J. Piasecki, Phys. Rev. E 51, 5535 (1995).
  * Avelino  _et al._ [2005] P. P. Avelino, C. J. A. P. Martins, and J. C. R. E. Oliveira, Phys. Rev. D 72, 083506 (2005), arXiv:hep-ph/0507272 .
  * Bishop [2013] C. Bishop,  _Pattern Recognition and Machine Learning: All ”just the Facts 101” Material_, Information science and statistics (Springer (India) Private Limited, 2013).
  * Kawano [1990] L. Kawano, Phys. Rev. D 41, 1013 (1990).
  * Martins  _et al._ [2016] C. J. A. P. Martins, I. Y. Rybak, A. Avgoustidis, and E. P. S. Shellard, Phys. Rev. D 93, 043534 (2016), arXiv:1602.01322 [hep-ph] .
  * Avelino [2020] P. P. Avelino, Phys. Rev. D 101, 023514 (2020), arXiv:1910.07011 [astro-ph.CO] .
  * Allen and Cahn [1979] S. M. Allen and J. W. Cahn, Acta Metallurgica 27, 1085 (1979).
  * Bray [1994] A. Bray, Advances in Physics 43, 357 (1994), https://doi.org/10.1080/00018739400101505 .