#  Hubble-induced phase transitions on the lattice with applications to Ricci
reheating

Dario Bettoni  Asier Lopez-Eiguren  Javier Rubio

###### Abstract

Using 3+1 classical lattice simulations, we follow the symmetry breaking
pattern and subsequent non-linear evolution of a spectator field non-minimally
coupled to gravity when the post-inflationary dynamics is given in terms of a
stiff equation-of-state parameter. We find that the gradient energy density
immediately after the transition represents a non-negligible fraction of the
total energy budget, steadily growing to equal the kinetic counterpart. This
behaviour is reflected on the evolution of the associated equation-of-state
parameter, which approaches a universal value $1/3$, independently of the
shape of non-linear interactions. Combined with kination, this observation
allows for the generic onset of radiation domination for arbitrary self-
interacting potentials, significantly extending previous results in the
literature. The produced spectrum at that time is, however, non-thermal,
precluding the naive extraction of thermodynamical quantities like
temperature. Potential identifications of the spectator field with the
Standard Model Higgs are also discussed.

##  1 Introduction

The appearance of tachyonic instabilities is an ubiquitous phenomenon in the
physics of the very early Universe, being generically associated with the
amplification of quantum fluctuations. Typical examples include natural and
new inflation models [1, 2, 3, 4], preheating scenarios with trilinear
interactions [5, 6] or spontaneous symmetry breaking settings such as hybrid
inflation [7, 8, 9, 10].

On general grounds, a tachyonic instability emerges whenever a mass eigenvalue
in an interacting or self-interacting theory becomes negative. This transition
can be caused by a bare mass parameter in the action or effectively induced by
the expectation value of a given field or operator. In this context, models in
which the instability is triggered by a non-minimal coupling to gravity have
recently received considerable attention.

Hubble-induced phase transitions for spectator fields have been advocated to
address several phenomenological aspects such as the generation of the
observed matter-antimatter asymmetry [11], the formation of topological
defects with observable gravitational wave backgrounds [12, 13], the onset of
the hot big bang era [14, 15, 16] or the non-thermal production of dark matter
[17, 18]. The main assumption behind these proposals is the existence of a
cosmological expansion era with a stiff equation-of-state parameter ($w>1/3$)
prior to radiation domination.111For a review of the different expansion
histories in the early Universe, see for instance Ref. [19].

In many instances, the evolution of the spectator field during Hubble-induced
phase transitions is discussed in terms of an homogeneous component rolling
down an effective potential and eventually oscillating around the new
symmetry-breaking minima in a fully coherent fashion [14, 15, 16]. This naive
picture is, however, inaccurate [7, 8]. On the one hand, it fails to capture
one of the main features of spontaneous symmetry breaking: the production of
topological defects [12]. On the other hand, the dynamics of the system is not
dominated by an homogeneous component, but rather by the rapid amplification
of long wavelength fluctuations that become relevant well before the field
reaches the new minima of the potential [13]. Among other effects, this has
important consequences for the energy equipartition and thermalization of the
system.

In order to verify and assess the relevance of the aforementioned features, we
explore here, for the first time, the non-linear dynamics of Hubble-induced
phase transitions using $3+1$ classical lattice simulations. To this end, we
focus on a $Z_{2}$ symmetric spectator field non-minimally coupled to gravity,
albeit our results can be easily extended to other global symmetry breaking
patterns.

The present analysis differs in various aspects from previous studies in the
context of hybrid inflation, for which lattice simulations are also available
[7, 8, 9, 10]. First, the tachyonic mass leading to the spontaneous breaking
of the aforementioned $Z_{2}$ symmetry is time-dependent. Second, the system
is completely characterized by a single physical scale, the Hubble rate, which
controls both the dilution and the tachyonic enhancement of perturbations.
When numerical tools are brought to bear, we find also several interesting
outcomes:

  1. i)

The Hubble-induced tachyonic instability leads to the formation of regions
with positive and negative field values separated by domain wall
configurations. These domains melt periodically at every field oscillation to
reappear again just slightly distorted.

  2. ii)

The time-dependent character of the symmetry-breaking minima prevents the
usual stabilization around them, giving rise to a bi-modal distribution
pattern that survives for a large number of oscillations.

  3. iii)

Although the energy transfer to gradients is initially just a small fraction
of the total energy density, it grows monotonically at the expense of the
interaction and potential contributions. This leads to a final virialized
state in which the kinetic energy density of the spectator field equals its
gradient counterpart.

  4. iv)

The virialization process is associated with a smooth turbulent regime,
similar to that in Refs. [20, 21]. The continuous population of ultraviolet
(UV) modes by this mechanism translates into an effective equation of state
$w_{\chi}\simeq 1/3$. Combined with kination, this allows for the onset of
radiation domination for general interacting potentials. This contrasts with
the naive expectation following from the homogeneous approximation, where this
transition is only possible for a quartic interaction [16].

  5. v)

Even if the late-time evolution of the spectator coincides with that of a
radiation fluid, the associated spectrum of perturbations remains highly non-
thermal for an extended period of time, precluding the naive extraction of
thermodynamical quantities like temperature.

The paper is organized as follows. In Section 2, we introduce the basic
ingredients of the model, leaving for Section 3 the description of the rapid
amplification of the spectator field fluctuations through the Hubble-induced
tachyonic instability. This classicalization process is essential for the
analysis of the subsequent non-linear stage, which we study in Section 4 via
$3+1$ classical lattice simulations. The implications of our findings for the
onset of radiation domination are discussed in Section 5. Section 6 contains a
summary of our main theoretical assumptions and future research lines,
including a discussion on the possibility of identifying the spectator field
with the Standard Model Higgs. Finally, our conclusions are presented in
Section 7.

##  2 Non-minimally coupled spectator field

A cosmological expansion epoch with a stiff equation-of state parameter
($w>1/3$) appears naturally in non-oscillatory quintessential inflation models
aiming to provide a unified description of inflation and dark energy using a
single degree of freedom [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36]. In this type of scenarios, the absence of a potential minimum
triggers generically a kination or deflation period immediately after the end
of inflation. During this era, the total energy budget of the Universe is
dominated by the kinetic energy density of the inflaton field, such that the
global effective equation of state approaches $w\simeq 1$. Within this
context, let us consider the following $Z_{2}$ symmetric action for an
energetically subdominant spectator field $\chi$,

$S_{\chi}=\int
d^{4}x\sqrt{-g}\left[-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\chi\partial_{\mu}\chi-\frac{1}{2}\xi
R\chi^{2}-\frac{1}{4}\lambda\chi^{4}\right]\,,$ (2.1)

where we have intentionally omitted a potential mass contribution which, if
sufficiently small, will not play any essential role in the following
discussions [11, 12, 13]. Here $R$ stands for the Ricci scalar constructed out
of the metric $g_{\mu\nu}$ and $\xi$ and $\lambda$ are two dimensionless
coupling constants that we will assume to be positive-definite in what
follows. The dynamics of the field $\chi$ is characterized by the Klein–Gordon
equation,

$\Box\,\chi-\xi R\chi-\lambda\chi^{3}=0\,,$ (2.2)

together with the energy density and pressure following from the $00$ and $ii$
components of the associated energy-momentum tensor,

$T^{(\chi)}_{\mu\nu}=D_{\mu}\chi D_{\nu}\chi-
g_{\mu\nu}\left[\frac{1}{2}D^{\rho}\chi
D_{\rho}\chi+\frac{\lambda}{4}\,\chi^{4}\right]+\xi\Big{[}G_{\mu\nu}+\left(g_{\mu\nu}\Box-
D_{\mu}\,D_{\nu}\right)\Big{]}\chi^{2},$ (2.3)

with $D_{\mu}$ the covariant derivative, $G_{\mu\nu}$ the Einstein tensor and
$\Box=D^{\lambda}D_{\lambda}$ the d’Alambertian operator. Note that the
contraction of $T^{(\chi)}_{\mu\nu}$ with arbitrary timelike vectors $u^{\mu}$
and $u^{\nu}$ is not guaranteed to be positive definite at all times [37, 38,
39, 40], thus violating the weak energy condition
$T^{(\chi)}_{\mu\nu}u^{\mu}u^{\nu}\geq 0$ within this particular sector.

For a fixed flat Friedmann–Lemaître–Robertson–Walker metric $g_{\mu\nu}={\rm
diag}(-1,a^{2}(t)\,\delta_{ij})$ with scale factor $a(t)$, we have

$R=3(1-3w)H^{2}\,,\hskip 19.91692ptG_{00}=3H^{2}\,,\hskip
19.91692ptG_{ij}=3\,w\,H^{2}g_{ij}\,,\hskip
19.91692pt\Box\chi=-\ddot{\chi}-3H\dot{\chi}+\frac{\nabla^{2}}{a^{2}}\,,$
(2.4)

with $H=\dot{a}/a$ the Hubble rate, the dots denoting derivatives with respect
to the coordinate time $t$, and $w$ the global equation–of–state parameter.
Using these expressions, we can rewrite Eqs. (2.2) and (2.3) as

$\displaystyle\ddot{\chi}+3H\dot{\chi}-a^{-2}\nabla^{2}\chi+3\xi(1-3w)H^{2}\chi+\lambda\,\chi^{3}=0\,,$
(2.5)
$\displaystyle\rho_{\chi}=T_{00}=\frac{1}{2}\dot{\chi}^{2}+\frac{1}{2a^{2}}|\nabla\chi|^{2}+\frac{\lambda}{4}\,\chi^{4}+3\xi\left[H^{2}+H\partial_{t}-\frac{1}{3}\frac{\nabla^{2}}{a^{2}}\right]\chi^{2}\,,$
(2.6) $\displaystyle
p_{\chi}=\frac{1}{3a^{2}}\sum_{i}T_{ii}=\frac{1}{2}\dot{\chi}^{2}-\frac{1}{6a^{2}}|\nabla\chi|^{2}-\frac{\lambda}{4}\,\chi^{4}+3\xi\left[wH^{2}-\frac{1}{3}\left(\partial^{2}_{t}+2H\partial_{t}\right)+\frac{2}{9}\frac{\nabla^{2}}{a^{2}}\right]\chi^{2}\,.$
(2.7)

Since we are assuming the spectator field to be energetically subdominant, the
equation-of-state parameter in these expressions can be well-approximated by
that of the quintessential inflation field. Having this in mind, we can
distinguish two regimes:

  1. 1.

Inflation ($w=-1$): During this epoch, the Hubble-induced mass term in Eq.
(2.5) is positive definite, making the spectator field heavy for sufficiently
large values of the non-minimal coupling ($\xi\gtrsim 1/12$). This confines
the $\chi$ field to the origin of its effective potential and suppresses the
generation of isocurvature perturbations [11, 12].

  2. 2.

Kination ($w=1$): The onset of a kinetic dominated era soon after the end of
inflation triggers the spontaneous symmetry breaking of the $Z_{2}$ symmetry
and the subsequent evolution of $\chi$ towards large field values, $|\chi_{\rm
min}|=\left(6\xi/\lambda\right)^{1/2}H$.

To study the dynamics of the spectator field in between these cosmological
epochs, we will assume the inflation-to-kination transition to be
instantaneous as compared with the evolution of the scalar field and
parametrize the temporal dependence of the scale factor and the Hubble rate as

$a=a_{\rm kin}\left[1+3H_{\rm kin}\left(t-t_{\rm
kin}\right)\right]^{1/3}\,,\hskip 56.9055ptH=\frac{H_{\rm kin}}{1+3H_{\rm
kin}(t-t_{\rm kin})}\,,$ (2.8)

with $a_{\rm kin}$ and $H_{\rm kin}$ the values of these quantities at the
onset of kinetic domination, arbitrarily defined at $t=t_{\rm kin}$.
Performing now the change of variables

$Y\equiv\frac{a}{a_{\rm kin}}\frac{\chi}{\,\chi_{*}}\,,\hskip
42.67912pt\vec{y}\equiv\,a_{\rm kin}\chi_{*}\,\vec{x}\,,\hskip
42.67912ptz\equiv\,a_{\rm kin}\chi_{*}\,\tau\,,$ (2.9)

with $\chi_{*}\equiv\sqrt{6\xi}H_{\rm kin}$ and conformal time $\tau\equiv\int
dt/a$, the spectator field action (2.1) becomes

$S_{\chi}=\int
d^{3}{\vec{y}}\,dz\,\left[\frac{1}{2}(Y^{\prime})^{2}-\frac{1}{2}|\nabla
Y|^{2}\,+\,\frac{1}{2}M^{2}(z)Y^{2}-\frac{\lambda}{4}\,Y^{4}\right]\,,$ (2.10)

with the primes denoting derivatives with respect to the dimensionless
conformal time $z$ and the spatial derivatives understood now as taken with
respect to the new spatial coordinates $\vec{y}$. All the explicit temporal
dependence in (2.10) is included in the effective mass term

$M^{2}(z)\equiv(1-6\xi)(\mathcal{H}^{2}+\mathcal{H}^{\prime})=(4\nu^{2}-1)\mathcal{H}^{2}\,,$
(2.11)

constructed out of the comoving Hubble rate

$\mathcal{H}(z)\equiv\frac{a^{\prime}(z)}{a(z)}=\frac{1}{2(z+\nu)}\,,\hskip
42.67912pt\mathcal{H}^{\prime}(z)=-2\mathcal{H}^{2}(z)\,,$ (2.12)

with

$\nu\equiv\sqrt{\frac{3\xi}{2}}$ (2.13)

an integration constant ensuring that $z(t_{\rm kin})=0$. In terms of the new
variables, the time-dependent minima of the potential are located at

$|Y_{\rm min}|=\left(\frac{4\nu^{2}-1}{\lambda}\right)^{1/2}\mathcal{H}(z)\,,$
(2.14)

while Eqs. (2.5) and (2.6) become

$Y^{\prime\prime}-\nabla^{2}\,Y-M^{2}(z)Y+\lambda\,Y^{3}=0\,,\hskip
42.67912pt\rho_{\chi}=\mathcal{X}\,\rho_{Y}\,,\hskip
42.67912ptp_{\chi}=\mathcal{X}\,p_{Y}\,,$ (2.15)

with $\mathcal{X}\equiv\chi_{*}^{4}(a_{\rm kin}/a)^{4}$ an overall rescaling
factor. Here

$\rho_{Y}\equiv K+G+V+I_{1}\,,\hskip 42.67912ptp_{Y}\equiv
K-\frac{1}{3}G-V+I_{2}\,,$ (2.16)

stand respectively for the dimensionless energy density and pressure of the
$Y$-field with

$K=\frac{1}{2}\left(Y^{\prime}-\mathcal{H}Y\right)^{2}\,,\hskip
28.45274ptG=\frac{1}{2}|\nabla Y|^{2}\,,\hskip
28.45274ptV=\frac{\lambda}{4}Y^{4},\\\ $ (2.17)
$I_{1}=2\nu^{2}\left[\mathcal{H}(Y^{2})^{\prime}-\mathcal{H}^{2}Y^{2}-\frac{1}{3}\nabla^{2}Y^{2}\right]\,,\hskip
14.22636ptI_{2}=2\nu^{2}\left[\mathcal{H}(Y^{2})^{\prime}-\mathcal{H}^{2}Y^{2}-\frac{1}{3}(Y^{2})^{\prime\prime}+\frac{2}{9}\nabla^{2}Y^{2}\right]\,,$
(2.18)

the kinetic ($K$), gradient ($G$), potential ($V$) and interaction
contributions ($I_{1}$,$I_{2}$) following from the corresponding terms in Eqs.
(2.6) and (2.7). As anticipated after Eq. (2.3), the spectator energy density
is not positive definite for non-vanishing $\nu$. Note, however, that the
total energy density including the inflaton field and other potential matter
components is positive definite at all times, as required by the $G_{00}$
component in (2.4).

##  3 Linear evolution of the quantum system

Immediately after symmetry breaking, the action (2.10) can be well
approximated by that of a free scalar field with time dependent mass $M(z)$,
making the problem exactly solvable. In this section, we follow the canonical
formalism of Ref. [10], which has been applied to the present scenario in Ref.
[13], where we refer the reader for a more detailed discussion. Expanding the
$Y$ field in Fourier modes with wavenumber

$\kappa\equiv|\vec{\kappa}|\equiv\frac{|\vec{k}|}{a_{\rm kin}\chi_{*}}\,,$
(3.1)

and introducing position and momentum operators satisfying the usual equal-
time commutation relations
$\left[Y_{\vec{\kappa}}(z),\Pi_{\vec{\kappa}^{\prime}}(z)\right]=i\delta^{3}(\vec{\kappa}+\vec{\kappa}^{\prime})$,
the quantum Hamiltonian following from Eq. (2.10) in the Gaussian
approximation becomes the sum of an infinite set of harmonic oscillators,

$H=\frac{1}{2}\int
d^{3}\kappa\,\left[\Pi_{\vec{\kappa}}(z)\,\Pi^{\dagger}_{\vec{\kappa}}(z)+\omega_{\kappa}^{2}(z)\,Y_{\vec{\kappa}}(z)\,Y^{\dagger}_{\vec{\kappa}}(z)\right]\,,$
(3.2)

with time-dependent dispersion relation

$\omega_{\kappa}^{2}(z)\equiv\kappa^{2}-M(z)^{2}\,.$ (3.3)

As clearly illustrated by this expression, the tachyonic instability induced
by kinetic domination affects not only the zero mode of the scalar field but
also all subhorizon modes below the time-dependent mass $M^{2}(z)$, i.e.,
those within a sliding window $\kappa_{\rm
min}(z)\lesssim\kappa\lesssim\kappa_{\rm max}(z)$, with

$\kappa_{\rm min}(z)=\mathcal{H}(z)\quad\quad{\rm and}\quad\quad\kappa_{\rm
max}=(4\nu^{2}-1)^{1/2}\kappa_{\rm min}\,.$ (3.4)

The evolution of this Gaussian quantum system is determined by the Heisenberg
equations,

$\frac{d}{dz}\begin{pmatrix}\Pi_{\vec{\kappa}}(z)\vspace{2mm}\cr
Y_{\vec{\kappa}}(z)\end{pmatrix}=\begin{pmatrix}0&-\omega^{2}_{\kappa}(z)\vspace{2mm}\cr
1&0\end{pmatrix}\begin{pmatrix}\Pi_{\vec{\kappa}}(z)\vspace{2mm}\cr
Y_{\vec{\kappa}}(z)\end{pmatrix}\,,$ (3.5)

being all physical information encoded in the two-mode correlation functions
at equal or different times,

$\langle
v^{I}_{\vec{\kappa}}(z)\,v^{J}_{\vec{\kappa}^{\prime}}(z^{\prime})\rangle=\Sigma^{IJ}_{\kappa}(z,z^{\prime})\,\delta^{3}({\vec{\kappa}}+{\vec{\kappa}^{\prime}})\,,$
(3.6)

with
$v_{\vec{\kappa}}\equiv\left(\Pi_{\vec{\kappa}}(z),Y_{\vec{\kappa}}(z)\right)^{T}$
and the brackets indicating a quantum average over the initial state. Using
Eq. (3.5), we can express the correlation matrix
$\Sigma^{IJ}_{\kappa}(z,z^{\prime})$ in (3.6) in terms of its value at any
reference time $z_{r}$, namely [41, 42, 43, 10]

$\Sigma_{\kappa}(z,z^{\prime})={\mathbf{M}}_{\kappa}(z)\,\Sigma_{\kappa}(z_{r},z_{r})\,{\mathbf{M}}^{T}_{\kappa}(z^{\prime})\,,\hskip
14.22636pt{\mathbf{M}}_{\kappa}(z)=\begin{pmatrix}\sqrt{\frac{2}{\kappa}}\,\textrm{Re}\,g_{\kappa}(z)&\sqrt{2\kappa}\,\textrm{Im}\,g_{\kappa}(z)\vspace{2mm}\cr-\sqrt{\frac{2}{\kappa}}\,\textrm{Im}\,f_{\kappa}(z)&\sqrt{2\kappa}\,\textrm{Re}\,f_{\kappa}(z)\end{pmatrix}\,,$
(3.7)

with $f_{\kappa}$ a solution of the Schrödinger-like differential equation

$f_{\kappa}^{\prime\prime}+\omega^{2}_{\kappa}(z)\,f_{\kappa}=0\,,$ (3.8)

with initial conditions $f_{\kappa}(z_{r})$ and $f^{\prime}_{\kappa}(z_{r})$
and $g_{\kappa}\equiv if^{\prime}_{\kappa}$. Restricting ourselves to the
large $\nu$ limit (cf. Section 6) and assuming the system to be initially in a
vacuum state with $f_{\kappa}(z_{r}=0)=1/\sqrt{2\kappa}$ and
$f^{\prime}_{\kappa}(z_{r}=0)=-i\sqrt{\kappa/2}$, the solution of Eq. (3.8)
takes the compact form [13]

$f_{\kappa}(z)=\sqrt{z+\nu}\left[A_{\kappa}\,\mathcal{J}_{\nu}(\kappa(z+\nu))-B_{\kappa}\,\mathcal{Y}_{\nu}(\kappa(z+\nu))\right]\,,$
(3.9)

with $\mathcal{J_{\nu}}$ and $\mathcal{Y_{\nu}}$ the Bessel’s functions of the
first and second kind,

$A_{\kappa}={\cal Y}_{\nu}(\kappa\nu)\,\delta\,,\hskip
42.67912ptB_{\kappa}=\,{\cal
J}_{\nu}(\kappa\nu)\,\delta-\frac{f_{\kappa}(0)}{\sqrt{\nu}{\cal
Y}_{\nu}(\kappa\nu)}\,,$ (3.10)

and

$\delta=\frac{\pi
f_{\kappa}(0)}{4\sqrt{\nu}}\left[1-2\nu+2\nu\left(\kappa\frac{{\cal
Y}_{\nu-1}(\kappa\nu)}{{\cal
Y}_{\nu}(\kappa\nu)}-\frac{f^{\prime}_{\kappa}(0)}{f_{\kappa}(0)}\right)\right]\,.$
(3.11)

Using these expressions, we can easily evaluate the equal-time correlation
matrix

$\Sigma_{\kappa}(z,z)=\begin{pmatrix}|g_{\kappa}(z)|^{2}&F_{\kappa}(z)-\frac{i}{2}\vspace{3mm}\cr
F_{\kappa}(z)+\frac{i}{2}&|f_{k}(z)|^{2}\end{pmatrix}\,,$ (3.12)

with

$F_{\kappa}(z)\equiv\frac{1}{2}\,\langle\Pi_{\vec{\kappa}}(z)\,Y_{\vec{\kappa}}^{\dagger}(z)+Y_{\vec{\kappa}}(z)\,\Pi_{\vec{\kappa}}^{\dagger}(z)\rangle={\rm
Im}\,(f_{\kappa}^{*}g_{\kappa})$ (3.13)

a WKB phase directly related to the Heisenberg uncertainty principle [41, 43,
42, 10]

$\Delta
Y_{\kappa}^{2}\,\Delta\Pi_{\kappa}^{2}=|F_{\kappa}(z)|^{2}+\frac{1}{4}\geq{\frac{1}{4}}\,\Big{|}\langle[Y_{\kappa}(z),\
\Pi^{\dagger}_{\kappa}(z)]\rangle\Big{|}^{2}\,.$ (3.14)

For $|F_{\kappa}(z)|\gg 1$, the commutator $\langle\\{Y_{\kappa}(z),\
\Pi_{\kappa}^{\dagger}(z)\\}\rangle$ in this expression exceeds significantly
the anticommutator
$\langle|[Y_{\kappa}(z),\Pi_{\kappa}^{\dagger}(z)]|\rangle$,

$\langle\\{Y_{\kappa}(z),\
\Pi_{\kappa}^{\dagger}(z)\\}\rangle\gg\langle|[Y_{\kappa}(z),\Pi_{\kappa}^{\dagger}(z)]|\rangle=\hbar\,,$
(3.15)

meaning that the properties of the system can be well-approximated by those of
a classical random field (see Fig. 1).

The equal-time correlation matrix (3.12) allows the computation of several
related quantities characterizing the Gaussian random field [13]. Among those
to be considered in the subsequent analysis, we can highlight the two-point
correlation function obtained by Fourier transforming the diagonal entry
$|f_{\kappa}(z)|^{2}$,

$\zeta(\vec{y},z)\equiv\left\langle
Y(\vec{y},z)Y(0,z)\right\rangle=\int\frac{d^{3}\kappa}{(2\pi)^{3}}\,e^{i\vec{\kappa}\cdot\vec{y}}\,|f_{\kappa}(z)|^{2}\,.$
(3.16)

For a highly-symmetric spherical configuration and large $\nu$ values, the
momentum integral in this expression can be computed explicitly, getting [13]

$\zeta(r,z)\simeq\zeta(0,z)\,G_{1}(\kappa_{*}r)\,,$ (3.17)

with

$\zeta(0,z)\equiv Y^{2}_{\rm
rms}(z)=\left(\frac{2\nu-1}{8\pi\nu}\right)^{2}\left(1+\frac{z}{\nu}\right)^{2\nu+1}\kappa_{*}^{2}(z)$
(3.18)

the square of the time-dependent dispersion determining the root mean-square
perturbation (rms), $\kappa_{*}(z)\equiv 2\sqrt{\nu+1}\,\kappa_{\rm min}(z)$ a
typical momentum scale and

$G_{1}(\kappa_{*}r)\equiv\sqrt{\frac{\pi}{2}}\frac{1}{\kappa_{*}r}\exp\left(-\frac{1}{2}\kappa^{2}_{*}r^{2}\right)\text{erfi}\left(\frac{\kappa_{*}r}{\sqrt{2}}\right)\,,$
(3.19)

a shape function, with erfi the imaginary error function [44].

##  4 Non-linear evolution of the classicalized system

As we have seen in the previous section, the quantum evolution of the
spectator field drives the low momenta modes in Fourier space into a highly
occupied state where quantum averages can be well approximated by classical
ones [10, 13]. This allows us to study the subsequent non-linear evolution
using just classical lattice simulations, instead of more involved
$2$-particle irreducible techniques [45, 46, 47].

###  4.1 Lattice discretization

Our approach relies on standard numerical methods common to lattice
implementations of classical field dynamics.222For more details, we refer the
reader to Ref. [48, 49, 50] where the same lattice implementation has been
used. We discretize the classical nonlinear equations of motion in both space
and time, solving them in a finite-volume cubic lattice of size $L$ and
$N^{3}$ points uniformly distributed, with periodic boundary conditions
$Y_{Njk}=Y_{0jk}$, $Y_{iNk}=Y_{i0k}$ and $Y_{ijN}=Y_{ij0}$. The lattice
parameters $L$ and $N$ are chosen in such a way that all the relevant scales
are sufficiently covered during the simulation time. In particular, the
infrared (IR) lattice cutoff $\kappa_{\rm lattice,IR}=2\pi/L$ is chosen to be
much smaller than the minimum value of the characteristic momentum scale
during the tachyonic regime, $\kappa_{\rm lattice,IR}\ll\kappa_{*}(z_{\rm
nl})$, while the UV lattice cutoff $\kappa_{\rm
lattice,UV}=\sqrt{3}/2N\,\kappa_{\rm lattice,IR}$ is taken to exceed the
maximum momentum excited by the tachyonic instability, $\kappa_{\rm
lattice,UV}\simeq{\rm few}\times\kappa_{\rm max}(0)$.

A set of classical initial realizations is then evolved using the discretized
equations of motion, with the expansion of the Universe assumed to be kination
dominated. This evolution is completely deterministic, being the random
character restricted only to the initial conditions inherited from the quantum
treatment in Section 3. The numerical integration is started at a fixed time
$z_{i}$ sufficiently advanced to guarantee that a large fraction of the modes
have become classical but well before the time $z_{\rm nl}$,

$\left(1+\frac{z_{\rm
nl}}{\nu}\right)^{\nu+1/2}\simeq\frac{8\pi\nu^{2}}{(2\nu-1)\sqrt{3\lambda(\nu+1)}}\left(1-\frac{1}{4\nu^{2}}\right)^{1/2}\,,$
(4.1)

at which non linearities become relevant and backreaction starts to affect the
analytical results [13].

Figure 1: The amplitude of the WKB phase $|F_{\kappa}|$ in the correlation
matrix (3.12) as a function of the dimensionless momentum $\kappa$ evaluated
at three different times. For the benchmark point $\nu=10$, $\lambda=10^{-4}$,
the quantum to classical transition takes place before the linear
approximation ceases to be valid at $z_{\rm nl}(\nu,\lambda)=11$. The shaded
areas mark the regions $|F_{\kappa}|\leq 10$ (darker) and $|F_{\kappa}|\leq
100$ (lighter). The dashed vertical lines indicate the values $\kappa_{\rm
min}(z)$ and $\kappa_{\rm max}(z)$ in (3.4), with the same color coding as
$F_{\kappa}(z)$.

The amplitude of all modes that have not become classical by the
initialization time $z_{i}$ (namely those with $|F_{\kappa}|\lesssim 10$, cf.
Fig. 1) is set to zero in our simulations. The remaining Fourier components
are initialized according to a Rayleigh distribution with uniform random phase
$\theta_{k}\in[0,2\pi]$ and dispersion $\sigma_{k}^{2}=|f_{k}|^{2}$, i.e.

$P(|Y_{k}|)\,d|Y_{k}|\,d\theta_{k}=\exp\Big{(}\\!-{\frac{|Y_{k}|^{2}}{\sigma_{k}^{2}}}\Big{)}\,{\frac{d|Y_{k}|^{2}}{\sigma_{k}^{2}}}\,{\frac{d\theta_{k}}{2\pi}}\,,\hskip
28.45274ptY_{k}^{\prime}=\frac{F_{k}(z_{i})}{|f_{k}(z_{i})|^{2}}\ Y_{k}\,.$
(4.2)

Since the short distance variations of the field fluctuations are not relevant
for the thermodynamics description we will be interested in, we will make use
of the ergodic postulate to identify ensemble averages of a given operator
${\cal O}(\vec{y},z)$ with volume averages over the lattice volume $V=L^{3}$,

$\langle{\cal O}(z)\rangle=\frac{1}{V}\int_{\,V}d^{3}y\,\,{\cal
O}(\vec{y},z)\,.$ (4.3)

In some cases, we will perform also suitable coarse-grainings on time with
resolution $z_{0}$ 333Note that $z_{0}$ is not necessarily a constant
throughout the time evolution. and/or averages over $J$ realizations of the
stochastic initial conditions,

$\langle{\cal
O}(z)\rangle=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{z_{0}}\int^{z}_{z-z_{0}}dz\,\,\frac{1}{V}\int_{\,V}d^{3}y\,\,{\cal
O}^{(i)}(\vec{y},z)\,.$ (4.4)

Finally, we will define occupation numbers as [51, 52]

${n}_{\vec{\kappa}}(t)=\left(|\Pi_{\vec{\kappa}}|^{2}\,|Y_{\vec{\kappa}}|^{2}\right)^{\frac{1}{2}}\,.$
(4.5)

For non-tachyonic frequencies this definition coincides with that following
from the correspondence principle. Indeed, in the classical limit, the
occupation number of a mode is given by
${n}_{\vec{\kappa}}(t)=A/(2\pi)=\frac{1}{2}|\Pi_{\vec{\kappa}}|_{\rm
max}||Y_{\vec{\kappa}}|_{\rm max}$ with $A$ the area of the ellipse encircled
by the trajectory in phase space. Taking into account that the modes are
oscillating at late times and assuming them to be weakly coupled, we can write
$|\Pi_{\vec{\kappa}}|_{\rm max}=(2\,|{\Pi}_{\vec{\kappa}}|^{2})^{1/2}$ and
$|Y_{\vec{\kappa}}|_{\rm max}=(2|Y_{\vec{\kappa}}|^{2})^{1/2}$, recovering Eq.
(4.5).

###  4.2 Results

Figure 2: Spectator field growth for the benchmark point $\nu=10$,
$\lambda=10^{-4}$. We display 2-dimensional sections of a fraction $V/8$ of
the 3-dimensional simulation volume $V$ at times
$z=5,\,10,\,15,\,20,\,25,\,60$, with the $z$-axis indicating the corresponding
field value at each point in this surface. The black segment in the second
panel indicates the typical scale of the highest peaks.

Figure 3: Domain wall configurations for the benchmark point $\nu=10$,
$\lambda=10^{-4}$, obtained by plotting the 3-dimensional coordinates of the
grid points with vanishing $Y$ at times $z=5,\,10,\,15,\,20,\,25,\,60$. As
before, we display only $1/8$ of the total simulation volume.

To facilitate the comparison with the analytical results of Ref. [13], and
without lack of generality, we concentrate in what follows on a benchmark
point $\nu=10$, $\lambda=10^{-4}$, for which $z_{\rm nl}\simeq 11$.

Our simulations are carried out on a lattice with $N=1024$ and $L=0.185N$. The
typical output is summarized in Fig. 2, where we present several 2-dimensional
slices of the 3-dimensional lattice volume at different times.444Computer
generated movies covering the whole simulation volume can be found at this
URL.  Although the initial conditions leading to these snapshots are random,
the qualitative features displayed there are universal up to small spatial and
temporal shifts. As anticipated above, an almost homogeneous Universe (1st
panel) becomes rapidly separated into regions with positive and negative field
values and typical size $L\lesssim\kappa_{*}^{-1}$, in agreement with the
analytical result in Eq. (3.18) (2nd panel and 3rd panel). These macroscopic
regions are melted at every field oscillation (4th panel) to reappear again
slightly distorted (5th panel). As time proceeds, additional fluctuations are
generated at smaller and smaller scales, giving eventually rise to a rather
diffuse structure (6th panel). The domain walls separating regions of positive
and negative $Y$ values can be explicitly detected. To this end, we plot in
Fig. 3 the 3D coordinates of the grid points where $|Y|$ is close to zero. The
displayed panels mirror the behaviour of Fig. 2. In particular, we observe a
resilient formation and dissolution of domains.

A better characterization of the above symmetry-breaking picture can be
achieved by considering the following observables:

  1. 1.

Probability distribution: The total simulation volume containing the field
value $Y$ at a given time $z$ is displayed in Fig. 4. This probability
distribution is initially peaked at $Y=0$, with a dispersion relation given by
$\sigma_{k}^{2}=|f_{k}|^{2}$. As time passes by, the distribution evolves
towards the effective minima of the potential at large field values,
eventually splitting in two. Contrary to other symmetry breaking scenarios [7,
8, 9, 10], the time-dependent character of the minima prevents the
stabilization around them, giving rise to an oscillatory bi-modal pattern that
survives for quite a number of oscillations.

Figure 4: Evolution of the probability distribution after symmetry breaking
for the benchmark point $\nu=10$, $\lambda=10^{-4}$. We display snapshots at
times $z=5,\,10,\,15,\,20,\,25,\,60$. The vertical green-dashed lines signal
the location of the time-dependent minima, as given in Eq (2.14).

  2. 2.

Time-dependent dispersion: The root mean-square perturbation $Y_{\rm rms}$ is
displayed in Fig. 5 as a function of the dimensionless conformal time $z$ and
the number of $e$-folds of kination, $N_{e}\equiv\ln(a/a_{\rm kin})$. Note
that, as argued in Refs. [7, 8, 13], the picture of a homogeneously
oscillating scalar field is incorrect when applied to spontaneous symmetry
breaking. Indeed, when the initial value of the homogeneous component is close
to zero, the classical rolling of the field is dominated by its quantum
fluctuations, in agreement with the analytical result (3.18).

Figure 5: Evolution of the volume-averaged lattice $Y_{\rm rms}$ for the
benchmark point $\nu=10$, $\lambda=10^{-4}$ as compared to the analytical
result (3.18), with $N_{e}$ the number of $e$-folds of kination.

  3. 3.

Spectral distributions: The evolution of the occupation numbers (4.5) is shown
in Fig. 6. We observe a rapid growth of infrared fluctuations at the early
stages of symmetry breaking followed by a slower motion towards the UV due to
the potential interactions. Note that the displayed behaviour is far from
thermal. Indeed, thermal equilibrium can be only achieved once the occupation
numbers in the infrared have significantly decreased to
$n_{\kappa}\sim\lambda^{-1/4}$; a limit that notably exceeds our simulation
time. On top of that, it is well known that classical lattice simulations
cannot completely describe the approach to equilibrium due to the Rayleigh-
Jeans divergence555In the continuum temperature goes to zero; on the lattice
it is cut-off dependent. [53].

A simple characterization of the UV cascade can be obtained by computing the
average wavenumber [54]

$\kappa_{c}\equiv\frac{\int_{\kappa>\kappa_{0}}d^{3}\kappa\,|\kappa|\,|\Pi_{\vec{\kappa}}|^{2}}{\int_{\kappa>\kappa_{0}}d^{3}\kappa\,|\Pi_{\vec{\kappa}}|^{2}}\,,$
(4.6)

where we have intentionally excluded the infrared modes filled at the onset of
the simulation, namely those below $\kappa_{0}=0.8$. As shown in Fig. 7, the
momentum scale $\kappa_{c}$ is always smaller than the UV lattice cutoff,
confirming that all relevant scales are well resolved throughout the entire
simulation time.

Figure 6: Time evolution of the occupation numbers (4.5) for the benchmark
point $\nu=10$, $\lambda=10^{-4}$. Different colors correspond to different
times, as indicated in the figure. Figure 7: Time evolution of the UV average
momentum (4.6), with $N_{e}$ the number of $e$-folds of kination.

  4. 4.

Energy budget: The evolution of the kinetic, gradient, potential and
interaction contributions to the total spectator energy density is displayed
in Fig. 8. After a rapid growth at early times and a violent oscillatory
regime completely dominated by the interaction term, the system approaches an
asymptotic state characterized by a large kinetic energy density. The smooth
UV cascade in Fig. 6 translates then into a steady flow of energy from the
interaction and potential terms to a slowly growing gradient component. As
time passes by, the system relaxes towards a virialized state with
$\braket{\textrm{K}}\simeq\braket{\textrm{G}}+2\braket{\textrm{V}}$, in
agreement with the expectation for a quartic potential [55, 54].

Figure 8: Evolution of the volume-averaged kinetic, gradient, potential and
interaction contributions to the total spectator energy density for the
benchmark point $\nu=10$, $\lambda=10^{-4}$, with $N_{e}$ the number of
$e$-folds of kination. In order to facilitate the comparison, we display the
absolute value of the interaction term $I_{1}$, since this is not positive
definite at all times.

  5. 5.

Equation-of-state parameter: The behaviour of the different energy density
contributions is reflected also on the evolution of the associated equation-
of-state parameter

$w_{\chi}\equiv\frac{\braket{p_{\chi}}}{\braket{\rho_{\chi}}}\,,$ (4.7)

which we display in Fig. 9. Since this quantity is rapidly oscillating as
compared to the characteristic expansion rate, we perform here an additional
coarse-graining in time, along the lines of Eq. (4.4). As shown in the left
panel of the figure, and in agreement with the analytical estimate in Appendix
A, the averaged equation of state is initially negative, becoming positive
only during the tachyonic-to-oscillatory transition. The small discrepancy
between the lattice and the analytical estimate can be traced back to the fact
that the latter assumes exact homogeneity and complete irrelevance of the non-
linear potential terms while, albeit small initially, gradients are always
present in the simulation and the potential terms can contribute. As shown in
the right panel, $w_{\chi}$ approaches a value rather close to that of a
radiation fluid ($w_{\chi}=1/3$) at late times. This behaviour is
intrinsically related to the virialization process displayed in Fig. 8 and
completely independent of whether the system is already thermalized or not.

Figure 9: (Left) Evolution of the volume- and time-averaged equation-of-state
parameter (4.7) immediately after symmetry breaking, with $N_{e}$ the number
of $e$-folds of kination. (Right) Evolution of the same quantity at the latest
stages of our simulation. The red horizontal line corresponds to the radiation
value $1/3$.

  6. 6.

Non-Gaussianities: To characterize the level of non-gaussianity in our
scenario, we consider the so-called excess kurtosis, customarily defined as

$\Delta K\equiv\frac{\langle Y^{4}\rangle}{\langle Y^{2}\rangle^{2}}-3\,.$
(4.8)

Figure 10: Time-averaged evolution of the excess kurtosis (4.8) for the
benchmark point $\nu=10$, $\lambda=10^{-4}$, with $N_{e}$ the number of
$e$-folds of kination.

The above results should be understood as robust against variations of the
initializing time $z_{i}$ and the lattice cutoffs. In particular, the temporal
evolution of the different energy components is almost insensitive to the
chosen initial condition and dominated by the modes lying far from the edges,
where cutoff effects are subdominant.

##  5 Implications for the onset of radiation domination

The existence of a kination dominated era in quintessential inflation
scenarios has strong implications for the onset of radiation domination [56,
33]. In particular, even if the energy density transferred to the spectator
field,

$\rho_{\chi}(z)=\mathcal{X}\,\rho_{Y}(z)=16\nu^{2}H_{\rm
kin}^{4}\,\left(\frac{a_{\rm kin}}{a(z)}\right)^{4}\rho_{Y}(z)\,,$ (5.1)

is initially very small, it will eventually become the dominant energy
component due to the rapid decrease of the background energy density during
this period,

$\rho_{B}(z)=3M_{P}^{2}H^{2}(z)\equiv\mathcal{X}\rho^{(r)}_{B}(z)\,\,,\hskip
28.45274pt\rho^{(r)}_{B}(z)=\frac{3}{16\nu^{4}}\left(\frac{M_{P}}{H_{\rm
kin}}\right)^{2}\left(\frac{a_{\rm kin}}{a(z)}\right)^{2}\,.$ (5.2)

The comparison between these two quantities for the benchmark point $\nu=10$,
$\lambda=10^{-4}$ is displayed in Fig. 11 for different values of the Hubble
rate at the onset of kinetic domination. As shown there, even for tiny initial
ratios $\rho_{\chi}(z_{i})/\rho_{B}(z_{i})\sim{\cal O}(10^{-9})$, the Universe
enters radiation domination in less than $10$ $e$-folds of kination.

Figure 11: Comparison between the background energy density $\rho_{B}(z)$ in
Eq. (5.2) for different $H_{\rm kin}$ scales and that transferred to the
spectator field via the Hubble-induced symmetry breaking for the benchmark
point $\nu=10$, $\lambda=10^{-4}$, with $N_{e}$ the number of $e$-folds of
kination.

It is important to emphasize at this point that the onset of radiation
domination is not related to the shape of the potential, as one could naively
expect by considering the homogeneous approximation [16], but rather to the
late-time hierarchy $\braket{V}\ll\braket{K},\braket{G}$ observed in Fig. 8.
Indeed, for any virialized late-time scenario,
[$1/2\langle{\dot{\chi}^{2}}\rangle=1/2\langle{(\nabla\chi/a)^{2}}\rangle+n\langle{V}\rangle$],
involving a monomial potential $V\propto\chi^{2n}$ with $n\geq 1$ and
negligible Hubble-induced interactions, one has [57, 58]

$w_{\chi}=\frac{1}{3}+\frac{2}{3}\frac{(n-2)}{(n+1)+\langle(\nabla\chi/a)^{2}\rangle/\langle
V\rangle}\,,$ (5.3)

and therefore $w_{\chi}\to 1/3$ if $\langle(\nabla\chi/a)^{2}\rangle\gg\langle
V\rangle$, in clear contrast with the value $w_{\chi}\to(n-1)/(n+1)$ obtained
in the homogeneous limit $\langle(\nabla\chi/a)^{2}\rangle\ll\langle V\rangle$
[59, 60]. This can be explicitly verified by performing, for instance, a
lattice simulation for a sextic potential
$\frac{\sigma}{6\Lambda^{2}}\chi^{6}$, with $\sigma$ a coupling constant and
$\Lambda$ a cutoff scale assumed to exceed $H_{\rm kin}$. In terms of the
rescaled variables (2.9), the corresponding Klein-Gordon equation and
potential energy density take the form

$Y^{\prime\prime}-\nabla^{2}\,Y-M^{2}(z)Y+\sigma_{\rm eff}\,\,{\cal
H}(z)\,Y^{5}=0\,,\hskip 28.45274ptV\equiv\frac{1}{6}\,\sigma_{\rm
eff}\,\,{\cal H}(z)\,Y^{6}\,,$ (5.4)

with

$\sigma_{\rm eff}\equiv\sigma(2\nu)^{3}\left(\frac{H_{\rm
kin}}{\Lambda}\right)^{2}$ (5.5)

an effective dimensionless coupling constant.

Figure 12: Energy distribution (left) and volume- and time-averaged equation
of state (4.7) (right) for the sextic scenario in (5.4). The asymptotic
behaviour of $w_{\chi}$ differs clearly from the homogeneous expectation
$w_{\chi}=1/2$, approaching a radiation-like value $w_{\chi}=1/3$. This
illustrates the importance of accounting for the production of IR
perturbations and their subsequent UV cascade when studying Ricci reheating
scenarios.

The evolution of the corresponding equation-of-state parameter $w_{\chi}$ for
fiducial values $\nu=10$, $\sigma=0.01$ and $\Lambda=100\,H_{\rm kin}$
($\sigma_{\rm eff}=0.008$) is displayed in Fig. 12. As expected, this quantity
approaches a radiation-like value $w_{\chi}=1/3$ rather than the
$w_{\chi}=1/2$ one following from the homogeneous approximation. Combined with
kination, this observation allows for the generic onset of radiation
domination for general interacting potentials, significantly extending the
results of Ref. [16].

##  6 Assumptions and future research lines

We have performed several approximations to reduce the physical scenario to a
baseline model that could be treated in a simple way with analytical and
numerical techniques. Relaxing any of these assumptions opens new avenues for
model building.

  1. 1.

Symmetry group: We have assumed the spectator field to display a $Z_{2}$
symmetry, understanding this choice as a mere proof of concept of a much more
general setting. In particular, we expect our results regarding the impact of
perturbations on symmetry breaking and virialization to be extensible to more
involved symmetry groups, provided that these are global. The extension to
local symmetries could, on the other hand, lead to new phenomenology,
certainly worthy to explore.

  2. 2.

Expansion of the Universe: We have assumed the inflation-to-kination
transition to be instantaneous as compared to the temporal evolution of the
spectator field. In practice, this corresponds to a quench approximation

$R(t<t_{\rm kin})>0\,,\hskip 28.45274ptR(t_{\rm kin}<t<t_{r})<0\,,$ (6.1)

with $t_{r}$ the symmetry restoration time. On top of that, we assumed the
spectator energy density to be subdominant during kination, such that the
Ricci scalar can be well approximated by an homogeneous function of time
depending only on the inflaton equation of state. While this is certainly a
reasonable assumption for the main stages described in this paper, it should
be eventually taken into account for a precise determination of the onset of
radiation domination.

  3. 3.

Large $\nu$ approximation: In order to minimize the impact of horizon effects,
we have focused on a large $\nu$ limit leading to a reasonable separation
between ${\cal H}(z)$ and the typical momentum of the spectator field
fluctuations, $\kappa_{*}(z)\equiv 2\sqrt{\nu+1}\,{\cal H}(z)$. This makes the
integrated contribution of superhorizon modes reentering the horizon less and
less important as compared with that of subhorizon scales never amplified by
the inflationary dynamics, justifying the choice of a vacuum initial state
with zero occupation numbers. We emphasize, however, that the asymptotic
virialization observed in our simulations possesses an universal character, as
we have explicitly verified by considering alternative initial conditions.

  4. 4.

Parameter space: We have restricted ourselves to model parameters leading to a
maximum $10\%$ correction to the reduced Planck mass, i.e. $\xi\chi_{\rm
max}^{2}/M_{P}^{2}\leq 0.1$ [13]. Combined with the above two assumptions,
this allows us to plausibly ignore metric fluctuations. Even within this self-
consistency regime, our investigations are still limited to a small portion of
the parameter space. Therefore, we cannot a priori exclude the existence of
specific combinations of $\nu$ and $\lambda$ leading to different dynamics. We
emphasize, however, that, given the decaying nature of the non-minimal
coupling to gravity, the asymptotic behavior described in this paper is
expected to be robust against variations of the model parameters. What could
be affected at most is the duration of the tachyonic phase and the amount of
defects that can be generated in the intermediate stages.

  5. 5.

Couplings: We have assumed the spectator field to be decoupled or weakly
coupled to an unspecified inflationary sector. Unless forbidden by any
symmetry, this requires some level of fine-tuning. On top of that, we have
ignored potential couplings to other matter components that must necessarily
exist to recover the Standard Model content and other relevant species as dark
matter [61], implicitly neglecting potential depletion and combined preheating
effects [62, 63, 64, 65] that could modify the approach to equilibrium.

An interesting possibility along this line of thought is the potential
identification of the spectator field with the Standard Model Higgs, as first
suggested in Ref. [14] within the context of kinetic domination. This
association would lead to a rather minimalistic heating scenario not involving
additional degrees of freedom beyond the Standard Model content. One of the
recurrent arguments for restricting the parameter space of this appealing
possibility [14, 16] is the apparent incompatibility of particle physics data
with the absolute stability of the Standard Model vacuum below the Planck
scale [66, 67]. We note, however, that the analyses of Ref. [14, 16] disregard
implicitly the ambiguous relation between the Montecarlo reconstructed top
quark mass and the actual top pole mass/Yukawa coupling entering in the
renormalization group equations. As argued in Ref. [68], once this source of
uncertainty is taken into account, it is perfectly possible for the Standard
Model vacuum to be absolutely stable all the way up till the Planck scale.
Indeed, the recent LHC measurements of the top pole mass point intriguingly
towards this scenario [69, 70].

Given the absolute stability of the Standard Model vacuum, the depletion of
the spectator Higgs field in our setting is expected to proceed along the
lines of Refs. [71, 72, 73, 74, 75], with potential differences associated to
the bimodal pattern in Fig. 4, not accounted for in these papers. The
thermalization of the resulting Standard Model species is estimated to be
almost instantaneous [76, 14, 77, 78]

##  7 Conclusions

In this work, we have made a step forward towards the complete
characterization of symmetry breaking in the context of Hubble-induced phase
transitions, setting the stage for getting more reliable predictions to be
eventually confronted with data. As a case of study, we investigated the full
symmetry-breaking dynamics of a $Z_{2}$ spectator field non-minimally coupled
to gravity. The transition from the unbroken to the broken phase is induced by
a change of sign in the Ricci scalar as the Universe evolves from inflation to
an era in which the equation-of-state parameter is stiff. In order to follow
the non-linear evolution of the spectator field until the development of an
equilibrium configuration, we made use of classical lattice simulations, being
the corresponding initial conditions analytically obtained from a complete
quantum field theory description of the early dynamics.

Our simulations allowed us to identify three different stages on the Hubble-
induced symmetry breaking process. First, when the background acquires a stiff
equation-of-state parameter immediately after the end of inflation, it induces
a tachyonic instability for the spectator field fluctuations, which become
amplified on infrared scales. This translates into the formation of rather
homogeneous regions in real space separated by domain wall configurations.
Second, as the field distribution splits and goes beyond the fading away
minima of the effective potential, there appears a bi-modal oscillatory
pattern that stays in place for a large number of oscillations. This stage is
initially dominated by the non-minimal coupling to gravity, which become
weaker and weaker as the Universe expands. Finally, we observe a long
virialization stage, where the total energy density becomes democratically
distributed among the kinetic and gradient components and the spectrum
develops a turbulent cascade towards the UV. The equation-of-state of the
spectator field approaches then a radiation-like value $w_{\chi}=1/3$,
independently of the shape of the potential. The memory loss induced by this
asymptotic stage does not imply that all the traces of the temporary Hubble-
induced symmetry breaking are completely erased. Among other consequences, the
temporal existence of topological defects may lead to the production of a
detectable gravitational waves’ backgrounds [12] or baryogenesis [11],
depending on the details of the symmetry breaking potential. We postpone the
lattice characterization of these crucial aspects to a further publication.

## Acknowledgments

We thank Mark Hindmarsh for useful discussions on the physics and simulation
of topological defects. DB acknowledges support from the “Atracción del
Talento Científico” en Salamanca programme, from Junta de Castilla y León and
Fondo Europeo de desarrollo regional through the project SA0096P20 and from
project PGC2018-096038-B-I00 by Spanish Ministerio de Ciencia, Innovación y
Universidades. JR acknowledges the support of the Fundação para a Ciência e a
Tecnologia (Portugal) through the CEECIND/01091/2018 grant and Project No.
UIDB/00099/2020. ALE (ORCID ID 0000-0002-1696-3579) is supported by the
National Science Foundation grant PHY-1820872. This work has been possible
thanks to the computational resources on the i2Basque academic network. The
authors acknowledge the Tufts University High Performance Compute Cluster
(https://it.tufts.edu/high-performance-computing) which was also utilized for
the research reported in this paper.

##  Appendix A The homogeneous approximation

For the sake of completeness and in order facilitate the comparison with
previous studies in the literature [14, 15, 16], we discuss here the evolution
of the zero mode in the absence of quantum fluctuations. In this classical and
homogeneous limit, the solution of the equation of motion (2.15) can be
analytically obtained for both the tachyonic and pure quartic regimes.
Assuming initial conditions $Y_{\rm kin}$ and $dY/dz|_{0}=Y_{\rm kin}/\nu$ and
matching the solutions at the time

$\frac{z_{\rm b}}{\nu}=\left(\frac{8}{\lambda Y_{\rm
kin}^{2}}\frac{2\nu-1}{2\nu+1}\right)^{1/(2\nu+3)}-1\,,$ (A.1)

at which the quadratic and quartic terms in the $Y$-field potential become
equal, we get [13]

$Y(z)\simeq\begin{cases}\frac{2\nu+1}{4\nu}Y_{\rm
kin}\left(1+\frac{z}{\nu}\right)^{\nu+1/2}\hskip 68.2866pt{\rm for}\hskip
14.22636ptz\leq z_{b}\,,\\\ Y_{\rm
osc}\,\mathrm{cn}\left[\sqrt{\lambda}\,Y_{\rm osc}\,(z-z_{\rm
osc}),\frac{1}{\sqrt{2}}\right]\hskip 34.1433pt{\rm for}\hskip
14.22636ptz>z_{b}\,,\end{cases}$ (A.2)

where in the first expression we have neglected a subdominant decaying mode.
Here $\mathrm{cn}$ stands for the Jacobi elliptic cosine,

$\displaystyle Y_{\rm
osc}=\frac{Y(z_{b})}{\sqrt{2}}\left[1+\left(1+\frac{2}{\beta^{2}}\frac{2\nu+1}{2\nu-1}\right)^{1/2}\right]^{1/2}$
(A.3)

denotes the amplitude of the field at the time

$\displaystyle z_{\rm osc}=z_{b}+\frac{1}{\,\beta\sqrt{\lambda}Y_{\rm
osc}}\arccos\left(\frac{Y(z_{b})}{Y_{\rm osc}}\right)$ (A.4)

at which the oscillatory part reaches its maximum ($z_{\rm osc}>z_{b}$) and
$\beta=2\pi/T_{x}=0.8472$ is a numerical factor associated with the leading
frequency of oscillation in the elliptic cosine series [13].

Although the solution (A.2) will break down as soon as gradients start to
develop, it can be understood a relatively good description of the complete
scenario at very early times, when all field fluctuations are still very close
to the maximum of the effective potential. In particular, we can use it to
estimate the value of the spectator field equation-of-state parameter (4.7) at
the onset of kinetic domination. Taking into account Eqs. (2.16) and (2.17) at
zero gradients, we get666If the decaying mode for the tachyonic regime in
(A.2) is not omitted, we rather get  $w_{\chi,{\rm
ini}}=\frac{3+12\nu^{2}-32\nu^{4}}{3+36\nu^{2}}\,,$ which for the benchmark
value $\nu=10$ translates into $w_{\chi,{\rm ini}}=-88.5$. Note, however, that
the growing mode becomes completely dominant already before the time $z_{i}$
at which we start our lattice simulations.

$w_{\chi,{\rm ini}}=1-\frac{4\nu}{3}\,,$ (A.5)

which translates into a value $w_{\chi,{\rm ini}}=-12.3$ for $\nu=10$, in
reasonable agreement with our lattice results and Ref. [16].

## References

  * [1] D. Cormier and R. Holman, _Spinodal inflation_ , _Phys. Rev. D_ 60 (1999) 041301, [hep-ph/9812476].
  * [2] S. Tsujikawa and T. Torii, _Spinodal effect in the natural inflation model_ , _Phys. Rev. D_ 62 (2000) 043505, [hep-ph/9912499].
  * [3] J. Rubio and E. S. Tomberg, _Preheating in Palatini Higgs inflation_ , _JCAP_ 04 (2019) 021, [1902.10148].
  * [4] A. Karam, E. Tomberg and H. Veermäe, _Tachyonic Preheating in Palatini $R^{2}$ Inflation_, 2102.02712.
  * [5] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso and D. Podolsky, _Preheating with trilinear interactions: Tachyonic resonance_ , _JCAP_ 0607 (2006) 006, [hep-ph/0602144].
  * [6] A. A. Abolhasani, H. Firouzjahi and M. M. Sheikh-Jabbari, _Tachyonic Resonance Preheating in Expanding Universe_ , _Phys. Rev. D_ 81 (2010) 043524, [0912.1021].
  * [7] G. N. Felder, J. Garcia-Bellido, P. B. Greene, L. Kofman, A. D. Linde and I. Tkachev, _Dynamics of symmetry breaking and tachyonic preheating_ , _Phys. Rev. Lett._ 87 (2001) 011601, [hep-ph/0012142].
  * [8] G. N. Felder, L. Kofman and A. D. Linde, _Tachyonic instability and dynamics of spontaneous symmetry breaking_ , _Phys. Rev._ D64 (2001) 123517, [hep-th/0106179].
  * [9] E. J. Copeland, S. Pascoli and A. Rajantie, _Dynamics of tachyonic preheating after hybrid inflation_ , _Phys. Rev._ D65 (2002) 103517, [hep-ph/0202031].
  * [10] J. Garcia-Bellido, M. Garcia Perez and A. Gonzalez-Arroyo, _Symmetry breaking and false vacuum decay after hybrid inflation_ , _Phys. Rev._ D67 (2003) 103501, [hep-ph/0208228].
  * [11] D. Bettoni and J. Rubio, _Quintessential Affleck-Dine baryogenesis with non-minimal couplings_ , _Phys. Lett._ B784 (2018) 122–129, [1805.02669].
  * [12] D. Bettoni, G. Doménech and J. Rubio, _Gravitational waves from global cosmic strings in quintessential inflation_ , _JCAP_ 1902 (2019) 034, [1810.11117].
  * [13] D. Bettoni and J. Rubio, _Hubble-induced phase transitions: Walls are not forever_ , _JCAP_ 2001 (2020) 002, [1911.03484].
  * [14] D. G. Figueroa and C. T. Byrnes, _The Standard Model Higgs as the origin of the hot Big Bang_ , _Phys. Lett._ B767 (2017) 272–277, [1604.03905].
  * [15] K. Dimopoulos and T. Markkanen, _Non-minimal gravitational reheating during kination_ , _JCAP_ 1806 (2018) 021, [1803.07399].
  * [16] T. Opferkuch, P. Schwaller and B. A. Stefanek, _Ricci Reheating_ , _JCAP_ 1907 (2019) 016, [1905.06823].
  * [17] M. Fairbairn, K. Kainulainen, T. Markkanen and S. Nurmi, _Despicable Dark Relics: generated by gravity with unconstrained masses_ , _JCAP_ 1904 (2019) 005, [1808.08236].
  * [18] L. Laulumaa, T. Markkanen and S. Nurmi, _Primordial dark matter from curvature induced symmetry breaking_ , _JCAP_ 08 (2020) 002, [2005.04061].
  * [19] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , 2006.16182.
  * [20] R. Micha and I. I. Tkachev, _Turbulent thermalization_ , _Phys. Rev._ D70 (2004) 043538, [hep-ph/0403101].
  * [21] R. Micha and I. I. Tkachev, _Relativistic turbulence: A Long way from preheating to equilibrium_ , _Phys. Rev. Lett._ 90 (2003) 121301, [hep-ph/0210202].
  * [22] C. Wetterich, _Cosmology and the Fate of Dilatation Symmetry_ , _Nucl. Phys._ B302 (1988) 668–696, [1711.03844].
  * [23] C. Wetterich, _The Cosmon model for an asymptotically vanishing time dependent cosmological ’constant’_ , _Astron. Astrophys._ 301 (1995) 321–328, [hep-th/9408025].
  * [24] P. J. E. Peebles and A. Vilenkin, _Quintessential inflation_ , _Phys. Rev._ D59 (1999) 063505, [astro-ph/9810509].
  * [25] B. Spokoiny, _Deflationary universe scenario_ , _Phys. Lett._ B315 (1993) 40–45, [gr-qc/9306008].
  * [26] P. Brax and J. Martin, _Coupling quintessence to inflation in supergravity_ , _Phys. Rev._ D71 (2005) 063530, [astro-ph/0502069].
  * [27] C. Wetterich, _Variable gravity Universe_ , _Phys. Rev._ D89 (2014) 024005, [1308.1019].
  * [28] C. Wetterich, _Inflation, quintessence, and the origin of mass_ , _Nucl. Phys._ B897 (2015) 111–178, [1408.0156].
  * [29] M. W. Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, _Variable gravity: A suitable framework for quintessential inflation_ , _Phys. Rev._ D90 (2014) 023512, [1402.6661].
  * [30] A. Agarwal, R. Myrzakulov, M. Sami and N. K. Singh, _Quintessential inflation in a thawing realization_ , _Phys. Lett._ B770 (2017) 200–208, [1708.00156].
  * [31] C.-Q. Geng, C.-C. Lee, M. Sami, E. N. Saridakis and A. A. Starobinsky, _Observational constraints on successful model of quintessential Inflation_ , _JCAP_ 1706 (2017) 011, [1705.01329].
  * [32] K. Dimopoulos and C. Owen, _Quintessential Inflation with $\alpha$-attractors_, _JCAP_ 1706 (2017) 027, [1703.00305].
  * [33] J. Rubio and C. Wetterich, _Emergent scale symmetry: Connecting inflation and dark energy_ , _Phys. Rev._ D96 (2017) 063509, [1705.00552].
  * [34] K. Dimopoulos, L. Donaldson Wood and C. Owen, _Instant preheating in quintessential inflation with $\alpha$-attractors_, _Phys. Rev._ D97 (2018) 063525, [1712.01760].
  * [35] Y. Akrami, R. Kallosh, A. Linde and V. Vardanyan, _Dark energy, $\alpha$-attractors, and large-scale structure surveys_, _JCAP_ 1806 (2018) 041, [1712.09693].
  * [36] C. García-García, E. V. Linder, P. Ruíz-Lapuente and M. Zumalacárregui, _Dark energy from $\alpha$-attractors: phenomenology and observational constraints_, _JCAP_ 1808 (2018) 022, [1803.00661].
  * [37] L. H. Ford, _Cosmological constant damping by unstable scalar fields_ , _Phys. Rev. D_ 35 (1987) 2339.
  * [38] L. H. Ford and T. A. Roman, _Classical scalar fields and violations of the second law_ , _Phys. Rev. D_ 64 (2001) 024023, [gr-qc/0009076].
  * [39] J. D. Bekenstein, _Nonsingular General Relativistic Cosmologies_ , _Phys. Rev. D_ 11 (1975) 2072–2075.
  * [40] E. E. Flanagan and R. M. Wald, _Does back reaction enforce the averaged null energy condition in semiclassical gravity?_ , _Phys. Rev. D_ 54 (1996) 6233–6283, [gr-qc/9602052].
  * [41] D. Polarski and A. A. Starobinsky, _Semiclassicality and decoherence of cosmological perturbations_ , _Class. Quant. Grav._ 13 (1996) 377–392, [gr-qc/9504030].
  * [42] J. Lesgourgues, D. Polarski and A. A. Starobinsky, _Quantum to classical transition of cosmological perturbations for nonvacuum initial states_ , _Nucl. Phys._ B497 (1997) 479–510, [gr-qc/9611019].
  * [43] C. Kiefer and D. Polarski, _Emergence of classicality for primordial fluctuations: Concepts and analogies_ , _Annalen Phys._ 7 (1998) 137–158, [gr-qc/9805014].
  * [44] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_.  Dover, New York, ninth dover printing, tenth gpo printing ed., 1964.
  * [45] J. M. Cornwall, R. Jackiw and E. Tomboulis, _Effective Action for Composite Operators_ , _Phys. Rev. D_ 10 (1974) 2428–2445.
  * [46] J. Berges and J. Cox, _Thermalization of quantum fields from time reversal invariant evolution equations_ , _Phys. Lett. B_ 517 (2001) 369–374, [hep-ph/0006160].
  * [47] A. Arrizabalaga, J. Smit and A. Tranberg, _Tachyonic preheating using 2PI-1/N dynamics and the classical approximation_ , _JHEP_ 10 (2004) 017, [hep-ph/0409177].
  * [48] A. Lopez-Eiguren, J. Lizarraga, M. Hindmarsh and J. Urrestilla, _Cosmic Microwave Background constraints for global strings and global monopoles_ , _JCAP_ 07 (2017) 026, [1705.04154].
  * [49] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren and J. Urrestilla, _Scaling Density of Axion Strings_ , _Phys. Rev. Lett._ 124 (2020) 021301, [1908.03522].
  * [50] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren and J. Urrestilla, _Approach to scaling in axion string networks_ , _Phys. Rev. D_ 103 (2021) 103534, [2102.07723].
  * [51] M. Salle, J. Smit and J. C. Vink, _Thermalization in a Hartree ensemble approximation to quantum field dynamics_ , _Phys. Rev. D_ 64 (2001) 025016, [hep-ph/0012346].
  * [52] M. Salle, J. Smit and J. C. Vink, _Staying thermal with Hartree ensemble approximations_ , _Nucl. Phys. B_ 625 (2002) 495–511, [hep-ph/0012362].
  * [53] J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Basin of attraction for turbulent thermalization and the range of validity of classical-statistical simulations_ , _JHEP_ 05 (2014) 054, [1312.5216].
  * [54] C. Destri and H. J. de Vega, _Ultraviolet cascade in the thermalization of the classical phi**4 theory in 3+1 dimensions_ , _Phys. Rev. D_ 73 (2006) 025014, [hep-ph/0410280].
  * [55] D. Boyanovsky, C. Destri and H. J. de Vega, _The Approach to thermalization in the classical phi**4 theory in (1+1)-dimensions: Energy cascades and universal scaling_ , _Phys. Rev. D_ 69 (2004) 045003, [hep-ph/0306124].
  * [56] G. N. Felder, L. Kofman and A. D. Linde, _Inflation and preheating in NO models_ , _Phys. Rev._ D60 (1999) 103505, [hep-ph/9903350].
  * [57] K. D. Lozanov and M. A. Amin, _Equation of State and Duration to Radiation Domination after Inflation_ , _Phys. Rev. Lett._ 119 (2017) 061301, [1608.01213].
  * [58] K. D. Lozanov and M. A. Amin, _Self-resonance after inflation: oscillons, transients and radiation domination_ , _Phys. Rev. D_ 97 (2018) 023533, [1710.06851].
  * [59] M. S. Turner, _Coherent Scalar Field Oscillations in an Expanding Universe_ , _Phys. Rev._ D28 (1983) 1243.
  * [60] M. C. Johnson and M. Kamionkowski, _Dynamical and Gravitational Instability of Oscillating-Field Dark Energy and Dark Matter_ , _Phys. Rev. D_ 78 (2008) 063010, [0805.1748].
  * [61] N. Bernal, J. Rubio and H. Veermäe, _Boosting Ultraviolet Freeze-in in NO Models_ , _JCAP_ 06 (2020) 047, [2004.13706].
  * [62] J. Garcia-Bellido, D. G. Figueroa and J. Rubio, _Preheating in the Standard Model with the Higgs-Inflaton coupled to gravity_ , _Phys. Rev._ D79 (2009) 063531, [0812.4624].
  * [63] J. Rubio, _Higgs inflation and vacuum stability_ , _J. Phys. Conf. Ser._ 631 (2015) 012032, [1502.07952].
  * [64] J. Repond and J. Rubio, _Combined Preheating on the lattice with applications to Higgs inflation_ , _JCAP_ 1607 (2016) 043, [1604.08238].
  * [65] J. Fan, K. D. Lozanov and Q. Lu, _Spillway Preheating_ , _JHEP_ 05 (2021) 069, [2101.11008].
  * [66] F. Bezrukov, M. Y. Kalmykov, B. A. Kniehl and M. Shaposhnikov, _Higgs Boson Mass and New Physics_ , _JHEP_ 10 (2012) 140, [1205.2893].
  * [67] G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori et al., _Higgs mass and vacuum stability in the Standard Model at NNLO_ , _JHEP_ 08 (2012) 098, [1205.6497].
  * [68] F. Bezrukov and M. Shaposhnikov, _Why should we care about the top quark Yukawa coupling?_ , _J. Exp. Theor. Phys._ 120 (2015) 335–343, [1411.1923].
  * [69] ATLAS collaboration, G. Aad et al., _Measurement of the top-quark mass in $t\bar{t}+1$-jet events collected with the ATLAS detector in $pp$ collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 11 (2019) 150, [1905.02302].
  * [70] CMS collaboration, A. M. Sirunyan et al., _Measurement of $\mathrm{t\bar{t}}$ normalised multi-differential cross sections in pp collisions at $\sqrt{s}=13$ TeV, and simultaneous determination of the strong coupling strength, top quark pole mass, and parton distribution functions_, _Eur. Phys. J. C_ 80 (2020) 658, [1904.05237].
  * [71] K. Enqvist, T. Meriniemi and S. Nurmi, _Generation of the Higgs Condensate and Its Decay after Inflation_ , _JCAP_ 10 (2013) 057, [1306.4511].
  * [72] K. Enqvist, S. Nurmi and S. Rusak, _Non-Abelian dynamics in the resonant decay of the Higgs after inflation_ , _JCAP_ 10 (2014) 064, [1404.3631].
  * [73] D. G. Figueroa, J. Garcia-Bellido and F. Torrenti, _Decay of the standard model Higgs field after inflation_ , _Phys. Rev._ D92 (2015) 083511, [1504.04600].
  * [74] K. Enqvist, S. Nurmi, S. Rusak and D. Weir, _Lattice Calculation of the Decay of Primordial Higgs Condensate_ , _JCAP_ 1602 (2016) 057, [1506.06895].
  * [75] D. G. Figueroa, A. Rajantie and F. Torrenti, _Higgs field-curvature coupling and postinflationary vacuum instability_ , _Phys. Rev. D_ 98 (2018) 023532, [1709.00398].
  * [76] F. Bezrukov, J. Rubio and M. Shaposhnikov, _Living beyond the edge: Higgs inflation and vacuum metastability_ , _Phys. Rev. D_ 92 (2015) 083512, [1412.3811].
  * [77] P. B. Arnold, G. D. Moore and L. G. Yaffe, _Effective kinetic theory for high temperature gauge theories_ , _JHEP_ 01 (2003) 030, [hep-ph/0209353].
  * [78] A. Kurkela and G. D. Moore, _Thermalization in Weakly Coupled Nonabelian Plasmas_ , _JHEP_ 12 (2011) 044, [1107.5050].


# Suppression of stacking order with doping in 1T-TaS2-xSex

Sharon S. Philip  University of Virginia, Charlottesville, VA 22904, USA
Despina Louca  University of Virginia, Charlottesville, VA 22904, USA  J. C.
Neuefeind  Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA  Matthew B.
Stone  Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA  A. I.
Kolesnikov  Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA

###### Abstract

In 1T-TaS2-2xSe2x, the charge density wave (CDW) state features a star of
David lattice that expands across layers as the system becomes commensurate on
cooling. The layers can also order along the c-axis and different stacking
orders have been proposed. Using neutron scattering on powder samples, we
compared the stacking order previously observed in 1T-TaS2 as the system is
doped with Se. While at low temperature, a 13c layer sequence stacking was
observed in TaS2, this type of ordering was not evident with doping. Doping
with Se results in a nearly commensurate state with the Mott state suppressed
which may be linked to the absence of the layer stacking.

††preprint: APS/123-QED

##  I Introduction

Quasi-two dimensional (2D) in nature, transition metal dichalcogenides (TMDs)
1T-MX2 (M = Ti, Ta and X = S, Se, Te) are prone to electronic instabilities
Wilson  _et al._ (1975). 1T-Ta(S/Se)2 exhibits an incredibly rich phase
diagram with multiple charge density wave (CDW) transitions emerging as a
function of temperature and upon doping. In 1T-TaS2-2xSe2x, macroscopic
behaviors such as CDW and superconductivity Sun (2013); Liu  _et al._ (2016);
Wang  _et al._ (2020); Sipos  _et al._ (2008) have been observed, and more
recently, a quantum spin liquid (QSL) has been proposed in 1T-TaS2 as well Law
and Lee (2017). In the typical Peierls model for CDW order Rossnagel (2011),
the instability of the coupled electron-lattice system brings a structural
phase transition that is driven by strong electron-phonon coupling Chen  _et
al._ (2023); Sun  _et al._ (2018); C. J. Butler (2020); Zhang and Wu (2023).
The CDW formation can bring electron localization where displacements along
phonon modes lower the total electronic energy by opening up a gap at the
Fermi level, EF Warawa  _et al._ (2023); Hansen  _et al._ (2023). This
scenario, although applicable to simple one-dimensional systems, does not
fully describe the case of 1T-TaS2-xSex where the CDW behavior is intertwined
with the opening of a Mott gap Manzke  _et al._ (1989). The origin of the CDW
has been highly debated in TMDs Wegner  _et al._ (2020); Chatterjee  _et al._
(2015); van Wezel  _et al._ (2010). The Fermi surface nesting scenario most
often does not apply. Existing models of the CDW order are broadly classified
into three types: in one, it involves an excitonic condensation mechanism; in
two, it involves a Jahn-Teller-like distortion mechanism; and in three, it
involves a hybrid model, a combination of Jahn-Teller and exciton
condensationvan Wezel  _et al._ (2010).

1T-TaS2 exhibits a strong CDW instability and electronic localization that
leads to several interesting effects. Upon cooling from high temperatures,
three main phases form: the high temperature incommensurate CDW (ICDW), the
intermediate temperature nearly commensurate CDW (NCCDW) and the low
temperature commensurate CDW (CCDW) Ang  _et al._ (2015). The ICDW appears
below 540 K on cooling from the high temperature metallic state, with a
transition from the $P\overline{3}m1$ crystal symmetry shown in Fig. 1(a) to
the $P\overline{3}$ structure shown in Fig. 1(b). This transition leads to
displacements of Ta ions that gives rise to the well-known star of David
motifs. Upon cooling from the normal, high temperature metallic state,
systematic displacements of the transition metal Ta leads to a star of David
formation consisting of 13 Ta ions, in-plane. Domains of these formations
expand to a commensurate CDW phase on cooling. Distinct from other CDW
systems, in 1T-TaS2, the commensurate CDW state is accompanied by a metal-
insulator (MI) transition that has been proposed to arise either from Mott
localization or from disorder induced Anderson localization. Important to the
MI behavior are the orbital ordering and out of plane correlations, as well as
layer stacking order.

In the ICDW, the stars have limited ordering in-plane. Further cooling leads
to the ICDW becoming NCCDW at T=350 K, where the $\sqrt{13}\cdot\sqrt{13}$
structural modulation first appears with a 12o tilt relative to the original
ab-plane. An expansion of the star of David motifs occurs in-plane Spijkerman
_et al._ (1997). Below 180 K, the $\sqrt{13}\cdot\sqrt{13}$ structural
modulation persists with a rotation of 13.9o relative to the plane while the
CDW becomes commensurate. The steps in the CDW transitions coincide with the
kinks observed in the transport Fazekas and Tosatti (1979) as the system goes
from the metallic to the insulating state. On the other end of the phase
diagram, in 1T-TaSe2, the CCDW sets in at T $\approx$ 430 K, and the system
shows no MI transition. It remains metallic down to the lowest temperature
according to transport data Wilson  _et al._ (1974). Between these two ends,
superconductivity emerges upon doping that coexists with a broad NCCDW region
in the phase diagram Ang  _et al._ (2015). The coexistence of
superconductivity with CDW domains has been observed in other TMDs such as in
the 2H polytype and in other systems such as the cuprates Canfield  _et al._
(1998); Kvashnin  _et al._ (2020); Egami  _et al._ (2000); Lee  _et al._
(2014).

The electronic bands appear to undergo a continuous change with decreasing
temperature in going though the many transition steps Smith  _et al._ (1985);
Aiura  _et al._ (2004); Zwick  _et al._ (1998). In the absence of high
temperature angle resolved photoemission spectroscopy (ARPES) due to
resolution, there is no apparent nesting of the Fermi surface and a CDW gap is
not necessarily located at the $\Gamma$ point. Measurements suggested that the
gap appears elsewhere in k-space Rossnagel (2011); Fei  _et al._ (2022). The
domain-like CDW structures of the NCCDW and ICDW states in 1T-TaS2 are
discommensurate and semi-metallic, but when the CDW becomes commensurate, the
Fermi surface disappears and either a Mott-Hubbard localization or a disorder
induced Anderson localization sets in Fei  _et al._ (2022). Across the NCCDW-
CCDW boundary, the Fermi surface is continuously reduced. This effect is
convoluted by d-electron localization that opens up an energy gap. In the CCDW
phase, the gap is fully present, leading to a semiconducting state with about
a 200 meV bandgap Rossnagel (2011).

In the normal phase above 540 K, the Ta 5d band at the $\Gamma$ point should
be above EF Bovet  _et al._ (2004). As the system goes through the NCCDW
phase, this band becomes visible with ARPES. Further cooling to the CCDW
state, band folding is observed because of the smaller Brillouin zone between
180 and 160 K, and an abrupt energy shift occurs due to opening of the energy
gap Wang  _et al._ (2020). The loss of the Fermi surface continues with
further cooling while the CDW gap continues to grow. The first order
transition seen in the transport at 180 K on cooling is most likely due to a
Mott-Hubbard localization Ritschel  _et al._ (2015, 2018a). On warming, a
different behavior is observed where the resistivity exhibits a hysteresis,
with its value dropping 280 K, marking the CCDW-NCCDW transition. This has
been attributed to be due to changes in the c-axis stacking order Lee  _et
al._ (2019); Wang  _et al._ (2020).

We report on the nature of the layer stacking order with temperature and
doping. Earlier, we observed that the c-axis expands in the CCDW phase of
1T-TaS2 on warming but drops at the crossover between the CCDW-NCCDW
transition Philip  _et al._ (2023a). It has been suggested that the
localization of the d-electrons that brings the gap in the electronic
structure depends on the expansion of the c-axis Ritschel  _et al._ (2018b,
2015). This in turn is related to the c-axis stacking order where changes in
the interlayer coupling might drive the Mott transition. Neutron diffraction
measurements confirmed the presence of 13c stacking order that disappears on
warming across the CCDW-NCDW transition in 1T-TaS2. The appearance of the 13c
layer sequence is expected to drive the Mott localization. The 13c stacking
sequence was previously suggested in Ref. Scruby  _et al._ (1975) from X-ray
diffraction data down to 80 K. This study extends the data down to 2 K.
Moreover, from single crystal measurements, we previously identified a 3c
layer stacking as well, that commences in the ICDW state and continues to grow
through the NCCDW to CCDW crossover Philip  _et al._ (2023b) It is possible
that both the 3c and 13c coexist at low temperatures in 1T-TaS2. Our single
crystal data only reached 150 K. With doping, the neutron diffraction data
clearly indicate that the 13c structure is suppressed. Its signature
diffraction peak around 0.6 Å is not observed with doping. At the same time,
it is not clear what happens to the 3c stacking with doping. Further
experiments using single crystals are underway to elucidate the doping
dependence of the 3c order.

##  II Results and Discussion

Shown in Fig.1(a) is the hexagonal crystal structure of the high temperature
undistorted lattice. Layers of the transition metal are separated by the
chalcogen ion creating a quasi-2D lattice where weak interlayer interactions
are expected due to the van der Waals nature of the forces holding the layers
together. However, orbitals play an important role in this TMD and out of
plane electron correlations lead to layer ordering and a gap in the density of
states. The out-of-plane coupling is important to understand the electronic
characteristics of these materials where band structure calculations suggested
that opening a gap at the $\Gamma$ point depends on the orbital order and out-
of-plane stacking Ritschel  _et al._ (2015, 2018b). Also shown in Fig.1(a) is
the low temperature crystal structure in the CCDW phase where the high
temperature cell has undergone a $\sqrt{13}\cdot\sqrt{13}$ structural
expansion and a rotation of 13.9o relative to the primary axis. The star
formation is the result of Ta displacements towards the middle Ta ion. The
extent to which the star lattice spreads in the ab-plane depends on
temperature. The star clusters expand on cooling giving rise to large domains
in the CCDW state that are highly ordered, but become disordered on warming,
breaking up into domains with star formations separated by regions of
undistorted lattice. Three samples were measured using neutron scattering and
the diffraction data are shown in Figs. 1(c) and 1(d). Powder samples of TaS2,
TaSSe and TaSe2 were measured as a function of temperature. At 300 K, several
superlattice reflections are indicated that might be due to 3c stacking
structure. Similarly, at 2 K, the same superlattice reflections are observed
as indicated. However, the intensity of the peaks is too small to be discerned
from the powder diffraction data and single crystal experiments will help
elucidate their presence.

The reciprocal lattice vector $\mathbf{Q}$ was calculated using
$\mathbf{Q}=h\mathbf{a_{0}^{*}}+k\mathbf{b_{0}^{*}}+l\mathbf{c^{*}}+m_{1}\mathbf{q^{1}}+m_{2}\mathbf{q^{2}}$,
following the formalism introduced in Spijkerman  _et al._ (1997). The
modulation wave vectors
$\mathbf{q^{1}}=\sigma_{1}\mathbf{a_{0}^{*}}+\sigma_{2}\mathbf{b_{0}^{*}}$ and
$\mathbf{q^{2}}=-\sigma_{2}\mathbf{a_{0}^{*}}+(\sigma_{1}+\sigma_{2})\mathbf{b_{0}^{*}}$
for the CCDW phase were obtained based on the commensurate wave vector
parameters $\sigma_{1}$ and $\sigma_{2}$. The reciprocal lattice commensurate
wave vector can be described as
$\mathbf{q}_{ccdw}=(\sigma_{1}\mathbf{a_{0}^{*}}+\sigma_{2}\mathbf{b_{0}^{*}})$
and for $\sqrt{13}\cdot\sqrt{13}$ in-plane translation, $\sigma_{1}$ and
$\sigma_{2}$ values are 0.2308 and 0.0769 respectively. In the NCCDW phase,
$\sigma_{1}$ and $\sigma_{2}$ values are 0.2448 and 0.0681 respectively. Shown
in Fig. 1(e) is a plot of the diffraction pattern at very small momentum
transfers, Q. At 5 K, a superlattice reflection belonging to 13c ordering is
observed in 1T-TaS2 at low temperatures. This reflection only appears from the
13c stacking order and not from the 3c order, even though most of the higher
order reflections overlap between the two stacking models as shown in Fig.
1(f). The calculated positions of the satellite peaks corresponding to $3\ c$
stacking order ($\mathbf{c^{*}}=\mathbf{c_{0}^{*}}/3$) and $13\ c$ stacking
order ($\mathbf{c^{*}}=\mathbf{c_{0}^{*}}/13$) are shown. Also shown in Fig.
1(e) are data for TaSSe in the same region of momentum transfer. In the TaSSe
data, the reflection $\sim$ 0.6 Å-1 is notably absent which indicates that
there is no 13c ordering in the superconducting state. A similar measurement
was carried out for 1T-TaSe2 and even though the sample was not a single phase
of 1T, no evidence for the 0.6 Å-1 was observed. This indicates that stacking
order might only be present in 1T-TaS2.

Figure 1: (a) The high-temperature and low temperature crystal structure of
1T-TaX2. The lattice symmetry is the trigonal $P\overline{3}m1$ at high
temperature which becomes $P\overline{3}$ at low temperatures. The star is the
result of the Ta displacements. (b) A plot of the diffraction pattern at low
temperatures showing the presence of the 13c superlattice in 1T-TaS2 although
absent in 1T-TaSSe. The calculated peak positions of 1T-TaSSe corresponding to
3c and 13c stacking order is shown at the bottom. (c) The neutron powder
diffraction data collected at 300 K compared among 1T-TaS2, 1T-TaSSe and
1T-TaSe2. All data are fit well using the $P\overline{3}$ symmetry. The
diffraction peaks shift to the left with doping because Se is nominally a
larger ion than S. (d) The neutron powder diffraction data collected at 2 K
are shown. The arrows mark the positions of the CDW superlattice
reflections.(e) The diffraction data plotted at very low Q indicate a
superlattice peak corresponding to the 13c stacking order present in 1T-TaS2.
(f) A plot of the expected positions of 13c and 3c stacking order.

The temperature and composition dependence of the Ta and S/Se thermal factors,
$\langle U\rangle^{2}$, lattice constants and unit cell volume are plotted in
Fig. 2. As a function of composition, superconducting TaSSe has the largest
thermal factor for the Ta ion that continues to increase on warming. Shown in
Fig. 2(c) are the thermal factors for S and Se. Fig. 2(b) is a plot of the c/a
ratio. In the case of TaSSe and TaSe2, the c/a ratio is almost constant as a
function of temperature which indicates that the unit cell expands uniformly
in the a- and c-direction. However, in TaS2, the ratio drops between 200 and
300 K because of the contraction of the c-lattice constant as previously
observed in our earlier study Philip  _et al._ (2023a) and by others Petkov
_et al._ (2022). The contraction of the c-axis corresponds to the transition
from the commensurate to the nearly commensurate state. We observed that this
transition is coupled to the disappearance of the 13c ordering. Shown in Fig.
2(d) is the unit cell volume for the three compositions as a function of
temperature.

Figure 2: (a) The Ta atomic displacement ¡U¿2 is shown as a function of
temperature for the three compositions. Of the three, superconducting TaSSe
shows the largest thermal factors. Shown in (c) are the S and Se thermal
factors for the three compositions. In (b) is a plot of the c/a ratio for the
three compositions and in (d) is a plot of the unit cell volume. The c/a ratio
in 1T-TaS2 shows a decline on warming past 200K.

Fig. 3 shows the results from the pair density function (PDF) analysis of the
Ta displacements from the local structure at 2 K. The local structure is
obtained by Fourier transforming the diffraction data shown in Fig. 1, to
obtain the pair correlation function, G(r). Fitting of the G(r) with a local
model results in the distortions shown in the table of Fig. 3(a). Local Ta
distortions are listed for the 12 Ta ions shown in 3(b). The 13th center Ta
ion does not move by symmetry. This indicates that even after the transition
from the $P\overline{3}m1$ to the $P\overline{3}$ symmetry, locally the stars
are distorted due to displacements of Ta in the directions shown with the
arrows in the star lattice on the right. Moreover, the Ta ions are not all
displaced in a symmetric way. This implies that the local trigonal symmetry is
broken but that there is, nonetheless, long-range order of the star of David
motifs in-plane. Similar distortions were observed in all three compositions
with the results listed for 1T-TaS2.

Figure 3: (a) A list of the Ta distortions obtained from fitting the local
atomic structure. The Ta atoms make up the star of David. The center Ta ion
does not move by symmetry.

##  III Materials and Methods

Powders were prepared using solid-state reaction. The neutron powder
diffraction measurements were performed to investigate the structure through
the multiple CDW steps. The time-of-flight (TOF) neutron measurements were
carried out at the Nanoscale Ordered Materials Diffractometer (NOMAD/BL-1B)
and at SEQUOIA (BL-17), a direct geometry spectrometer, at the Spallation
Neutron Source (SNS) of Oak Ridge National Laboratory (ORNL) at temperatures
ranging from 1.8 to 500 K. The aluminium can was used for SEQUOIA measurements
and the empty can data were subtracted from the data. The reason SEQUOIA was
used is that it reaches very small momentum transfers, not accessible to
NOMAD. The diffraction data from NOMAD were analyzed using the Rietveld
refinement to obtain the unit cell parameters characterizing the crystal
structure Toby (2001), resulting in what is referred to as the average model.
The pair density function (PDF) analysis Proffen  _et al._ (2003); Egami and
Billinge (2003) provides information on the local arrangement of atoms in real
space without the assumption of periodicity. It was performed on the same
neutron diffraction data as the ones used for the Rietveld refinement. NOMAD
is a diffractometer with a large bandwidth of momentum transfer $Q$, and it
provides the total structure function $S(Q)$. The $S(Q)$ was Fourier
transformed into real-space to obtain the $G(r)$ Warren (1990); Peterson  _et
al._ (2000). The instrument background and empty sample container were
subtracted from the $S(Q)$ and the data were normalized by a vanadium rod. A
maximum $Q$ of 40 Å-1 was used.

##  IV Conclusions

The formation of hetero Layer stacking can be engineered to enable new
behaviors and new properties. For instance it has been theoretically proposed
that stacking of the honeycomb ferromagnet CrI3 has the potential to give rise
to ferroelectricity Ji  _et al._ (2023). Moreover, stacking in moire
superlattices can create polar domains because of local spontaneous
polarization Bennett and Remez (2022). Hexagonal boron nitride was shown to
exhibit ferroelectric switching in bilayers, leading to new concepts for
functional heterostructures Yasuda  _et al._ (2021). Similarly, ferromagnetic
heterostructures were demonstrated by stacking non-magnetic WS2 with
antiferromagnetic FePS3. At the interface, the FePS3 shows ferromagnetism Gong
_et al._ (2023).

Layer stacking in homostructure TMDs maybe similalry linked to the transport
behavior. 1T-TaS2 has an insulating CDW in contrast to other CDW
dichalcogenides that have metallic CDW’s. The reason for this is linked to
Mott-Hubbard electron-electron correlations. Every star of David contributes
one 5d electron to a half filled narrow conduction band. In the 13c layer
stacking, there is an odd number of electrons and in the presence of large
Coulomb repulsion acting on the layers, the Mott-Hubbard transition occurs Fei
_et al._ (2022). Density functional theory (DFT) calculations Lee  _et al._
(2019) have shown that the insulating phase and the MI transition originate
not from the 2D order of the stars of David but by the vertical order. Our
results confirm the significance of interlayer coupling and the insulating
property. Interlayer stacking order in the CCDW phase has been verified to be
a 13c repeat unit cell. This result contradicts the notion that the stacking
is partially disordered in the CCDW state. Hence it is less likely that
Anderson localization drives the MI transition. This result also contradicts
the bilayer stacking model.

##  V Acknowledgements

A portion of this research used resources at the Spallation Neutron Source, a
DOE Office of Science User Facility operated by Oak Ridge National Laboratory.
We thank Dr. John Schneeloch (University of Virginia) for valuable inputs on
the sample growth and Dr. Utpal Chatterjee for valuable discussions about the
ARPES data.

## References

  * Wilson  _et al._ (1975) J. A. Wilson, F. Di Salvo, and S. Mahajan, Advances in Physics 24, 117 (1975).
  * Sun (2013) Y. L. R. A. W. J. L. W. H. S. L. J. L. Y. P. Sun, Applied Physics Letters 102 (2013).
  * Liu  _et al._ (2016) Y. Liu, D.-F. Shao, L. Li, W. Lu, X. Zhu, P. Tong, R. Xiao, L. Ling, C. Xi, L. Pi, _et al._ , Physical Review B 94, 045131 (2016).
  * Wang  _et al._ (2020) Y. Wang, W. Yao, Z. Xin, T. Han, Z. Wang, L. Chen, C. Cai, Y. Li, and Y. Zhang, Nature communications 11, 1 (2020).
  * Sipos  _et al._ (2008) B. Sipos, A. F. Kusmartseva, A. Akrap, H. Berger, L. ForrÃ³, and E. TutÇ, Nature Materials 2008 7:12 7, 960 (2008).
  * Law and Lee (2017) K. T. Law and P. A. Lee, Proceedings of the National Academy of Sciences 114, 6996 (2017).
  * Rossnagel (2011) K. Rossnagel, Journal of Physics: Condensed Matter 23, 213001 (2011).
  * Chen  _et al._ (2023) V. Chen, H. R. Lee, C. Köroğlu, C. J. McClellan, A. Daus, and E. Pop, Nano Letters 23, 4095 (2023).
  * Sun  _et al._ (2018) K. Sun, S. Sun, C. Zhu, H. Tian, H. Yang, and J. Li, Science advances 4, eaas9660 (2018).
  * C. J. Butler (2020) T. H. . Y. I. C. J. Butler, M. Yoshida, Nature Communications 11 (2020).
  * Zhang and Wu (2023) W. Zhang and J. Wu, Materials Research Express 10, 046302 (2023).
  * Warawa  _et al._ (2023) K. Warawa, N. Christophel, S. Sobolev, J. Demsar, H. G. Roskos, and M. D. Thomson, arXiv preprint arXiv:2303.08558 (2023).
  * Hansen  _et al._ (2023) M. O. Hansen, Y. Palan, V. Hahn, M. D. Thomson, K. Warawa, H. G. Roskos, J. Demsar, F. Pientka, O. Tsyplyatyev, and P. Kopietz, arXiv preprint arXiv:2303.08580 (2023).
  * Manzke  _et al._ (1989) R. Manzke, T. Buslaps, B. Pfalzgraf, M. Skibowksi, and O. Anderson, Europhys. Lett. 8, 195 (1989).
  * Wegner  _et al._ (2020) A. Wegner, J. Zhao, J. Li, J. Yang, A. A. Anikin, G. Karapetrov, K. Esfarjani, D. Louca, and U. Chatterjee, Physical Review B 101, 195145 (2020).
  * Chatterjee  _et al._ (2015) U. Chatterjee, J. Zhao, M. Iavarone, R. Di Capua, J. Castellan, G. Karapetrov, C. Malliakas, M. G. Kanatzidis, H. Claus, J. Ruff, _et al._ , Nature communications 6, 6313 (2015).
  * van Wezel  _et al._ (2010) J. van Wezel, P. Nahai-Williamson, and S. S. Saxena, Europhysics Letters 89, 47004 (2010).
  * Ang  _et al._ (2015) R. Ang, Z. C. Wang, C. L. Chen, J. Tang, N. Liu, Y. Liu, W. J. Lu, Y. P. Sun, T. Mori, and Y. Ikuhara, Nature Communications 2015 6:1 6, 1 (2015).
  * Spijkerman  _et al._ (1997) A. Spijkerman, J. L. de Boer, A. Meetsma, G. A. Wiegers, and S. van Smaalen, Physical review B 56, 13757 (1997).
  * Fazekas and Tosatti (1979) P. Fazekas and E. Tosatti, Philosophical Magazine B 39, 229 (1979).
  * Wilson  _et al._ (1974) J. A. Wilson, F. J. Di Salvo, and S. Mahajan, Phys. Rev. Lett. 32, 882 (1974).
  * Canfield  _et al._ (1998) P. C. Canfield, P. L. Gammel, and D. J. Bishop, Physics Today 51, 40 (1998), https://pubs.aip.org/physicstoday/article-pdf/51/10/40/8312569/40_1_online.pdf .
  * Kvashnin  _et al._ (2020) Y. Kvashnin, D. VanGennep, M. Mito, S. A. Medvedev, R. Thiyagarajan, O. Karis, A. N. Vasiliev, O. Eriksson, and M. Abdel-Hafiez, Phys. Rev. Lett. 125, 186401 (2020).
  * Egami  _et al._ (2000) T. Egami, Y. Petrov, and D. Louca, J. Superconductivity 13, 709 (2000).
  * Lee  _et al._ (2014) J. Lee, S. Demura, M. B. Stone, K. Iida, G. Ehlers, C. R. Dela Cruz, M. Matsuda, K. Deguchi, Y. Takano, Y. Mizuguchi, _et al._ , Physical Review B 90, 224410 (2014).
  * Smith  _et al._ (1985) N. V. Smith, S. D. Kevan, and F. J. Disalvo, J. Phys. C: Solid State Phys 18, 3175 (1985).
  * Aiura  _et al._ (2004) Y. Aiura, I. Hase, K. Yagi-Watanabe, H. Bando, K. Ozawa, K. Tanaka, R. Kitagawa, S. Maruyama, T. Iwase, Y. Nishihara, K. Horiba, O. Shiino, M. Oshima, M. Nakatake, M. Kubota, and K. Ono, Phys. Rev. B 69, 245123 (2004).
  * Zwick  _et al._ (1998) F. Zwick, H. Berger, I. Vobornik, G. Margaritondo, L. Forró, C. Beeli, M. Onellion, G. Panaccione, A. Taleb-Ibrahimi, and M. Grioni, Phys. Rev. Lett. 81, 1058 (1998).
  * Fei  _et al._ (2022) Y. Fei, Z. Wu, W. Zhang, and Y. Yin, AAPPS Bulletin 32, 20 (2022).
  * Bovet  _et al._ (2004) M. Bovet, D. Popović, F. Clerc, C. Koitzsch, U. Probst, E. Bucher, H. Berger, D. Naumović, and P. Aebi, Phys. Rev. B 69, 125117 (2004).
  * Ritschel  _et al._ (2015) T. Ritschel, J. Trinckauf, K. Koepernik, B. Büchner, M. v. Zimmermann, H. Berger, Y. Joe, P. Abbamonte, and J. Geck, Nature physics 11, 328 (2015).
  * Ritschel  _et al._ (2018a) T. Ritschel, H. Berger, and J. Geck, Physical Review B 98 (2018a), 10.1103/PHYSREVB.98.195134.
  * Lee  _et al._ (2019) S. H. Lee, J. S. Goh, and D. Cho, Physical Review Letters 122, 106404 (2019).
  * Philip  _et al._ (2023a) S. S. Philip, J. C. Neuefeind, M. B. Stone, and D. Louca, Physical Review B 107, 184109 (2023a).
  * Ritschel  _et al._ (2018b) T. Ritschel, H. Berger, and J. Geck, Physical Review B 98, 195134 (2018b).
  * Scruby  _et al._ (1975) C. B. Scruby, P. M. Williams, and G. Parry, Philosophical Magazine 31, 255 (1975).
  * Philip  _et al._ (2023b) S. S. Philip, J. C. Neuefeind, M. B. Stone, and D. Louca, Phys. Rev. B 107, 184109 (2023b).
  * Petkov  _et al._ (2022) V. Petkov, J. Peralta, B. Aoun, and Y. Ren, Journal of Physics: Condensed Matter 34, 345401 (2022).
  * Toby (2001) B. H. Toby, Journal of applied crystallography 34, 210 (2001).
  * Proffen  _et al._ (2003) T. Proffen, S. Billinge, T. Egami, and D. Louca, Zeitschrift für Kristallographie-Crystalline Materials 218, 132 (2003).
  * Egami and Billinge (2003) T. Egami and S. J. Billinge,  _Underneath the Bragg peaks: structural analysis of complex materials_ (Elsevier, 2003).
  * Warren (1990) B. E. Warren,  _X-ray Diffraction_ (Courier Corporation, 1990).
  * Peterson  _et al._ (2000) P. Peterson, M. Gutmann, T. Proffen, and S. Billinge, Journal of Applied Crystallography 33, 1192 (2000).
  * Ji  _et al._ (2023) J. Ji, G. Yu, C. Xu, and H. Xiang, Physical Review Letters 130, 146801 (2023).
  * Bennett and Remez (2022) D. Bennett and B. Remez, npj 2D Materials and Applications 6, 7 (2022).
  * Yasuda  _et al._ (2021) K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Science 372, 1458 (2021).
  * Gong  _et al._ (2023) C. Gong, P. Zhang, T. Norden, Q. Li, Z. Guo, A. Chaturvedi, A. Najafi, S. Lan, X. Liu, Y. Wang, _et al._ , Nature Communications 14, 3839 (2023).