††thanks: A footnote to the article title

# Algorithm to Produce a Density Field with Given Two,Three, and Four-Point
Correlation Functions

Zachary Slepian  Department of Astronomy, University of Florida, Gainesville,
FL 32601

###### Abstract

Here we show how to produce a 3D density field with a given set of higher-
order correlation functions. Our algorithm enables producing any desired two-
point, three-point, and four-point functions, including odd-parity for the
latter. We also generalize the algorithm to i) N-point correlations with
$N>4$, ii) dimensions other than 3, and iii) beyond scalar quantities. This
algorithm should find use in verifying analysis pipelines for higher-order
statistics in upcoming galaxy redshift surveys such as DESI, Euclid, Roman,
and Spherex, as well as intensity mapping.

††preprint: APS/123-QED

##  I  Introduction

For a long time, it has been known how to realize a density field with a given
pair correlation, either power spectrum or 2PCF [1]. This has been useful in
3D large-scale structure (LSS) analyses as well as in CMB analyses. In the CMB
domain, it has also been known how to realize a given bispectrum (triple
correlator of CMB temperature anisotropies) [2]. Furthermore, for LSS, it has
been known how to realize certain specific bispectra, such as those produced
by foorms of primordial non-Gaussianity during inflation.

Yet, a general solution to the problem of realizing arbitrary NPCFs for a 3D
density field has yet to be presented. This, if done, would be quite useful in
testing higher-order correlation analysis pipelines, as is indeed desirable
given the wealth of upcoming 3D LSS data.

It is not quite so problematic to realize a density with a given 2PCF and
3PCF, as this can simply be done using N-body simulations. However, even here,
the simulations are expensive, and each gives only a realization of the true
desired NPCF. Furthermore, non-standard physics, such as parity violation in
the 4PCF, is extremely expensive to simulate, and since it is very model-
dependent, one must simulate many models, each at high expense.

On the observational side, surveys such as DESI have a density-dependent
fiber-assignment effect. This could adversely effect studies of higher-order
correlations. While some of this can be solved by running fiber assignment on
N-body simulations, the currently-standard scheme for fixing fiber assignment
errors (probability inverse pair weighting) may require a prohibitively large
number of simulations if it is to be extended to the 3PCF and beyond. Thus,
producing simulations with a given 3PCF, etc. without the high expense of
performing N-body runs would perhaps bring within reach the number of
realizations needed for this.

In this work, we present a method to realize any desired set of NPCFs about a
given primary galaxy at $\vec{x}$. This may be done about many primaries so
that the spheres of galaxies around each out to the maximal radius to which
one wishes to measure the NPCF do not overlap. One could then simply run an
NPCF algorithm such as [3] by flagging only the primaries desired as those to
be used. This would be a suitable setup for testing the effects outlined
above. The problem of optimal sphere-packing in 3D has also been solved, so
one could optimally pack as many primaries as possible into a given volume.

This work is structured as follows. In Sec. II we review the CMB bispectrum
approach; in Sec. III we extend it to the 3D LSS 4PCF; in Sec. IV to the 3PCF,
and in Sec. V we include 2PCF. We discuss generalization to higher-order
correlators in Sec. VI, higher dimensions of space and arbitrary spatial
geometry in Sec. VII, and conclude in Sec. VIII.

##  II  CMB Maps with a Given Power Spectrum and Bispectrum

###  II.1 CMB correlators

In the CMB context, one decomposes the temperature fluctuations around the
mean as

$\displaystyle\frac{\Delta(T)}{\bar{T}}=\sum_{\ell m}a_{\ell m}Y_{\ell
m}(\hat{n}),$ (1)

with $\Delta T=T(\hat{n})-\bar{T}$, $\hat{n}$ the line of sight, and $\bar{T}$
the mean CMB temperature. One may then compute correlators of the $a_{\ell
m}$. The CMB fluctuations have statistically isotropic correlations; that is,
the pair correlation, triple correlations, etc., are the same, on average, in
any direction one looks, depending only on the relative direction vectors
between the pair or triple, etc. of CMB points. This statistically isotropy
means that the correlators must be invariant under rotations about the
observer.

We have the CMB power spectrum

$\displaystyle\left<a_{\ell m}a_{\ell m}^{*}\right>=C_{\ell}.$ (2)

There are only correlations at equal $\ell$ because of rotation-invariance
about the observer; unequal $\ell$ correlators vanish under rotation-averaging
due to orthogonality of the spherical harmonics. The power spectrum must be
$m$-independent since $m$ is the $z$-component of $\ell$, and hence dependent
on choice of $z$ axis, which would not be rotation-invariant.

The CMB bispectrum is

$\displaystyle\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}a_{\ell_{3}m_{3}}\right>=B_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right),$ (5)

[4, 5] where the matrix is a Wigner 3-$j$ symbol. Given rotation-invariance
about the observer, the Wigner-Eckart theorem implies that the CMB bispectrum
must be proportional to a 3-$j$ symbol.

One conventionally defines the reduced CMB bispectrum via

$\displaystyle
B_{\ell_{1}\ell_{2}\ell_{3}}=\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{0}&{0}&{0}\\\ \end{array}\right)b_{\ell_{1}\ell_{2}\ell_{3}},$ (8)

with

$\displaystyle\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}\equiv\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}}.$
(9)

Combining Eqs. 5 and 8 shows that

$\displaystyle\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}a_{\ell_{3}m_{3}}\right>=\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}}^{m_{1}m_{2}m_{3}}b_{\ell_{1}\ell_{2}\ell_{3}},$
(10)

with the Gaunt integral

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}}^{m_{1}m_{2}m_{3}}\equiv\int
d\Omega\,Y_{\ell_{1}m_{1}}(\hat{n})Y_{\ell_{2}m_{2}}(\hat{n})Y_{\ell_{3}m_{3}}(\hat{n})$
$\displaystyle=\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{0}&{0}&{0}\\\
\end{array}\right)\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right).$ (15)

Reading Eq. 10 and Eq. 15 together highlights two symmetries: the 3-$j$ symbol
with $m$’s, as previously discussed, enforces rotation invariance about the
observer. The 3-$j$ symbol with zeros enforces that
$\ell_{1}+\ell_{2}+\ell_{3}$ is even; this means the CMB bispectrum is parity-
even. These symmetries are reviewed in [6] Sec. IIA. and B.

The CMB trispectrum [6]

$\displaystyle\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}a_{\ell_{3}m_{3}}a_{\ell_{4}m_{4}}\right>=\sum_{\ell_{12},m_{12}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{12}}\\\
{m_{1}}&{m_{2}}&{-m_{12}}\\\ \end{array}\right)$ (18)
$\displaystyle\qquad\times\left(\begin{array}[]{ccc}{\ell_{12}}&{\ell_{3}}&{\ell_{4}}\\\
{m_{12}}&{m_{3}}&{m_{4}}\\\
\end{array}\right)Q_{\ell_{1}\ell_{2}(\ell_{12})}^{\ell_{3}\ell_{4}},$ (21)

where relative to [6] we have switched the upper and lower indices of $Q$,
rewritten $LM\to\ell_{12}m_{12}$ to highlight the coupling scheme used, and
subscripted $\ell_{12}$ where [6] would have it appear like an argument of
$Q$. As discussed in [6], there are several choices for how to couple the
“external” $\ell m$’s, which make no physical difference and are related by
angular momentum recoupling with Wigner 6-$j$ symbols (discussion below his
Eq. 16 and in Appendix A). Here we adopt the coupling scheme of [7], which
matches that used in the body of [6]. We refer to $\ell_{12}$ as an
“intermediate” and parentheses will distinguish intermediates throughout.

###  II.2 Smith-Zaldarriaga Algorithm

It has for some years been of interest to produce a CMB map that has a given
pair (power spectrum) and triple (bispectrum) correlation function. The
problem of realizing a given power spectrum is simple; one just draws $a_{\ell
m}$ as an appropriate Gaussian Random Field (GRF). This problem was also
solved in 3D, as appropriate for large-scale structure, with the well-known
constrained realization algoorithm of [1]. For the CMB bispectrum [2] gives an
algorithm, which we recap below.

[2] defines the functional

$\displaystyle
T[a]\equiv\frac{1}{6}\sum_{\ell_{i}m_{i}}B_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (24) $\displaystyle\qquad\times
a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}a_{\ell_{3}m_{3}}$ (25)

(their eq. 55), where $a=\\{a_{\ell_{i}m_{i}}\\}$ is shorthand for the set of
harmonic coefficients appearing on the righthand side. We note that since the
sum is over all $\ell_{i}m_{i}$, a particular given $\ell m$ pair may appear
three times in the sum; we may have $\ell_{1}m_{1}=\ell m$,
$\ell_{2}m_{2}=\ell m$, or $\ell_{3}m_{3}=\ell m$.

They also define the operator

$\displaystyle\nabla_{\ell m}T[a]\equiv\frac{\partial T[a]}{\partial a_{\ell
m}^{*}}=\frac{1}{2}\sum_{\ell_{i}m_{i}}B_{\ell\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell}&{\ell_{2}}&{\ell_{3}}\\\
{m}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (28)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times
a_{\ell_{2}m_{2}}^{*}a_{\ell_{3}m_{3}}^{*}$ (29)

(their eq. 56). The factor of $1/2$ comes from accounting for the fact that in
the sum in Eq. 24, a particular pair $\ell m$ can appear three times, and the
derivative with respect to this $\ell m$ will pick up each of these, yielding
$3\times 1/6=1/2$. We note that their Eq. 56 seems to have a missing
$\sum_{\ell_{i}m_{i}}$ on its righthand side, but this error does not affect
their Eq. 84 because the expectation values they have there
($<a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}a_{\ell_{2}m_{2}}>$ and same for the
pair with subscripts $\ell_{3}m_{3}$ and $\ell_{3}^{\prime}m_{3}^{\prime}$)
would eliminate the sum in any case. We can be sure this sum is missing from
the righthand side because without it, the number of free indices does not
match between the lefthand and righthand sides.

They then show that if one takes a GRF $a_{\ell m}$ and perturbs it by the
gradient above evaluated at the inverse power spectrum acting on the $a_{\ell
m}$, at leading order in the perturbation, the resulting bispectrum will be
$B_{\ell_{1}\ell_{2}\ell_{3}}$. Writing the perturbed $a_{\ell m}$ as

$\displaystyle a_{\ell m}^{\prime}\equiv a_{\ell m}+\frac{1}{3}\nabla_{\ell
m}T[C^{-1}_{\ell^{\prime}}a_{\ell^{\prime}m^{\prime}}],$ (30)

we may compute

$\displaystyle\left<a_{\ell_{1}m_{1}}^{\prime}a_{\ell_{2}m_{2}}^{\prime}a_{\ell_{3}m_{3}}^{\prime}\right>=\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}\frac{1}{3}\nabla_{\ell_{3}m_{3}}T[C^{-1}_{\ell_{3}^{\prime}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}]\right>$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+5\;{\rm perms.}$ (31)

Terms involving three unperturbed $a_{\ell m}$ vanish because $a_{\ell m}$ is
a GRF and so has vanishing bispectrum. we have dropped terms involving two or
three $\nabla_{\ell m}$ as being sub-leading. The five permutations are the
five other possible orderings of $\ell_{i}m_{i}$.

Applying the definition Eq. 29, the first term of Eq. II.2 becomes

$\displaystyle\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}\frac{1}{3}\nabla_{\ell_{3}m_{3}}T[C^{-1}_{\ell_{3}^{\prime}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}]\right>=$
$\displaystyle\qquad\frac{1}{6}\sum_{\ell_{i}^{\prime}}B_{\ell_{1}\ell_{2}^{\prime}\ell_{3}^{\prime}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}^{\prime}}&{\ell_{3}^{\prime}}\\\
{m_{1}}&{m_{2}^{\prime}}&{m_{3}^{\prime}}\\\ \end{array}\right)$ (34)
$\displaystyle\qquad\times\left<C_{\ell_{2}^{\prime}}^{-1}a_{\ell_{2}m_{2}}a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}C_{\ell_{3}^{\prime}}^{-1}a_{\ell_{3}m_{3}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}^{*}\right>.$

Since $a_{\ell m}$ is a GRF, the expectation value may be pulled apart into
two pairwise expectations:

$\displaystyle\left<C_{\ell_{2}^{\prime}}^{-1}a_{\ell_{2}m_{2}}a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}C_{\ell_{3}^{\prime}}^{-1}a_{\ell_{3}m_{3}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}^{*}\right>=C_{\ell_{2}^{\prime}}^{-1}\left<a_{\ell_{2}m_{2}}a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}\right>$
$\displaystyle\times
C_{\ell_{3}^{\prime}}^{-1}\left<a_{\ell_{3}m_{3}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}^{*}\right>=\delta^{\rm
K}_{\ell_{2}\ell_{2}^{\prime}}\delta^{\rm K}_{m_{2}m_{2}^{\prime}}\delta^{\rm
K}_{\ell_{3}\ell_{3}^{\prime}}\delta^{\rm K}_{m_{3}m_{3}^{\prime}},$ (35)

with $\delta^{\rm K}$ the Kronecker delta, unity if its subscripts are equal
and zero otherwise. Inserting this result in Eq. II.2 we have

$\displaystyle\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}\frac{1}{3}\nabla_{\ell_{3}m_{3}}T[C^{-1}_{\ell_{3}^{\prime}}a_{\ell_{3}^{\prime}m_{3}^{\prime}}]\right>=$
$\displaystyle\qquad\frac{1}{6}B_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right).$ (38)

The same process may be applied to the other five permutations in Eq. II.2, so
we obtain

$\displaystyle\left<a_{\ell_{1}m_{1}}^{\prime}a_{\ell_{2}m_{2}}^{\prime}a_{\ell_{3}m_{3}}^{\prime}\right>=B_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right),$ (41)

i.e. Eq. 84 of [2]. Recalling that the CMB power spectrum $C_{\ell}$ used here
was arbitrary (and cancelled out of the CMB bispectrum obtained above), we see
that this algorithm enables realizing a set of CMB $a_{\ell m}$ with any given
power spectrum and bispectrum.

We now seek to generalize this algorithm to 3D LSS. Let us notice that the CMB
bispectrum actually involves four points: the observer, and the three points
on the CMB spherical shell. This geometry is exactly the same as that used in
the spherical-harmonic decomposition NPCF algorithms [8, 9, 3, 10], but in
that case, the “primary” galaxy is set to be at the origin and the
“secondaries” are on spherical shells with different radii. In particular, the
geometry for the CMB bispectrum is the same as that for the 4PCF around a
given primary.

Let us briefly review this. The NPCF algorithms of the series of papers noted
above expands the density on each shell about a primary galaxy (at $\vec{x}$)
as

$\displaystyle\delta(r,\hat{r};\vec{x})=\sum_{\ell m}a_{\ell
m}(r;\vec{x})_{\ell m}(\hat{r}).$ (42)

The NPCFs are then expanded in the isotropic basis of [7]. In particular, the
4PCF $\zeta$ is

$\displaystyle\zeta(r_{1},r_{2},r_{3};\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})=$
$\displaystyle\qquad\left<\hat{\zeta}(r_{1},r_{2},r_{3};\hat{r}_{1},\hat{r}_{2},\hat{r}_{3};\vec{x})\right>_{\vec{x}}$
$\displaystyle=\sum_{\ell_{i}}\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3}).$
(43)

Hat denotes the estimate of $\zeta$ about a specific primary galaxy at
$\vec{x}$, and $\left<\cdots\right>_{\vec{x}}$ the average over all primaries,
and in the third line, we have performed the expansion into the isotropic
basis $\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}$ [7]. The expansion coefficients
$\zeta_{\ell_{1}\ell_{2}\ell_{3}}$ are then

$\displaystyle\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})=\left<\hat{\zeta}_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3};\vec{x})\right>_{\vec{x}}$
$\displaystyle=\sum_{m_{i}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (46)
$\displaystyle\qquad\times\left<a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})a_{\ell_{3}m_{3}}(r_{3};\vec{x})\right>_{\vec{x}}.$
(47)

The similarity of the last line of Eq. II.2 to the CMB bispectrum is notable.
It leads us to ask if the [2] bispectrum algorithm, if we generalize it so
that each $a_{\ell m}$ may lie on a different spherical shell, could offer a
means of realizing an arbitrary 4PCF about a given primary.

##  III  Realizing a Given 4PCF

We first make a comment about permutations. For the CMB, there is no way to
distinguish the $\ell_{i}m_{i}$, so we end up with the sum over permutations
of them. Once we are working on spherical shells of different radii $r_{i}$,
we may let the $r_{i}$ give a unique ordering to the $\ell_{i}$. Indeed, this
is needed to uniquely define the coupling scheme used for the isotropic basis
functions in which the harmonic-based algorithms measure NPCFs, as detailed in
[7].

For the 4PCF, let us choose the convention that we consider only
configurations of three galaxies around a given primary where the relative
distances $r_{i}$ satisfy $r_{1}\leq r_{2}\leq r_{3}$. This then gives unique
meaning to $\ell_{1},\ell_{2}$, etc. $\ell_{1}$ is the $\ell$ associated with
the harmonic decomposition of the pattern of secondaries on the smallest shell
of a set of three, $\ell_{2}$ of the second-smallest, etc. This convention
eliminates the need for permutations over the $\ell_{i}$, rendering the LSS
problem actually simpler than the CMB problem.

By analogy with Sec. II, let us define

$\displaystyle\mathcal{T}[a]\equiv\sum_{\ell_{i}m_{i},r_{i}}\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (50)
$\displaystyle\qquad\qquad\times
a_{\ell_{1}m_{1}}(r_{1})a_{\ell_{2}m_{2}}(r_{2})a_{\ell_{3}m_{3}}(r_{3}).$
(51)

We note that any particular given value of $\ell m$ may appear in three
different positions: $\ell_{1}m_{1}=\ell m$, $\ell_{2}m_{2}=\ell m$, and
$\ell_{3}m_{3}=\ell m$, just as in the CMB algorithm. However, the ordering by
$r_{i}$ means there is no need to consider permutations.

We also define

$\displaystyle\nabla_{\ell
m}\mathcal{T}\equiv\frac{\partial\mathcal{T}[a]}{\partial(a_{\ell
m}^{*}(r))}=\sum_{\ell_{i}m_{i},r_{i}}\zeta_{\ell\ell_{2}\ell_{3}}(r,r_{2},r_{3})\left(\begin{array}[]{ccc}{\ell}&{\ell_{2}}&{\ell_{3}}\\\
{m}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (54)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times
a_{\ell_{2}m_{2}}^{*}(r_{2})a_{\ell_{3}m_{3}}^{*}(r_{3}).$ (55)

The derivative is intended to be for a specific $\ell m$ and shell $r$.

Let us now consider perturbing the $a_{\ell m}$ on each shell (denoted by $r$)
about a primary galaxy at $\vec{x}$ as

$\displaystyle a_{\ell m}^{\prime}(r;\vec{x})=a_{\ell
m}(r;\vec{x})+\frac{1}{3}\nabla_{\ell
m}\mathcal{T}[C_{\ell^{\prime}}^{-1}a_{\ell^{\prime}m^{\prime}}(r;\vec{x})].$
(56)

Again, by assumption, $a_{\ell m}(r)$ is a GRF with power spectrum on that
shell $C_{\ell}(r)$. We note that the $r$ inside
$a_{\ell^{\prime}m^{\prime}}(r)$ is indeed intended to be unprimed.
Importantly, the perturbation on a given shell depends only on that shell’s
$r$ and $\ell m$. Going back to the definition Eq. 54 we see the perturbation
is a sum over all $r_{2},r_{3}$ and $\ell_{2}m_{2},\ell_{3}m_{3}$.

We now compute the leading-order 4PCF coefficient about the primary, that is,
Eq. II.2 prior to averaging over $\vec{x}$. Let us first focus on the last
line of Eq. II.2, and write it prior to averaging over $\vec{x}$. We have,
dropping terms with three unperturbed $a_{\ell m}$ (as they will vanish after
averaging), and dropping sub-leading terms, that

$\displaystyle
a_{\ell_{1}m_{1}}^{\prime}(r_{1};\vec{x})a_{\ell_{2}m_{2}}^{\prime}(r_{1};\vec{x})a_{\ell_{3}m_{3}}^{\prime}(r_{1};\vec{x})=$
$\displaystyle
a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})\frac{1}{3}\nabla_{\ell_{3}m_{3}}\mathcal{T}[C_{\ell_{3}^{\prime}}^{-1}a_{\ell_{3}^{\prime}m_{3}^{\prime}}(r_{3};\vec{x})]$
$\displaystyle+2\;{\rm perms.}$ (57)

Using the definition Eq. 54 we find

$\displaystyle
a_{\ell_{1}m_{1}}^{\prime}(r_{1};\vec{x})a_{\ell_{2}m_{2}}^{\prime}(r_{2};\vec{x})a_{\ell_{3}m_{3}}^{\prime}(r_{3};\vec{x})=$
$\displaystyle
a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})\,\frac{1}{3}\sum_{\ell_{i}^{\prime}m_{i}^{\prime},r_{i}^{\prime}}\zeta_{\ell_{1}^{\prime}\ell_{2}^{\prime}\ell_{3}}(r_{1}^{\prime},r_{2}^{\prime},r_{3})$
$\displaystyle\qquad\times\left(\begin{array}[]{ccc}{\ell_{1}^{\prime}}&{\ell_{2}^{\prime}}&{\ell_{3}}\\\
{m_{1}^{\prime}}&{m_{2}^{\prime}}&{m_{3}}\\\
\end{array}\right)C_{\ell_{1}^{\prime}}^{-1}(r_{1}^{\prime})C_{\ell_{2}^{\prime}}^{-1}(r_{2}^{\prime})$
(60) $\displaystyle\qquad\times
a_{\ell_{1}^{\prime}m_{1}^{\prime}}^{*}(r_{1}^{\prime};\vec{x})a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}(r_{2}^{\prime};\vec{x})+2\;{\rm
perms.}$ (61)

Now, inserting this in Eq. II.2 and taking the expectation value, we find

$\displaystyle\left<a_{\ell_{1}m_{1}}^{\prime}(r_{1};\vec{x})a_{\ell_{2}m_{2}}^{\prime}(r_{2};\vec{x})a_{\ell_{3}m_{3}}^{\prime}(r_{3};\vec{x})\right>_{\vec{x}}=$
$\displaystyle\frac{1}{3}\sum_{\ell_{i}^{\prime}m_{i}^{\prime},r_{i}^{\prime}}\zeta_{\ell_{1}^{\prime}\ell_{2}^{\prime}\ell_{3}}(r_{1}^{\prime},r_{2}^{\prime},r_{3})$
$\displaystyle\qquad\times\left(\begin{array}[]{ccc}{\ell_{1}^{\prime}}&{\ell_{2}^{\prime}}&{\ell_{3}}\\\
{m_{1}^{\prime}}&{m_{2}^{\prime}}&{m_{3}}\\\
\end{array}\right)C_{\ell_{1}^{\prime}}^{-1}(r_{1}^{\prime})C_{\ell_{2}^{\prime}}^{-1}(r_{2}^{\prime})$
(64)
$\displaystyle\qquad\times\left<a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})a_{\ell_{1}^{\prime}m_{1}^{\prime}}^{*}(r_{1}^{\prime};\vec{x})a_{\ell_{2}^{\prime}m_{2}^{\prime}}^{*}(r_{2}^{\prime};\vec{x})\right>_{\vec{x}}$
$\displaystyle\qquad+2\;{\rm perms.}$ (65)

Now, the $a_{\ell m}$ are, by construction, not correlated across shells, so
we must have $r_{1}=r_{1}^{\prime},r_{2}=r_{2}^{\prime}$, and since on each
shell, the field is a GRF, we have Kronecker deltas setting
$\ell_{1}m_{1}=\ell_{1}^{\prime}m_{1}^{\prime}$ and
$\ell_{2}m_{2}=\ell_{2}^{\prime}m_{2}^{\prime}$ as in the CMB case. Thus the
sum over $\ell_{i}^{\prime}m_{i}^{\prime},r_{i}^{\prime}$ drops out and when
the conditions above are satisfied, we simply get
$C_{\ell_{1}}(r_{1})C_{\ell_{2}}(r_{2})$. We thus have

$\displaystyle\left<a_{\ell_{1}m_{1}}^{\prime}(r_{1};\vec{x})a_{\ell_{2}m_{2}}^{\prime}(r_{2};\vec{x})a_{\ell_{3}m_{3}}^{\prime}(r_{3};\vec{x})\right>_{\vec{x}}=$
$\displaystyle\qquad\frac{1}{3}\,\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (68)
$\displaystyle\qquad+2\;{\rm perms.}$
$\displaystyle\qquad=\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right).$ (71)

Inserting this result for the average of the three $a_{\ell m}^{\prime}$ in
Eq. II.2, we have

$\displaystyle\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})\sum_{m_{i}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\
\end{array}\right)\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{3}}\\\
{m_{1}}&{m_{2}}&{m_{3}}\\\ \end{array}\right)$ (76)
$\displaystyle=\zeta_{\ell_{1}\ell_{2}\ell_{3}}(r_{1},r_{2},r_{3})$ (77)

using the orthogonality identity 34.3.18 in [11] for the 3-$j$ symbols.

Thus, we see that producing a GRF $a_{\ell m}(r)$ on each shell (independent
on each shell), and perturbing on each shell as in Eq. 56, will produce the
desired 4PCF, with coefficients $\zeta_{\ell_{1}\ell_{2}\ell_{3}}$ when
expanded in the isotropic basis of [7]. We emphasize that the perturbation on
a given shell $r$ depends only on that shell’s $r$; it is a sum over all the
other shells’ $r_{i}$. Thus, it can give the correct desired 4PCF at all
combinations of $r_{i}$ shells; it is not just tuned to give the right 4PCF at
only one combination.

Finally, we emphasize that no 3-$j$ symbol with all zero $m$ ever appeared;
thus there is no restriction on the sum of the $\ell_{i}$, and so we can
realize a parity-odd 4PCF of the kind considered in [12, 13, 14].

##  IV  Realizing a Given 3PCF

The harmonic-based 3PCF algorithm [8, 9, 10] essentially sits around a primary
galaxy at $\vec{x}$ and computes angular cross power spectra of two different
spherical shells, $r_{1}$, $r_{2}$. Thus, it is exactly like the CMB angular
power spectrum problem with the observer being analogous to the primary
galaxy, and if we generalize to two shells rather than just one. This provokes
wondering whether simply promoting the unperturbed $a_{\ell m}$ of Sec. III to
have a non-vanishing power spectrum across shells would enable realizing a
desired 3PCF, i.e. allowing $<a_{\ell m(r_{1})a_{\ell m}^{*}(r_{2})}>\neq 0$.
However, this will cause the 4PCF algorithm of Sec. III to break down, as
there we needed to assume that the $a_{\ell m}$ are independent on each shell
so that only the desired $r_{i}$ shells get picked out from the sum in Eq. 54.
We might consider removing this sum from the definition 54 and simply defining
that quantity at a desired set of three shells, but then our 4PCF algorithm
would only reproduce the desired 4PCF at one set of shells at a time. So we
need a different strategy.

Let us consider adding a sum of independent GRFs $b^{(j)}$ to the perturbed
$a^{\prime}_{\ell m}$ of Eq. 56. Let us define them so that at fixed $j$, the
GRF $b^{(j)}$ on each shell is simply a scaled version of some global, shell-
independent GRF $\mathcal{B}^{(j)}$, scaled by a universal function of radius
$f^{(j})(r)$:

$\displaystyle b_{\ell m}^{(j)}(r;\vec{x})\equiv\mathcal{B}_{\ell
m}^{(j)}f^{(j)}(r),$ (78)

where the power spectrum of $\mathcal{B}_{\ell m}^{(j)}$ is
$C_{\mathcal{B},\ell}^{(j)}$. We require that the field has mean zero, i.e.
$<b_{\ell m}^{(j)}>=0$.

We have

$\displaystyle a_{\ell m}^{\prime}(r;\vec{x})=a_{\ell
m}(r;\vec{x})+\frac{1}{3}\nabla_{\ell
m}\mathcal{T}[C_{\ell^{\prime}}^{-1}a_{\ell^{\prime}m^{\prime}}(r;\vec{x})]$
$\displaystyle\qquad+\sum_{j}(C_{\mathcal{B},\ell}^{(j)})^{-1/2}b_{\ell
m}^{(j)}(r;\vec{x}).$ (79)

In the isotropic basis of two-argument functions
$\mathcal{P}_{\ell\ell}(\hat{r}_{1},\hat{r}_{2})$ [7] (or, up to a scaling and
phase, the Legendre basis [8, 9]), the 3PCF is

$\displaystyle\zeta(r_{1},r_{2};\hat{r}_{1}\cdot\hat{r}_{2})=\sum_{\ell}\zeta_{\ell}(r_{1},r_{2})\mathcal{P}_{\ell\ell}(\hat{r}_{1},\hat{r}_{2}).$
(80)

The expansion coefficient is

$\displaystyle\zeta_{\ell}(r_{1},r_{2})=\left<\hat{\zeta}_{\ell}(r_{1},r_{2};\vec{x})\right>_{\vec{x}}$
$\displaystyle=\sum_{m}\left<a_{\ell m}(r_{1};\vec{x})a_{\ell
m}^{*}(r_{2};\vec{x})\right>_{\vec{x}}.$ (81)

We now compute

$\displaystyle\left<a^{\prime}_{\ell m}(r_{1};\vec{x})(a^{\prime}_{\ell
m}(r_{2};\vec{x}))^{*}\right>_{\vec{x}}$
$\displaystyle=\left<\sum_{j_{1}j_{2}}b^{(j_{1})}(r_{1};\vec{x})b^{(j_{2})}(r_{2};\vec{x})\right>_{\vec{x}}$
$\displaystyle=\sum_{j_{1}j_{2}}f^{(j_{1})}(r_{1})f^{(j_{2})}(r_{2})(C_{\mathcal{B},\ell}^{(j_{1})})^{-1/2}(C_{\mathcal{B},\ell}^{(j_{2})})^{-1/2}$
$\displaystyle\qquad\times\left<\mathcal{B}^{(j_{1})}(\vec{x})\mathcal{B}^{(j_{2})}(\vec{x})\right>_{\vec{x}}.$
(82)

Several terms vanished under the averaging. First, the unperturbed $a_{\ell
m}$ is by assumption an independent GRF on each shell, so cross-shell
correlators of it vanish:
$<a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}^{*}(r_{2};\vec{x})>=0$ for
$r_{1}\neq r_{2}$. $a_{\ell m}$ is also uncorrelated with $b_{\ell m}$, so any
terms $<ab>$ vanish. At leading order, we can only have terms with one factor
of $\nabla_{\ell m}$ in $a^{\prime}_{\ell m}(r_{1};\vec{x})(a^{\prime}_{\ell
m}(r_{2};\vec{x}))^{*}$. While each $\nabla_{\ell m}$ contains two $a_{\ell
m}$, which can have a non-zero expectation since they are on the same shell,
these factors of $\nabla_{\ell m}$ will also be multiplied by either an
$a_{\ell m}$ or a $b_{\ell m}$. In the former case, we then have an
expectation value of three $a_{\ell m}$, which vanishes because it is an odd
correlator of a GRF. In the latter case, $b_{\ell m}$ is independent of
$a_{\ell m}$ and we recall that $b_{\ell m}$ has zero mean on each shell. Thus
$<b_{\ell m}><\nabla_{\ell m}>=0$.

We now further simplify Eq. 82. By assumption $\mathcal{B}^{(j_{1})}$ and
$\mathcal{B}^{(j_{2})}$ are uncorrelated unless $j_{1}=j_{2}$. Focusing on the
last line of Eq. 82 and setting $j_{1}=j_{2}$ we have

$\displaystyle(C_{\mathcal{B},\ell}^{(j_{1})})^{-1/2}(C_{\mathcal{B},\ell}^{(j_{2})})^{-1/2}\left<\mathcal{B}^{(j_{1})}(\vec{x})\mathcal{B}^{(j_{2})}(\vec{x})\right>_{\vec{x}}=$
$\displaystyle(C_{\mathcal{B},\ell}^{(j_{1})})^{-1/2}(C_{\mathcal{B},\ell}^{(j_{1})})^{-1/2}\left<\mathcal{B}^{(j_{1})}(\vec{x})\mathcal{B}^{(j_{1})}(\vec{x})\right>_{\vec{x}}$
$\displaystyle=(C_{\mathcal{B},\ell}^{(j_{1})})^{-1}C_{\mathcal{B},\ell}^{(j_{1})}=1.$
(83)

Inserting this result in Eq. 82, we find

$\displaystyle\left<a^{\prime}_{\ell m}(r_{1};\vec{x})(a^{\prime}_{\ell
m}(r_{2};\vec{x}))^{*}\right>_{\vec{x}}$
$\displaystyle=\left<\sum_{j_{1}j_{2}}b^{(j_{1})}(r_{1};\vec{x})b^{(j_{2})}(r_{2};\vec{x})\right>_{\vec{x}}$
$\displaystyle=\sum_{j_{1}}f^{(j_{1})}(r_{1})f^{(j_{1})}(r_{2}).$ (84)

Inserting this into Eq. IV, we have

$\displaystyle\zeta_{\ell}(r_{1},r_{2})=\sum_{m}\sum_{j_{1}}f^{(j_{1})}(r_{1})f^{(j_{1})}(r_{2})$
$\displaystyle=(2\ell+1)\sum_{j_{1}}f^{(j_{1})}(r_{1})f^{(j_{1})}(r_{2}).$
(85)

Thus, we see that, if the 3PCF is factorizable, it may be obtained in this way
if we simply set

$\displaystyle
f^{(j_{1})}(r_{1})f^{(j_{1})}(r_{2})=(2\ell+1)^{-1}\zeta_{\ell}(r_{1},r_{2}).$
(86)

We do notice that the functional form of the dependence on $r_{1}$ and $r_{2}$
must be the same. Even for the real-space tree level 3PCF (e.g. [15]), this
gives an issue as even prior to cyclic summing, the dipole term (with respect
to the angle between $r_{1}$ and $r_{2}$ in the triangle) has the form

$\displaystyle\zeta_{{\rm pc},1}(r_{1},r_{2})=$
$\displaystyle\qquad\xi^{[1+]}(r_{1})\xi^{[1-]}(r_{2})+\xi^{[1+]}(r_{2})\xi^{[1-]}(r_{1}),$
(87)

with subscript “pc” denoting “pre-cyclic” and where

$\displaystyle\xi^{[1\pm]}(r)\equiv\int\frac{k^{2}dk}{2\pi^{2}}\;j_{1}(kr)k^{\pm
1}P(k)$ (88)

and $P(k)$ the power spectrum. Moreover, once redshift-space distortions (RSD)
are included, the pre-cyclic bispectrum begins to include not just wave-
numbers $k_{1}$ and $k_{2}$, but also $k_{3}$ [16], leading to
$r_{3}=[r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}(\hat{r}_{1}\cdot\hat{r}_{2})]^{1/2}$
dependence in the pre-cyclic 3PCF [17]. However, this is suppressed by three
orders of magnitude relative to the $r_{1},r_{2}$-only terms, as shown in [17]
(compare Figs. 1 and 2 there). Thus it is not important to model in a
realization of the 3PCF to match that in tree-level standard perturbation
theory.

Nonetheless, cyclically summing the pre-cyclic 3PCF model and then re-
projecting onto the isotropic basis in $\hat{r}_{1},\hat{r}_{2}$ still leads
to non-factorizable behavior. We use a trick similar to that of [2] to address
this. There, [2] writes any non-factorizable CMB bispectrum as an integral
over a factorizable integrand. [2] then notes that if the integral is
performed numerically with quadrature, that integral becomes a sum over the
integrand evaluated at the quadrature points, and each term in that sum is
factorizable. This then gives the form required for their algorithm, where the
CMB bispectrum is represented as a sum of factorizable pieces.

For us, we observe that, if the post-cyclic 3PCF is written in the form Eq.
80, each of the coefficients $\zeta_{\ell}(r_{1},r_{2})$ will be a sum of
terms derived from projecting the pre-cyclic 3PCF at $\vec{r}_{2},\vec{r}_{3}$
and $\vec{r}_{3},\vec{r}_{1}$ back onto the isotropic basis at
$\hat{r}_{1},\hat{r}_{2}$ (e.g. discussed in Sec. 3 of [17]). Each such term
will look like

$\displaystyle\xi^{[L]}(r_{1})\xi^{[L]}(r_{3})\mathcal{P}_{L}(\hat{r}_{1},\hat{r}_{3})=$
$\displaystyle\sum_{J}\zeta_{J}(r_{1},r_{2})\mathcal{P}_{J}(\hat{r}_{1},\hat{r}_{2})$
(89)

with coefficients

$\displaystyle\zeta_{J}(r_{1},r_{2})=\xi^{[L]}(r_{1})\sum_{L_{1}}\mathcal{C}_{LJL_{1}}^{2}\left(\begin{array}[]{ccc}{L}&{J}&{L_{1}}\\\
{0}&{0}&{0}\\\ \end{array}\right)^{2}$ (92)
$\displaystyle\qquad\times\int\frac{k^{2}dk}{2\pi^{2}}\;P(k)j_{L_{1}}(kr_{1})j_{J}(kr_{2})$
(93)

([17] Eq. 30); we recall that $\mathcal{C}_{LJL_{1}}$ is defined in Eq. 9.

If this integral is performed by numerical quadrature, then we will have a sum
over the quadrature points in $k$, and at each such point, the integrand is
factorizable.

The spherical Bessel functions, at a given $r$, are band-limited functions of
$k$, since they are simply the projection of a plane wave along a line that is
not necessarily perpendicular to the wave-fronts; thus their frequencies are
set by that of the underlying plane wave. Put mathematically, at each $r$,
they are just sums of cosines and sines of $kr$ weighted by rational functions
of $kr$. Thus, given some $r_{\rm max}$ representing the maximal distance out
to which we wish to measure our 3PCF, we may set an exact sampling scheme for
them. Standard tricks like the Poisson summation formula may be employed to
accelerate the convergence of the sum by which we numerically evaluate the
integral.

The power spectrum also should be considered when choosing quadrature points;
it is a smooth envelope, with a small additional “rider” of BAO wiggles
beginning at $k\simeq 0.01\;h/{\rm Mpc}$. Since for numerical integrations,
$P(k)$ is always smoothed by a function such as $\exp[-k^{2}\sigma^{2}]$ with
$\sigma\simeq 1\;{\rm Mpc}/h$ to avoid Gibbs’ phenomenon ringing, $P(k)$
effectively has roughly compact support, further reducing the number of
quadrature points needed. Often a logarithmic grid is employed. Typically of
order a hundred or so points is sufficient.

Conveniently, $r_{\rm max}$ we usually employ for 3PCF measurements (e.g. [18,
19, 20] is roughly twice the BAO scale (i.e. $\sim 200\;{\rm Mpc}/h$), meaning
that the frequency of the sBFs will at most be twice or so the BAO wiggles’
frequency. Now, the two sBFs can resonate (consider multiplying out the sines
and cosine and reducing the products to sums using the double-angle formulae),
giving another factor of two. Thus, a sampling four times as fine as that
needed for the BAO should suffice.

Typically a numerical integration might need to range from $k_{\rm min}\simeq
10^{-4}\;h/{\rm Mpc}$ to $k_{\rm max}\sim 3\;h/{\rm Mpc}$, and sample
logarithmically. The BAO in the transfer function has of order 10 peaks from
$0.01\;h/{\rm Mpc}$ to $0.1\;h/{\rm Mpc}$ ([21] Fig. 3); the power spectrum is
the transfer function’s square and so with resonance of the oscillations, can
have twice that. We then need to sample at roughly four times the frequency of
power spectrum peak for the reasons in the prior paragraph. Thus, we should
have roughly 80 points from $0.01\;h/{\rm Mpc}$ to $0.1\;h/{\rm Mpc}$. If we
sample uniformly per dex, and have about 5 dex from $k_{\rm min}$ to $k_{\rm
max}$, this implies 400 quadrature points. Thus, the integral required by Eq.
92 should be well-approximated by a sum over 400 terms, each of which is
factorizable. This is comparable in order of magnitude to the number of points
used for the Planck forecast in [2] (discussion following their Eq. 44), and
likely can be improved upon significantly using e.g. Poisson summation.

###  IV.1  3PCF Realization Does Not Affect 4PCF

We now show that our addition of terms in $b_{\ell m}$ to $a_{\ell
m}^{\prime}$ as in Eq. 79 does not alter our 4PCF realization, and so one can
produce a single density field with both a desired 3PCF and a desired 4PCF.

We simply need to compute Eq. III again, but using Eq. 79 for $a_{\ell
m}^{\prime}$ (and working at leading order in $b_{\ell m}$ and $\nabla_{\ell
m}$, as always). Only correlators involving at least one $b_{\ell m}$ will be
“new” terms, and any with more than one $b_{\ell m}$ we drop as sub-leading
($b_{\ell m}$ is perturbative). For the same reason, we drop any terms with
more than one factor of $\nabla$. We can have terms as $baa$ and $ba\nabla$.
However, the latter is also sub-leading (because both $b$ and $\nabla$ are
perturbative), so it should also be dropped. Thus, we are left just with a
term:

$\displaystyle\left<a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})b_{\ell_{3}m_{3}}(r_{3};\vec{x})\right>_{\vec{x}}=$
$\displaystyle\left<a_{\ell_{1}m_{1}}(r_{1};\vec{x})a_{\ell_{2}m_{2}}(r_{2};\vec{x})\right>_{\vec{x}}\left<b_{\ell_{3}m_{3}}(r_{3};\vec{x})\right>_{\vec{x}}=0$
(94)

since the mean of $b_{\ell m}$ is zero (noted below Eq. 78). Thus, the
perturbation needed to produce a desired 3PCF will not impact our production
of a desired 4PCF.

##  V  Including the 2PCF

We now discuss how to include a desired 2PCF about the primary galaxy at
$\vec{x}$. For a GRF on each shell $a_{\ell m}(r)$, the 2PCF between the
primary and that shell will simply be

$\displaystyle\xi_{0}(r;\vec{x})=\delta(\vec{x})a_{00}(r;\vec{x}),$ (95)

where $\delta(\vec{x})$ is the density fluctuation associated with the primary
and only $a_{00}$ contributes because the 2PCF is an isotropic average.

We recall that the 4PCF is actually independent of the angular power spectrum
$C_{\ell}(r)$ of the $a_{\ell m}(r)$. Thus we may give $a_{00}$ some desired
“global” radial dependence without affecting the 4PCF, as the inverse power
spectrum there will cancel it out. Thus, we may realize any desired 2PCF along
with any 3PCF and 4PCF. We do note that there is an effect on the 3PCF
monopole when $r_{1}=r_{2}$, but fortunately that is, up to scaling, the 2PCF
squared (at tree-level in PT), and so modifying $a_{00}$ as outlined above is
acceptable.

##  VI  Extending to NPCFs, $N>4$

The fundamental requirements of the algorithm of this work are simply that the
NPCFs be estimated around a primary galaxy by decomposing the density field on
spherical shells in spherical harmonics. Thus, the algorithm trivially extends
to NPCFs for $N>4$. We note this is useful because [3] indeed can measure the
5PCF and 6PCF in the isotropic basis of [7].

Here we write down the analog of $\mathcal{T}$ (Eq. 51) required for realizing
a given 5PCF as well; going to $N>5$ is straightforward but the expressions
are lengthy.

For extending to a 5PCF, we would have four $a_{\ell m}$ about the primary.
Denoting the 5PCF analog of $\mathcal{T}$ as $\mathcal{T}^{(5)}$, we have

$\displaystyle\mathcal{T}^{(5)}[a]\equiv\sum_{\ell_{i}m_{i},r_{i}}\sum_{\ell_{12},m_{12}}\left(\begin{array}[]{ccc}{\ell_{1}}&{\ell_{2}}&{\ell_{12}}\\\
{m_{1}}&{m_{2}}&{-m_{12}}\\\ \end{array}\right)$ (98)
$\displaystyle\qquad\times\left(\begin{array}[]{ccc}{\ell_{12}}&{\ell_{3}}&{\ell_{4}}\\\
{m_{12}}&{m_{3}}&{m_{4}}\\\
\end{array}\right)\zeta_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(r_{1},r_{2},r_{3},r_{4})$
(101) $\displaystyle\qquad\qquad\times
a_{\ell_{1}m_{1}}(r_{1})a_{\ell_{2}m_{2}}(r_{2})a_{\ell_{3}m_{3}}(r_{3})a_{\ell_{4}m_{4}}(r_{4}).$

We recall that $\ell_{12}$ is an “intermediate” (discussion in Sec. II and in
[7]).

Defining $\nabla_{\ell m}$ analogously to Eq. 54, we have

$\displaystyle\nabla_{\ell
m}\mathcal{T}^{(5)}\equiv\frac{\partial\mathcal{T}[a]}{\partial(a_{\ell
m}^{*}(r))}=\sum_{\ell_{i}m_{i},r_{i}}\sum_{\ell_{12},m_{12}}\left(\begin{array}[]{ccc}{\ell}&{\ell_{2}}&{\ell_{12}}\\\
{m}&{m_{2}}&{-m_{12}}\\\ \end{array}\right)$ (104)
$\displaystyle\qquad\qquad\times\left(\begin{array}[]{ccc}{\ell_{12}}&{\ell_{3}}&{\ell_{4}}\\\
{m_{12}}&{m_{3}}&{m_{4}}\\\
\end{array}\right)\zeta_{\ell\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(r,r_{2},r_{3},r_{4})$
(107) $\displaystyle\qquad\qquad\times
a_{\ell_{2}m_{2}}^{*}(r_{2})a_{\ell_{3}m_{3}}^{*}(r_{3})a_{\ell_{4}m_{4}}^{*}(r_{4}).$
(108)

The computation of the 5PCF about a primary proceeds analogously to that in
Sec. III. We note that at leading order, adding a term like this (with a pre-
factor of $1/4$ to account for the four different times one factor of
$\nabla_{\ell m}$ can enter the correlation) will not affect the lower-order
NPCFs.

##  VII  Extending to Dimensions Greater than 3 and non-scalar fields

[22] shows how to extend the harmonic-based algorithm to arbitrary spatial
dimension and arbitrary spatial geometry, via hyperspherical harmonics and
generalizations of the Wigner 3-$j$ symbols. Since this algorithm requires
only the same factorization properties as the 3D, flat space version,
extension of the approach in this paper for realizing NPCFs is
straightforward. Indeed, formally, it will look the same, but the 3-$j$
symbols should be rewritten as Clebsch-Gordan (CG) coefficients, in bra-ket
notation; these are then agnostic as to dimension. In 2D, the CG coefficient
becomes a Kronecker delta; in 3D, it is proportional to a 3-$j$ symbol; in 4D,
it involves a Wigner 9-$j$ symbol (discussed below Eq. 17 of [22]). Then the
coupling coefficients may be built up as in Eq. 22 of [22]; these should
replace the 3-$j$ symbols in our expressions for the 3PCF, 4PCF, etc. when
they are expanded in the isotropic basis.

The algorithm here could also straightforwardly be expanded to non-scalar
fields, such as fields with spin, using the spin-weighted spherical harmonics.
Since these have the same factorization properties as the harmonics used here,
the extension will be straightforward.

Finally, we note that optimal sphere packing has also been solved in higher
dimensions, so one can pack as many primaries as possible in a given volume in
these extensions as well as in the 3D case.

##  VIII  Conclusions

We have presented an algorithm that enables realization, about a given primary
galaxy, of a density field with any desired 2PCF, 3PCF, 4PCF, and even 5PCF.
This should be very useful in the context of testing pipelines for upcoming
NPCF analyses in large-scale structure surveys such as DESI, Euclid, Roman,
Spherex, etc.

We note that the 4PCF algorithm in particular would be enabling for testing
the parity-odd analysis pipeline to search for parity-violation using the
galaxy 4PCF, as initially described in [12] and carried out in [13, 14] using
[23] for the covariance matrix.

###  VIII.1 Second-level heading: Formatting

## References

  * Hoffman and Ribak [1991] Y. Hoffman and E. Ribak, Constrained Realizations of Gaussian Fields: A Simple Algorithm, ApJL 380, L5 (1991).
  * Smith and Zaldarriaga [2011] K. M. Smith and M. Zaldarriaga, Algorithms for bispectra: forecasting, optimal analysis and simulation, MNRAS 417, 2 (2011), arXiv:astro-ph/0612571 [astro-ph] .
  * Philcox  _et al._ [2022] O. H. E. Philcox, Z. Slepian, J. Hou, C. Warner, R. N. Cahn, and D. J. Eisenstein, ENCORE: an O (Ng2) estimator for galaxy N-point correlation functions, MNRAS 509, 2457 (2022), arXiv:2105.08722 [astro-ph.IM] .
  * Spergel and Goldberg [1999] D. N. Spergel and D. M. Goldberg, Microwave background bispectrum. I. Basic formalism, PRD 59, 103001 (1999), arXiv:astro-ph/9811252 [astro-ph] .
  * Goldberg and Spergel [1999] D. M. Goldberg and D. N. Spergel, Microwave background bispectrum. II. A probe of the low redshift universe, PRD 59, 103002 (1999), arXiv:astro-ph/9811251 [astro-ph] .
  * Hu [2001] W. Hu, Angular trispectrum of the cosmic microwave background, PRD 64, 083005 (2001), arXiv:astro-ph/0105117 [astro-ph] .
  * Cahn and Slepian [2020] R. N. Cahn and Z. Slepian, Isotropic N-Point Basis Functions and Their Properties, arXiv e-prints , arXiv:2010.14418 (2020), arXiv:2010.14418 [astro-ph.CO] .
  * Slepian and Eisenstein [2015a] Z. Slepian and D. J. Eisenstein, Computing the three-point correlation function of galaxies in O(N2̂) time, MNRAS 454, 4142 (2015a), arXiv:1506.02040 [astro-ph.CO] .
  * Slepian and Eisenstein [2016] Z. Slepian and D. J. Eisenstein, Accelerating the two-point and three-point galaxy correlation functions using Fourier transforms, MNRAS 455, L31 (2016), arXiv:1506.04746 [astro-ph.CO] .
  * Sunseri  _et al._ [2023] J. Sunseri, Z. Slepian, S. Portillo, J. Hou, S. Kahraman, and D. P. Finkbeiner, SARABANDE: 3/4 point correlation functions with fast Fourier transforms, RAS Techniques and Instruments 2, 62 (2023).
  * Olver, F and Clark, C [2010] B. Olver, F, Lozier, D and Clark, C,  _NIST Handbook of Mathematical Functions_ (NIST and Cambridge University Press, 2010).
  * Cahn  _et al._ [2023] R. N. Cahn, Z. Slepian, and J. Hou, Test for Cosmological Parity Violation Using the 3D Distribution of Galaxies, PRL 130, 201002 (2023).
  * Hou  _et al._ [2023] J. Hou, Z. Slepian, and R. N. Cahn, Measurement of parity-odd modes in the large-scale 4-point correlation function of Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey twelfth data release CMASS and LOWZ galaxies, MNRAS 522, 5701 (2023).
  * Philcox [2022] O. H. E. Philcox, Probing parity violation with the four-point correlation function of BOSS galaxies, PRD 106, 063501 (2022), arXiv:2206.04227 [astro-ph.CO] .
  * Slepian and Eisenstein [2015b] Z. Slepian and D. J. Eisenstein, On the signature of the baryon-dark matter relative velocity in the two- and three-point galaxy correlation functions, MNRAS 448, 9 (2015b), arXiv:1411.4052 [astro-ph.CO] .
  * Scoccimarro  _et al._ [1999] R. Scoccimarro, H. M. P. Couchman, and J. A. Frieman, The Bispectrum as a Signature of Gravitational Instability in Redshift Space, Astrophys. J. 517, 531 (1999), arXiv:astro-ph/9808305 [astro-ph] .
  * Slepian and Eisenstein [2017] Z. Slepian and D. J. Eisenstein, Modelling the large-scale redshift-space 3-point correlation function of galaxies, MNRAS 469, 2059 (2017), arXiv:1607.03109 [astro-ph.CO] .
  * Slepian  _et al._ [2017a] Z. Slepian, D. J. Eisenstein, J. R. Brownstein, C.-H. Chuang, H. Gil-Marín, S. Ho, F.-S. Kitaura, W. J. Percival, A. J. Ross, G. Rossi, H.-J. Seo, A. Slosar, and M. Vargas-Magaña, Detection of baryon acoustic oscillation features in the large-scale three-point correlation function of SDSS BOSS DR12 CMASS galaxies, MNRAS 469, 1738 (2017a), arXiv:1607.06097 [astro-ph.CO] .
  * Slepian  _et al._ [2017b] Z. Slepian, D. J. Eisenstein, F. Beutler, C.-H. Chuang, A. J. Cuesta, J. Ge, H. Gil-Marín, S. Ho, F.-S. Kitaura, C. K. McBride, R. C. Nichol, W. J. Percival, S. Rodríguez-Torres, A. J. Ross, R. Scoccimarro, H.-J. Seo, J. Tinker, R. Tojeiro, and M. Vargas-Magaña, The large-scale three-point correlation function of the SDSS BOSS DR12 CMASS galaxies, MNRAS 468, 1070 (2017b), arXiv:1512.02231 [astro-ph.CO] .
  * Moresco  _et al._ [2021] M. Moresco, A. Veropalumbo, F. Marulli, L. Moscardini, and A. Cimatti, C3: Cluster Clustering Cosmology. II. First Detection of the Baryon Acoustic Oscillations Peak in the Three-point Correlation Function of Galaxy Clusters, Astrophys. J. 919, 144 (2021), arXiv:2011.04665 [astro-ph.CO] .
  * Eisenstein and Hu [1998] D. J. Eisenstein and W. Hu, Baryonic Features in the Matter Transfer Function, Astrophys. J. 496, 605 (1998), arXiv:astro-ph/9709112 [astro-ph] .
  * Philcox and Slepian [2022] O. H. E. Philcox and Z. Slepian, Efficient computation of N-point correlation functions in D dimensions, Proceedings of the National Academy of Science 119, e2111366119 (2022), arXiv:2106.10278 [astro-ph.IM] .
  * Hou  _et al._ [2022] J. Hou, R. N. Cahn, O. H. E. Philcox, and Z. Slepian, Analytic Gaussian covariance matrices for galaxy N -point correlation functions, PRD 106, 043515 (2022), arXiv:2108.01714 [astro-ph.CO] .


# Peaceful coexistence of thermal equilibrium and the emergence of time

Tommaso Favalli<EMAIL_ADDRESS>QSTAR, INO-CNR, and LENS, Largo Enrico
Fermi 2, I-50125 Firenze, Italy  Università degli Studi di Napoli “Federico
II”, Via Cinthia 21, I-80126 Napoli, Italy  Augusto Smerzi
<EMAIL_ADDRESS>QSTAR, INO-CNR, and LENS, Largo Enrico Fermi 2, I-50125
Firenze, Italy

###### Abstract

We consider a quantum Universe composed by a small system $S$ and a large
environment. It has been demonstrated that, for the vast majority of randomly
chosen wave-functions of the Universe satisfying a total energy constraint,
the reduced density matrix of the system $S$ is given by the canonical
statistical distribution. Here, through the Page and Wootters mechanism, we
show that time and non-equilibrium dynamics can emerge as a consequence of the
entanglement between the system and the environment present in the (randomly
chosen) global wave-function of the Universe. The clock is provided by the
environment, which ticks the temporal evolution of $S$. The paradox of the
peaceful coexistence of statistical equilibrium and non-equilibrium dynamics
is solved by identifying the trace over the environment degrees of freedom
with the temporal trace over the entire history of the system $S$.

##  I Introduction

We consider a global quantum system, the “Universe”, composed by a small
system $S$ and a large environment $E$. It has been recently demonstrated
canonical ; popescu1  (see also zurek ) that, under a suitable constraint on
the total energy and for the overwhelming majority of randomly chosen pure
wave-functions of the Universe, the reduced density matrix of $S$ is
indistinguishable from a thermal canonical distribution, a property named in
canonical  as Canonical Typicality.

Back in 1983 Don N. Page and William K. Wootters (PaW) suggested that time can
be an emergent property of entanglement between subsystems in a globally
static Universe pagewootters ; wootters , a proposal that has recently
attracted much attention lloydmaccone ; esp1 ; esp2 ; vedral ;
vedraltemperature ; macconeoptempoarrivo ; wigner ; interacting ; simile ;
simile2 ; leonmaccone ; review ; nostro ; nuovo ; timedilation ;
scalarparticles ; dirac ; foti ; asimmetry  as a viable route for a new
description of space-time, including a new perspective for merging quantum
clocks and gravity brukner ; giacomini  (see Appendix A for a brief summary of
PaW theory).

The goal of this work is to show that canonical thermal equilibrium and
dynamics can coexist in a quantum Universe, with the environment providing the
clock for the evolution of the system $S$. The dynamics is governed by the
Schrödinger equation corrected by a non-local term which vanishes in the limit
of fixed total energy of the Universe. The paradox of the coexistence of
thermal equilibrium and non-trivial evolution of $S$ is solved because the
trace over the degrees of freedom of the environment coincides with a temporal
average over the entire life of $S$. The temporal dynamics of the system $S$
emerges by considering the relative states of $S$ (in Everett sense everett )
with respect to the states of the environment.

##  II Canonical Equilibrium Distribution

Our quantum Universe is composed by a small system $S$ weakly interacting with
a large environment $E$ and it is governed by the Hamiltonian

$\hat{H}=\hat{H}_{E}+\hat{H}_{S}$ (1)

where $\hat{H}_{E}$ and $\hat{H}_{S}$ are the Hamiltonians of the subsystems
$E$ and $S$ respectively, having dimensions $d_{E}\gg d_{S}$. The relatively
small interaction between the environment and the system $S$ is neglected
canonical ; popescu1 . In analogy with the standard derivation of the
canonical distribution of a subsystem given a global microcanonical
distribution, we impose the total energy in the interval
$\left[E,E+\delta\right]$, where $\delta\ll E,\Delta E^{(S)}$ (being $\Delta
E^{(S)}$ the typical spacing between energy levels of the system $S$) but
large enough to contain many eigenvalues of $\hat{H}_{E}$. This constraint
corresponds to the choice of
$\mathcal{H}_{U}\subseteq\mathcal{H}_{E}\otimes\mathcal{H}_{S}$. We consider,
following canonical , a Universe in a pure state
$\ket{\Psi}\in\mathcal{H}_{U}$. The state of the system $S$ is obtained after
tracing out the environment degrees of freedom:

$\hat{\rho}_{S}=Tr_{E}\left[\ket{\Psi}\bra{\Psi}\right].$ (2)

According to canonical ; popescu1 , for almost every pure state
$\ket{\Psi}\in\mathcal{H}_{U}$, the state of $S$ is well approximated by the
canonical distribution

$\hat{\rho}_{S}\approx
Tr_{E}\left[\hat{\Omega}_{U}\right]\approx\frac{1}{Z}e^{-\beta\hat{H}_{S}}$
(3)

where $\beta=dS(E)/dE$ is the inverse temperature with $S(E)$ the environment
entropy, $Z=Tr\left[e^{-\beta\hat{H}_{S}}\right]$ and
$\hat{\Omega}_{U}=d_{U}^{-1}\hat{P}_{U}$ is the equiprobable mixed state in
$\mathcal{H}_{U}$ with $d_{U}$ the dimension of the space $\mathcal{H}_{U}$
and $\hat{P}_{U}$ the projection on $\mathcal{H}_{U}$. Notice that equal
probabilities (and random phases) are assigned in this case to all the states
within $\hat{\Omega}_{U}$ which is thus maximally mixed in $\mathcal{H}_{U}$
cantyp . Equation (3) implies that the thermal state of the small subsystem
$S$ can be derived from a (randomly chosen) pure state
$\ket{\Psi}\in\mathcal{H}_{U}$ or from the maximally mixed state
$\hat{\Omega}_{U}$.

##  III Environment as a Clock

###  III.1 General framework

We are now going to merge Canonical Typicality and PaW theory. The key point
is to recognize the environment as a clock:

$\hat{H}_{E}\equiv\hat{H}_{C}.$ (4)

Notice that in the PaW framework a good clock has to have a Hilbert space
dimension larger than the dimension of the system $S$, otherwise it would no
longer be possible to relate each energy eigenstate of the system to an energy
eigenstate of the clock (see Appendix A and nostro ). Furthermore a good clock
has to interact weakly with the system $S$ or, in the ideal case, it should
not interact at all. These conditions coincide with those required for the
environment by Canonical Typicality.

In the original PaW theory the global state of the Universe is an eigenstate
of the total Hamiltonian with zero eigenvalue, i.e. $\hat{H}\ket{\Psi}=0$. In
this work we follow a slightly different path and we weakly relax the PaW
constraint considering the total energy of the Universe within the energy
shell $\left[E,E+\delta\right]$ where $\delta\ll E,\Delta E^{(S)}$ but large
enough to contain many energy eigenvalues of the clock $\equiv$ environment
$C$. In this framework we find a non-local Schrödinger-like evolution for the
relative state of $S$ that reduces to the usual Schrödinger dynamics for times
$t-t_{0}\ll 1/\delta$ (where $t_{0}$ is the initial time) or for all times in
the limit $\delta\rightarrow 0$.

We define the time states in the Hilbert space of the clock $\equiv$
environment using the approach developed in nostro ; pegg . We assume that the
Hamiltonian of the environment has non-degenerate energy eigenstates with
rational energy ratios

$\frac{E^{(C)}_{i}}{E^{(C)}_{1}}=\frac{A_{i}}{D_{i}},$ (5)

where $A_{i}$ and $D_{i}$ are integers with no common factors and
$E^{(C)}_{0}=0$. This implies that all energy values are integer multiples of
a (arbitrarily small) step ($\hslash=1$):

$E^{(C)}_{i}=r_{i}\frac{2\pi}{T}$ (6)

where $T=\frac{2\pi r_{1}}{E_{1}}$, $r_{i}=r_{1}\frac{A_{i}}{D_{i}}$ for $i>1$
(with $r_{0}=0$) and $r_{1}$ is equal to the lowest common multiple of the
values of $D_{i}$. In this space we define the states

$\ket{t_{m}}=\frac{1}{\sqrt{p+1}}\sum_{i=0}^{p}e^{-iE^{(C)}_{i}t_{m}}\ket{E^{(C)}_{i}}$
(7)

with $p+1=d_{C}$, $t_{m}=t_{0}+mT/(s+1)$, $m=0,1,2,...,s$ and $s+1\geq r_{p}$.
The number of states $\ket{t_{m}}$ is therefore greater than the number of
energy states in $\mathcal{H}_{C}$ and the $s+1$ values of $t_{m}$ are
uniformly distributed over $T$. These states are not orthogonal but provide an
overcomplete basis in $C$ with the resolution of the identity

$\frac{p+1}{s+1}\sum_{m=0}^{s}\ket{t_{m}}\bra{t_{m}}=\mathbb{1}_{C}.$ (8)

In order to obtain a continuous flow of time we can now consider the limit
$s\rightarrow\infty$ and define

$\ket{{t}}=\sum_{i=0}^{p}e^{-iE^{(C)}_{i}t}\ket{E^{(C)}_{i}}$ (9)

where $t$ can now take any real value from $t_{0}$ to $t_{0}+T$. In this limit
the resolution of the identity (8) becomes

$\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\ket{t}\bra{t}=\mathbb{1}_{C}$ (10)

and the states $\ket{t}$ provide again an overcomplete basis with
$\braket{t}{t^{\prime}}=\sum_{n}e^{iE^{(C)}_{n}(t-t^{\prime})}$. We notice
that, if we would consider non-rational ratios of energy levels, the
resolutions of the identity (8) and (10) are no longer exact and the time
states $\ket{t_{m}}$ and $\ket{t}$ do not provide an overcomplete basis in
$C$. However, since any real number can be approximated with arbitrary
precision by a ratio between two rational numbers, the residual terms in the
resolutions of the identity can be arbitrarily reduced nostro .

###  III.2 Random Universe and dynamics

Here we show that a Universe in a (randomly chosen) pure state is compatible
with the emergence of time and a non-trivial dynamical evolution of the system
$S$. The global state of the Universe is

$\ket{\Psi}=\sum_{j}\sum_{i}c_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}$ (11)

where we take the coefficients $c_{ij}$ distributed as in canonical . With
$\ket{\Psi}=\ket{\Phi}/||\ket{\Phi}||$ and
$\ket{\Phi}=\sum_{j}\sum_{i}\tilde{c}_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}$,
the real and imaginary parts of the coefficients
$\tilde{c}_{ij}=c_{ij}||\ket{\Phi}||$ are chosen as independent real Gaussian
random variables with mean zero and variance $1/2$ for the values of $i,j$
such that $E^{(C)}_{i}+E^{(S)}_{j}\in\left[E,E+\delta\right]$ and
$\tilde{c}_{ij}=0$ otherwise. Since $\delta\ll\Delta
E^{(S)}_{j}=E^{(S)}_{j+1}-E^{(S)}_{j}~{}\forall j$ and considering that the
spectrum of $\hat{H}_{C}$ is much denser than the spectrum of $\hat{H}_{S}$,
the constraint on the total energy implies that each level of the system is
coupled with several neighbour levels of the clock. The random choice of the
coefficients provides Canonical Typicality canonical : here we show that is
also sufficient to provide the temporal dynamics for the relative state (in
Everett sense everett ) of the system $S$.

The action of the global Hamiltonian $\hat{H}$ on the global state
$\ket{\Psi}$ gives

$\begin{split}\hat{H}\ket{\Psi}&=\left(\hat{H}_{C}+\hat{H}_{S}\right)\sum_{j}\sum_{i\in
I_{j}}c_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}=\\\ &=\sum_{j}\sum_{i\in
I_{j}}c_{ij}\left(E^{(C)}_{i}+E^{(S)}_{j}\right)\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}=\\\
&=E\ket{\Psi}+\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}\end{split}$ (12)

where $I_{j}$ is the set of the environment levels such that
$E^{(C)}_{i}\in\left[E-E^{(S)}_{j},E-E^{(S)}_{j}+\delta\right]$ and where we
have written $E^{(C)}_{i}+E^{(S)}_{j}=E+\Delta_{ij}$ with
$\Delta_{ij}\in\left[0,\delta\right]$. The relative state of the system $S$ at
a certain “time” $t$ is defined by

$\ket{\phi(t)}_{S}=\braket{t}{\Psi}$ (13)

(notice that (13) is still a pure state) and its time evolution can be easily
calculated:

$\begin{split}i\frac{\partial}{\partial
t}\ket{\phi(t)}_{S}&=\frac{\partial}{\partial
t}\sum_{k}^{d_{C}}\bra{E_{k}}e^{iE_{k}t}\ket{\Psi}=-\sum_{k}^{d_{C}}\bra{E_{k}}E_{k}e^{iE_{k}t}\ket{\Psi}=\\\
&=-\bra{t}\hat{H}_{C}\ket{\Psi}=\hat{H}_{S}\braket{t}{\Psi}-\bra{t}\hat{H}\ket{\Psi}=\\\
\\\ &=\left(\hat{H}_{S}-E\right)\ket{\phi(t)}_{S}-\bra{t}\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}\end{split}$ (14)

where we have used (9), $\hat{H}=\hat{H}_{C}+\hat{H}_{S}$, (13) and (12). By
defining the operator $\hat{\Delta}=\sum_{j}\sum_{i\in
I_{j}}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}\bra{E^{(C)}_{i}}\bra{E^{(S)}_{j}}$
and removing the term related to $E$ which gives an irrelevant phase factor in
the evolution of $S$, (14) becomes the time non-local Schrödinger-like
equation:

$i\frac{\partial}{\partial
t}\ket{\phi(t)}_{S}=\hat{H}_{S}\ket{\phi(t)}_{S}-\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt^{\prime}\hat{\Delta}(t,t^{\prime})\ket{\phi(t^{\prime})}_{S}$
(15)

where $\hat{\Delta}(t,t^{\prime})=\bra{t}\hat{\Delta}\ket{t^{\prime}}$. The
last term in the right-hand side of the equation is an integral operator
acting on $S$.

For times $t-t_{0}\ll 1/\delta$ (and so $t-t_{0}\ll 1/\Delta_{ij}$ for typical
$\Delta_{ij}$) (15) reduces to the ordinary Schrödinger equation. Indeed we
have:

$\begin{split}\ket{\phi(t)}_{S}&=\braket{t}{\Psi}=\bra{t_{0}}e^{i\hat{H}_{C}(t-t_{0})}\ket{\Psi}=e^{-i\hat{H}_{S}(t-t_{0})}\bra{t_{0}}e^{i\hat{H}(t-t_{0})}\ket{\Psi}\\\
&=e^{-i\hat{H}_{S}(t-t_{0})}\bra{t_{0}}\sum_{j}\sum_{i\in
I_{j}}c_{ij}e^{i(E^{(C)}_{i}+E^{(C)}_{j})(t-t_{0})}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}=\\\
&=e^{-i\hat{H}_{S}(t-t_{0})}\bra{t_{0}}\sum_{j}\sum_{i\in
I_{j}}c_{ij}e^{iE(t-t_{0})}e^{i\Delta_{ij}(t-t_{0})}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}\end{split}$
(16)

where we used (13), $\hat{H}=\hat{H}_{C}+\hat{H}_{S}$ and
$E^{(C)}_{i}+E^{(S)}_{j}=E+\Delta_{ij}$. For $t-t_{0}\ll 1/\delta$,
considering $e^{i\Delta_{ij}(t-t_{0})}\simeq 1$ and removing the irrelevant
global phase factor $e^{iE(t-t_{0})}$, (16) becomes

$\ket{\phi(t)}_{S}\simeq e^{-i\hat{H}_{S}(t-t_{0})}\ket{\phi(t_{0})}_{S}$ (17)

which provides the Schrödinger evolution for the system $S$. In Section III.D
we’ll briefly discuss the effect of the non-local term in equation (15) for
times $t-t_{0}\geq 1/\delta$.

Equation (15) can be explicitly solved: with
$E^{(C)}_{i}+E^{(S)}_{j}=\Delta_{ij}$ (the term related to $E$ has been
removed) we obtain (see Appendix B):

$\begin{split}\ket{\phi(t)}_{S}=\sum_{j}\alpha_{j}(t)e^{-iE^{(S)}_{j}t}\ket{E^{(S)}_{j}}\end{split}$
(18)

with

$\alpha_{j}(t)=\sum_{i\in I_{j}}c_{ij}e^{i\Delta_{ij}t}.$ (19)

Equation (18) provides an additional time dependence to the Schrödinger
evolution through the coefficients $\alpha_{j}(t)$. In the case $t-t_{0}\ll
1/\delta$ we have $\alpha(t)\simeq\sum_{i\in
I_{j}}c_{ij}e^{i\Delta_{ij}t_{0}}\equiv\alpha_{j}(t_{0})$ and the state (18)
becomes
$\ket{\phi(t)}_{S}\simeq\sum_{j}\alpha_{j}(t_{0})e^{-iE^{(S)}_{j}t}\ket{E^{(S)}_{j}}$
where we recognize again the Schrödinger evolution for the system $S$.

We emphasize that (15) does not preserve the norm of the state
$\ket{\phi(t)}_{S}$ over time. Indeed, calculating the scalar product
$\braket{\phi(t)}{\phi(t)}$ through (18) we obtain

$\begin{split}\braket{\phi(t)}{\phi(t)}&=\sum_{j}\left|\alpha_{j}(t)\right|^{2}=\sum_{j}\sum_{i\in
I_{j}}\sum_{k\in I_{j}}c_{ij}c^{*}_{kj}e^{i(\Delta_{ij}-\Delta_{kj})t}=\\\
&=\sum_{j}\sum_{i\in I_{j}}\sum_{k\in
I_{j}}|c_{ij}||c_{kj}|\cos((\Delta_{ij}-\Delta_{kj})t-\Delta\varphi^{(j)}_{ik})\end{split}$
(20)

where $c_{ij}=|c_{ij}|e^{i\varphi_{ij}}$, $c_{kj}=|c_{kj}|e^{i\varphi_{kj}}$
and $\Delta\varphi^{(j)}_{ik}=\varphi_{kj}-\varphi_{ij}$. However, the
corrections remain small when $t-t_{0}\ll 1/\delta$ (and vanish in the limit
$\delta\rightarrow 0$) being
$\braket{\phi(t)}{\phi(t)}\simeq\sum_{j}\sum_{i,k\in
I_{j}}c_{ij}c^{*}_{kj}e^{i(\Delta_{ij}-\Delta_{kj})t_{0}}=\sum_{j}|\alpha_{j}(t_{0})|^{2}$
that is different from the unity but (approximately) constant over time. We
notice that a similar problem arose, for different reasons, in interacting
where it was handled by introducing a new definition for the inner product.

To summarize, the environment can provide the clock for the dynamical
evolution of the system $S$. The state of $S$ conditioned to a certain value
of the clock through (13) consists of a pure state obeying a non-local
dynamical Schrödinger-like equation (that reduces to the standard Schrödinger
equation for $t-t_{0}\ll 1/\delta$). Nevertheless, after tracing out the
degrees of freedom of the clock $\equiv$ environment, we find the system in a
state of thermal equilibrium provided by the canonical distribution. This
compatibility is simply explained by the fact that the trace over the
environment degrees of freedom is equivalent to the trace over all times.
Indeed we have (see Appendix D):

$\hat{\rho}_{S}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\braket{t}{\Psi}\braket{\Psi}{t}\approx\frac{1}{Z}e^{-\beta\hat{H}_{S}}$
(21)

where again $\beta=dS(E)/dE$ is the inverse temperature with $S(E)$ the clock
$\equiv$ environment entropy and $Z=Tr\left[e^{-\beta\hat{H}_{S}}\right]$ as
in equation (3).

###  III.3 Initial conditions for $S$

The merging of Canonical Typicality and PaW imposes a constraint on the
allowed initial conditions of the subsystem $S$. The initial condition for the
state (18) is

$\ket{\phi(t_{0})}_{S}=\sum_{j}\alpha_{j}(t_{0})e^{-iE^{(S)}_{j}t_{0}}\ket{E^{(S)}_{j}}$
(22)

with

$\alpha_{j}(t_{0})=\sum_{i\in I_{j}}c_{ij}e^{i\Delta_{ij}t_{0}}.$ (23)

The reduced density matrix $\hat{\rho}_{S}$ of the subsystem $S$ is

$\begin{split}\hat{\rho}_{S}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\sum_{j}\sum_{i\in
I_{j}}\left|c_{ij}\right|^{2}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{j}}\end{split}$
(24)

where, following canonical , we have considered that the relevant energy
intervals for the clock $\equiv$ environment coupled with each level
$E^{(S)}_{j}$ of $S$ are pairwise disjoint (which is a consequence of
$\delta\ll\Delta E^{(S)}_{j}=E^{(S)}_{j+1}-E^{(S)}_{j}~{}\forall j$).
Canonical Typicality implies:

$\sum_{i\in I_{j}}\left|c_{ij}\right|^{2}\approx\frac{1}{Z}e^{-\beta
E^{(S)}_{j}}$ (25)

which constrains the initial conditions by selecting a set of allowed
$\alpha_{j}(t_{0})$ through (23). It is crucial to notice that the condition
(25) contrains the sum of the absolute values of the coefficients $c_{ij}$ and
therefore leaves a large margin of freedom on the possible values of
$\alpha_{j}(t_{0})$. In conclusion, Canonical Typicality states that the
reduced density matrix
$\hat{\rho}_{S}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]$ of the overwhelming
majority of the pure wave-functions
$\ket{\Psi}\in\mathcal{H}_{U}\subseteq\mathcal{H}_{C}\otimes\mathcal{H}_{S}$
is canonical. This means that the overwhelming majority of the randomly chosen
coefficients $c_{ij}$ satisfies the condition (25) which, in our framework, is
consistent with a non-trivial dynamical evolution of the system $S$.

###  III.4 A toy model

We look now at a simple example assuming that the subsystem $S$ consists of a
one-dimensional harmonic oscillator with Hamiltonian
$\hat{H}_{S}=\frac{\hat{P}^{2}}{2m}+\frac{1}{2}m\omega^{2}\hat{X}^{2}$. We
assume that the dynamics is confined among the two lowest energy levels of the
oscillator and therefore the global state of the Universe (11), satisfying the
constraint on the total energy, is $\ket{\Psi}=\sum_{i\in
I_{0}}c_{i0}\ket{E^{(C)}_{i}}\ket{0^{(S)}}+\sum_{i\in
I_{1}}c_{i1}\ket{E^{(C)}_{i}}\ket{1^{(S)}}$.

With (22) and $t_{0}=0$, we look at the (pure) initial state of the subsystem
$S$:
$\ket{\phi(0)}_{S}=\alpha_{0}(0)\ket{0^{(S)}}+\alpha_{1}(0)\ket{1^{(S)}}$,
where we set the initial values $\alpha_{0}(0)$ and $\alpha_{1}(0)$ according
to (23). Thanks to (18) and (19), for the conditioned state of $S$ at a
generic time $t$, we have

$\begin{split}\ket{\phi(t)}_{S}&=\alpha_{0}(t)e^{-iE^{(S)}_{0}t}\ket{0^{(S)}}+\alpha_{1}(t)e^{-iE^{(S)}_{1}t}\ket{1^{(S)}}=\\\
&=\sum_{i\in
I_{0}}c_{i0}e^{-i(E^{(S)}_{0}-\Delta_{i0})t}\ket{0^{(S)}}+\sum_{i\in
I_{1}}c_{i1}e^{-i(E^{(S)}_{1}-\Delta_{i1})t}\ket{1^{(S)}}.\end{split}$ (26)

We recall here that the state (26) is not normalized. To restore the
normalization (for $t\ll 1/\delta$ where the norm is approximately preserved)
we should divide (26) by $\sqrt{\alpha^{2}_{0}(0)+\alpha^{2}_{1}(0)}$.
However, for the sake of simplicity, we proceed with the calculation without
considering the normalization.

We look now at the time dependence of the expectation value of the position
operator $\langle\hat{X}\rangle_{t}=\bra{\phi(t)}\hat{X}\ket{\phi(t)}$ and
obtain (see Appendix E):

$\begin{split}\langle\hat{X}\rangle_{t}&=\sqrt{\frac{2}{m\omega}}|\alpha_{0}(t)||\alpha_{1}(t)|\cos(\omega
t-\Delta\phi(t))=\\\ &=\sqrt{\frac{2}{m\omega}}\sum_{i\in I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|\cos((\omega+(\Delta_{i0}-\Delta_{k1}))t-\Delta\varphi^{(0,1)}_{ik})\end{split}$
(27)

where we have considered $\alpha_{0}(t)=|\alpha_{0}(t)|e^{i\phi_{0}(t)}$,
$\alpha_{1}(t)=|\alpha_{1}(t)|e^{i\phi_{1}(t)}$,
$\Delta\phi(t)=\phi_{1}(t)-\phi_{0}(t)$, $c_{i0}=|c_{i0}|e^{i\varphi_{i0}}$,
$c_{k1}=|c_{k1}|e^{i\varphi_{k1}}$,
$\Delta\varphi^{(0,1)}_{ik}=\varphi_{k1}-\varphi_{i0}$ and
$E^{(S)}_{1}-E^{(S)}_{0}=\omega$. For times $t\ll
1/|\Delta_{i0}-\Delta_{k1}|$, up to first order of approximation in
$t(\Delta_{i0}-\Delta_{k1})$, (27) reduces to (see Appendix F):

$\langle\hat{X}\rangle_{t}\simeq\sqrt{\frac{2}{m\omega}}|\alpha_{0}(0)||\alpha_{1}(0)|\cos(\omega
t-\Delta\phi(0))-\sqrt{\frac{2}{m\omega}}\sum_{i\in I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|t(\Delta_{i0}-\Delta_{k1})\sin(\omega
t-\Delta\varphi^{(0,1)}_{ik})$ (28)

where $\alpha_{0}(0)=|\alpha_{0}(0)|e^{i\phi_{0}(0)}$,
$\alpha_{1}(0)=|\alpha_{1}(0)|e^{i\phi_{1}(0)}$,
$\Delta\phi(0)=\phi_{1}(0)-\phi_{0}(0)$ and where we used again (23). Equation
(28) indicates, as expected, that the expectation value
$\langle\hat{X}\rangle_{t}$ oscillates between
$\pm\sqrt{\frac{2}{m\omega}}|\alpha_{0}(0)||\alpha_{1}(0)|$ with frequency
$\omega$ (apart from small corrections) and this is not surprising since we
know that for $t\ll 1/\delta$ the system $S$ exhibits a Scrödinger-like
evolution.

If we trace over the clock $\equiv$ environment degrees of freedom, which
corresponds to a time average, for the overwhelming majority of the randomly
chosen coefficients $c_{ij}$ we obtain, thanks to (21), (24) and (25), the
canonical mixed density matrix for the subsystem $S$:

$\begin{split}\hat{\rho}_{S}&=\frac{1}{T}\int_{0}^{T}dt\braket{t}{\Psi}\braket{\Psi}{t}=\\\
&=\sum_{i\in I_{0}}|c_{i0}|^{2}\ket{0^{(S)}}\bra{0^{(S)}}+\sum_{i\in
I_{1}}|c_{i1}|^{2}\ket{1^{(S)}}\bra{1^{(S)}}\approx\\\
&\approx\frac{1}{Z}\left(e^{-\beta
E^{(S)}_{0}}\ket{0^{(S)}}\bra{0^{(S)}}+e^{-\beta
E^{(S)}_{1}}\ket{1^{(S)}}\bra{1^{(S)}}\right).\end{split}$ (29)

##  IV Discussion

###  IV.1 On conditional probabilities

An important point in the PaW mechanism concerns conditional probabilities. In
the original PaW proposal the probability to obtain the outcome $a$ when
measuring the observable $\hat{A}$ (with $\hat{A}\ket{a}=a\ket{a}$) on the
subspace $S$ at a certain clock time $t$ is:

$\begin{split}p(a\>on\>S\>|\>t\>on\>C)=\frac{p(a\>on\>S,\>t\>on\>C)}{p(t\>on\>C)}\end{split}$
(30)

that is the conditional probability of obtaining $a$ on $S$ given that the
clock $C$ shows $t$.

This aspect of the theory has been criticised by K. V. Kuchar kuchar  who
emphasized that the PaW mechanism is not able to provide the correct
propagators when considering multiple measurements. Indeed measurements of the
system at two times will give the wrong statistics because the first
measurement “collapses”the time state and freezes the system. Two possible
ways out of this problem have been proposed: the first in gppt  (which we call
GPPT theory in the following) with an experimental illustration in esp1  and
the second in lloydmaccone .

We focus now on the GPPT proposal. As pointed out in esp1  one of the main
ingredients in the GPPT theory is the averaging over the abstract coordinate
time (the “external time”) in order to eliminate any external time dependence
in the observables. So, in the GPPT proposal, the probability to obtain the
outcome $a$ when measuring the observable $\hat{A}$ on the subsystem $S$
conditioned to having the outcome $t$ on $C$ is given by gppt :

$p(a\>on\>S\>|\>t\>on\>C)=\frac{\int
d\theta\>Tr\left[\hat{P}_{a,t}(\theta)\hat{\rho}\right]}{\int
d\theta\>Tr\left[\hat{P}_{t}(\theta)\hat{\rho}\right]}$ (31)

where $\theta$ is the external time, $\hat{\rho}=\ket{\Psi}\bra{\Psi}$ is the
global state of the Universe,
$\hat{P}_{t}(\theta)=\hat{U}^{\dagger}(\theta)\hat{P}_{t}\hat{U}(\theta)$
(with $\hat{U}(\theta)=e^{-i\hat{H}\theta}$) is the projector relative to the
result $t$ for a clock measurement at external time $\theta$ and
$\hat{P}_{a,t}(\theta)=\hat{U}^{\dagger}(\theta)\hat{P}_{a,t}\hat{U}(\theta)$
is the projector relative to the result $a$ for a measurement on $S$ and $t$
for a measurement on $C$ at external time $\theta$ (we are working here in the
Heisenberg picture with respect to the external time $\theta$). The
generalization of equation (31) to multiple time measurements is given by gppt
:

$\begin{split}&p(a_{f}\>on\>S\>|\>t_{f}\>on\>C,a_{i},t_{i})=\\\ &=\frac{\int
d\theta\int
d\theta^{\prime}\>Tr\left[\hat{P}_{a_{f},t_{f}}(\theta)\hat{P}_{a_{i},t_{i}}(\theta^{\prime})\hat{\rho}\hat{P}_{a_{i},t_{i}}(\theta^{\prime})\right]}{\int
d\theta\int
d\theta^{\prime}\>Tr\left[\hat{P}_{t_{f}}(\theta)\hat{P}_{a_{i},t_{i}}(\theta^{\prime})\hat{\rho}\hat{P}_{a_{i},t_{i}}(\theta^{\prime})\right]}\end{split}$
(32)

which provides the conditional probability of obtaining $a_{f}$ on the system
$S$ at clock time $t_{f}$, given that a “previous”joint measurement of $S$ and
$C$ returns $a_{i}$, $t_{i}$.

The GPPT proposal needs a global state $\hat{\rho}$ commuting with the global
Hamiltonian and then it can be applied to our framework in the limit
$\delta\rightarrow 0$. We finally notice and emphasize that, in our framework,
the fact of having non-orthogonal time states does not constitute a problem in
the application of the GPPT theory. Indeed when we calculate the relative
state of the system $S$ to a certain clock value $t$ through (13), we find no
contributions from times $\neq t$ and interference phenomena are not present
even if the time states are not orthogonal nostro .

###  IV.2 Non-observable Universe as a clock

A condition to merge Canonical Typicality and PaW time is to have an
environment that is negligibly interacting with $S$. This is certainly not the
case in our everyday life where decoherence due to the surrounding environment
is often non-negligeable. A possible choice of a good clock non-interacting
with $S$ is the non-observable Universe (namely, a part of the global system
that is outside the light cone of $S$). In this respect, it is intriguing to
notice that the recent observations on the cosmic microwave background planck
together with the inflationary paradigm indicate that at the beginning of
cosmic inflation the Universe was in a pure state with highly-correlated
quantum fluctuations bianchi . Furthermore, it has been suggested that the
assumption that the observable and the non-observable Universe might be
entangled provides an argument in support of inflation vedraltemperature ,
vedralinflation . It is therefore somehow natural to speculate that the non-
observable Universe acts as a clock for the observable Universe. Indeed in
this framework the two requirements for a “good clock” are satisfied: the
clock and the system $S$ are non-interacting and, in addition, the dimension
of the clock is presumably larger than the dimension of $S$ (i.e. the the non-
observable Universe is bigger than the observable Universe). A very simple
estimate supports the consistency of the scenario. The spacing between the
energy levels in the clock space is:

$\delta
E^{(C)}_{i}=E^{(C)}_{i+1}+E^{(C)}_{i}=\frac{2\pi\hslash}{T}\left(r_{i+1}-r_{i}\right)$
(33)

where $r_{i+1}-r_{i}$ is an integer. Notice that $\delta
E^{(C)}_{min}=2\pi\hslash/T$ is the minimum energy step value, so that all
other energy values can be considered as integer multiples of this minimum
step. The $\delta E^{(C)}_{min}$ is inversely proportional to the value of
$T$, that is the time taken by the clock to return to its initial state
(notice that the framework we introduced in Section III leads to a cyclical
flow of time). So, considering the global system as the whole Universe, we can
relate $T$ to the current age of the Universe $T_{U}$ by assuming $T\geq
T_{U}\sim 13.8\times 10^{9}y\simeq 4.35\times 10^{17}s$ (and $t_{0}=0$ the
instant of the Big Bang). This leads to

$\delta
E^{(C)}_{min}=\frac{2\pi\hslash}{T}\leq\frac{2\pi\hslash}{T_{Universe}}\simeq
1.5\times 10^{-51}J.$ (34)

This upper limit for $\delta E^{(C)}_{min}$ is very small compared to other
energies on atomic scale and fits into our framework of constructing the
spectrum of the Hamiltonian $\hat{H}_{C}$ with integer multiples of a minimum
energy step.

##  V Conclusions

In this work we have merged Canonical Typicality and the PaW quantum time. We
consider a quantum Universe made by a small system $S$ and a large environment
which serves as a clock for $S$. Thanks to Canonical Typicality we know that
for almost all pure states in which the whole Universe can be, after tracing
over the environment the system $S$ is in a state of equilibrium described by
the canonical distribution. In the same scenario we find a Schrödinger-like
evolution corrected by a non-local term for the relative state of $S$ with
respect to the clock $\equiv$ environment. Canonical Typicality and dynamics
can coexist because in our protocol the action of tracing out the environment
is equivalent to tracing over all times: the trace over the environment
coincides with a temporal average.

## Acknowledgements

We acknowledge funding from the H2020 QuantERA ERA-NET Cofund in Quantum
Technologies projects QCLOCKS.

##  Appendix A Summary of Page and Wootters theory

We give here a brief review of the PaW theory. We consider the whole Universe
as being in a stationary state with zero eigenvalue (and therefore there is no
need for an external time), consistently with the the Wheeler-DeWitt equation

$\hat{H}\ket{\Psi}=0$ (35)

where $\hat{H}$ and $\ket{\Psi}$ are the Hamiltonian and the state of the
Universe respectively. We can then divide the Universe into two non-
interacting subsystems, the clock $C$ and the rest of the Universe $S$, and
thus the total Hamiltonian can be written as

$\hat{H}=\hat{H}_{C}+\hat{H}_{S}$ (36)

where $\hat{H}_{C}$ and $\hat{H}_{S}$ are the Hamiltonians acting on $C$ and
$S$ respectively. The condensed history of the system $S$ is written through
the entangled global stationary state
$\ket{\Psi}\in\mathcal{H}\subseteq\mathcal{H}_{C}\otimes\mathcal{H}_{S}$
(which satisfies the constraint (35)) as follows:

$\ket{\Psi}=\frac{1}{\sqrt{d_{C}}}\sum_{\tau}\ket{\tau}_{C}\otimes\ket{\phi_{\tau}}_{S}$
(37)

where $d_{C}$ is the dimension of the clock subspace and the states
$\left\\{\ket{\tau}_{C}\right\\}$ are eigenstates of the clock observable. We
notice here that a “good clock”must have $d_{C}\gg d_{S}$. Indeed, writing
$\ket{\Psi}$ in the energy eigenbasis, one finds that with $d_{C}\leq d_{S}$
it would not be possible to couple every energy state of $S$ to an energy
state of $C$ satisfying the constraint (35) and some states of $S$ would be
excluded from the dynamics nostro .

In this framework the relative state (in Everett sense everett ) of the
subsystem $S$ with respect to the clock $C$ can be obtained via conditioning

$\ket{\phi_{\tau}}_{S}=\frac{\braket{\tau}{\Psi}}{1/\sqrt{d_{C}}}.$ (38)

Note that, as mentioned before, equation (38) is the Everett relative state
definition of the subsystem $S$ with respect to the clock system $C$. As
pointed out in vedral , this kind of projection has nothing to do with a
measurement. Rather, $\ket{\phi_{\tau}}_{S}$ is a state of $S$ conditioned to
the clock $C$ in the state $\ket{\tau}_{C}$. Then, from equations (35), (36)
and (38), it is possible to derive the Schrödinger evolution for the relative
state of the subsystem $S$ with respect to the clock $C$ ($\hslash=1$):

$\ket{\phi_{\tau}}_{S}=e^{-i\hat{H}_{S}(\tau-\tau_{0})}\ket{\phi_{\tau_{0}}}_{S}$
(39)

where $\ket{\phi_{\tau_{0}}}_{S}=\sqrt{d_{C}}\braket{\tau_{0}}{\Psi}$, being
$\ket{\tau_{0}}_{C}$ the clock eigenstate taken as initial time.

We want to emphasize here that, as pointed out in lloydmaccone , in the PaW
approach the zero eigenvalue of $\hat{H}$ does not play a special role in
identifying the global state $\ket{\Psi}$. Indeed up to an irrelevant global
phase in the dynamics of $\ket{\phi_{\tau}}_{S}$, the global state
$\ket{\Psi}$ can be obtained also by imposing the constraint
$\hat{H}\ket{\Psi}=E\ket{\Psi}$ with real $E$. This consideration is useful
for the arguments we address in the main text.

##  Appendix B Proof of equation (18)

Here we show that the state

$\ket{\phi(t)}_{S}=\sum_{j}\sum_{i\in
I_{j}}c_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}$ (40)

is a solution for the equation

$i\frac{\partial}{\partial
t}\ket{\phi(t)}_{S}=\hat{H}_{S}\ket{\phi(t)}_{S}-\bra{t}\hat{\Delta}\ket{\Psi}$
(41)

(where $\hat{\Delta}=\sum_{j}\sum_{i\in
I_{j}}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}\bra{E^{(C)}_{i}}\bra{E^{(S)}_{j}}$)
and for equation (15). In the first case we substitute (40) in (41) thus
obtaining

$i\left(-i\sum_{j}\sum_{i\in
I_{j}}c_{ij}E^{(S)}_{j}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}+i\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}\right)=\\\
=\hat{H}_{S}\ket{\phi(t)}_{S}-\bra{t}\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}$ (42)
$\Rightarrow\hat{H}_{S}\ket{\phi(t)}_{S}-\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}=\hat{H}_{S}\ket{\phi(t)}_{S}-\bra{t}\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}.$ (43)

So we have:

$\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}=\bra{t}\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}\ket{E^{(C)}_{i}}\ket{E^{(S)}_{j}}$ (44)

that is an identity, considering that in the right-hand side of the equation
(44) $\braket{t}{E^{(C)}_{i}}=e^{iE^{(C)}_{i}t}$ and the global constraint on
energy gives $E^{(C)}_{i}=-E^{(S)}_{j}+\Delta_{ij}$ (remember that the total
energy term $E$ has been removed).

To verify that the (40) is a solution for (15) we have first to see how
$\hat{\Delta}(t,t^{\prime})$ acts on $\ket{\phi(t^{\prime})}_{S}$. We have

$\hat{\Delta}(t,t^{\prime})\ket{\phi(t^{\prime})}_{S}=\sum_{j}\sum_{k\in
I_{j}}\sum_{i\in
I_{j}}e^{iE^{(C)}_{k}t}\Delta_{kj}e^{-iE^{(C)}_{k}t^{\prime}}c_{ij}e^{iE^{(C)}_{i}t^{\prime}}\ket{E^{(S)}_{j}}.$
(45)

So, by substituting the state (40) in equation (15) and using (45), we obtain:

$\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt^{\prime}\sum_{j}\sum_{k\in
I_{j}}\sum_{i\in
I_{j}}e^{iE^{(C)}_{k}t}\Delta_{kj}e^{-iE^{(C)}_{k}t^{\prime}}c_{ij}e^{iE^{(C)}_{i}t^{\prime}}\ket{E^{(S)}_{j}}$
(46) $\Rightarrow\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}=\sum_{j}\sum_{i\in
I_{j}}c_{ij}\Delta_{ij}e^{-i(E^{(S)}_{j}-\Delta_{ij})t}\ket{E^{(S)}_{j}}$ (47)

where in the last step we used again the constraint on the energy
$E^{(C)}_{i}=-E^{(S)}_{j}+\Delta_{ij}$ and (see Appendix C for the proof):

$\int_{t_{0}}^{t_{0}+T}dt^{\prime}e^{i(E^{(C)}_{i}-E^{(C)}_{k})t^{\prime}}=T\delta_{i,k}.$
(48)

##  Appendix C Proof of equation (48)

We start here considering a generic state $\ket{\psi}\in\mathcal{H}_{C}$, that
is

$\ket{\psi}=\sum_{k=0}^{p}c_{k}\ket{E^{(C)}_{k}}$ (49)

where, we recall, $p+1=d_{C}$. We now apply in sequence the resolutions of the
identity $\mathbb{1}_{C}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\ket{t}\bra{t}$
and $\mathbb{1}_{C}=\sum_{n=0}^{p}\ket{E^{(C)}_{n}}\bra{E^{(C)}_{n}}$ thus
obtaining:

$\begin{split}\ket{\psi}&=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\ket{t}\braket{t}{\psi}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\sum_{k=0}^{p}c_{k}e^{itE^{(C)}_{k}}\ket{t}=\\\
\\\
&=\sum_{n=0}^{p}\ket{E^{(C)}_{n}}\bra{E^{(C)}_{n}}\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\sum_{k=0}^{p}c_{k}e^{itE^{(C)}_{k}}\ket{t}=\\\
\\\
&=\sum_{n=0}^{p}\sum_{k=0}^{p}c_{k}\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dte^{it(E^{(C)}_{k}-E^{(C)}_{n})}\ket{E^{(C)}_{n}}.\end{split}$
(50)

Since the state $\ket{\Psi}$ in equation (50) has to be equal to (49), we have
$\int_{t_{0}}^{t_{0}+T}dte^{it(E^{(C)}_{k}-E^{(C)}_{n})}=T\delta_{k,n}$.

##  Appendix D Proof of equation (21)

Here we prove equation (21), namely we show that, although the states
$\ket{t}$ are not orthogonal, we have

$\hat{\rho}_{S}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\braket{t}{\Psi}\braket{\Psi}{t}.$
(51)

Then, thanks to Canonical Typicality, we know that
$\hat{\rho}_{S}\approx\frac{1}{Z}e^{-\beta\hat{H}_{S}}$ where $\beta=dS(E)/dE$
is the inverse temperature (with $S(E)$ the entropy of $C$) and
$Z=Tr\left[e^{-\beta\hat{H}_{S}}\right]$ canonical . We start calculating the
partial trace of the global state $Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]$
through the energy basis in the subspace $C$:

$\begin{split}\hat{\rho}_{S}&=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\sum_{n}\braket{E^{(C)}_{n}}{\Psi}\braket{\Psi}{E^{(C)}_{n}}=\\\
&=\sum_{n}\sum_{j}\sum_{i\in I_{j}}\sum_{l}\sum_{k\in
I_{l}}c_{ij}c^{*}_{kl}\braket{E^{(C)}_{n}}{E^{(C)}_{i}}\braket{E^{(C)}_{k}}{E^{(C)}_{n}}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{l}}\end{split}$
(52)

where $I_{j}$ is the set of the environment levels such that
$E^{(C)}_{i}\in\left[E-E^{(S)}_{j},E-E^{(S)}_{j}+\delta\right]$. Now, being
$\delta\ll\Delta E^{(S)}_{j}~{}\forall j$, then the energy intervals for the
clock $\equiv$ environment coupled with each level $E^{(S)}_{j}$ of $S$ are
pairwise disjoint. So equation (52) becomes

$\hat{\rho}_{S}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\sum_{j}\sum_{i\in
I_{j}}\left|c_{ij}\right|^{2}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{j}}.$ (53)

Going instead to calculate the right-hand side of equation (51) we have

$\begin{split}\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\braket{t}{\Psi}\braket{\Psi}{t}&=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\sum_{j}\sum_{i\in
I_{j}}\sum_{l}\sum_{k\in
I_{l}}c_{ij}c^{*}_{kl}\braket{t}{E^{(C)}_{i}}\braket{E^{(C)}_{k}}{t}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{l}}=\\\
&=\sum_{j}\sum_{i\in I_{j}}\sum_{l}\sum_{k\in
I_{l}}c_{ij}c^{*}_{kl}\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dte^{-it(E^{(C)}_{i}-E^{(C)}_{k})}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{l}}\end{split}$
(54)

and, considering
$\int_{t_{0}}^{t_{0}+T}dte^{-it(E^{(C)}_{i}-E^{(C)}_{k})}=T\delta_{i,k}$ (see
Appendix C), we obtain

$\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\braket{t}{\Psi}\braket{\Psi}{t}=\sum_{j}\sum_{i\in
I_{j}}\left|c_{ij}\right|^{2}\ket{E^{(S)}_{j}}\bra{E^{(S)}_{j}}$ (55)

where we used again the fact that the relevant energy intervals for the clock
$\equiv$ environment coupled with each level $E^{(S)}_{j}$ of $S$ are pairwise
disjoint. Equation (55) that is the same of (53), so
$\frac{1}{T}\int_{t_{0}}^{t_{0}+T}dt\braket{t}{\Psi}\braket{\Psi}{t}=Tr_{C}\left[\ket{\Psi}\bra{\Psi}\right]=\hat{\rho}_{S}$.

##  Appendix E Proof of equation (27)

We calculate here the expectation value
$\langle\hat{X}\rangle_{t}=\bra{\phi(t)}\hat{X}\ket{\phi(t)}$, where
$\hat{X}=\sqrt{\frac{1}{2m\omega}}(\hat{a}+\hat{a}^{\dagger})$, considering
the relative state $\ket{\phi(t)}_{S}$ written as

$\ket{\phi(t)}_{S}=\sum_{i\in
I_{0}}c_{i0}e^{-i(E^{(S)}_{0}-\Delta_{i0})t}\ket{0^{(S)}}+\sum_{i\in
I_{1}}c_{i1}e^{-i(E^{(S)}_{1}-\Delta_{i1})t}\ket{1^{(S)}}.$ (56)

We have:

$\langle\hat{X}\rangle_{t}=\left(\sum_{i\in
I_{0}}c^{*}_{i0}e^{i(E^{(S)}_{0}-\Delta_{i0})t}\bra{0^{(S)}}+\sum_{i\in
I_{1}}c^{*}_{i1}e^{i(E^{(S)}_{1}-\Delta_{i1})t}\bra{1^{(S)}}\right)\times\\\
\times\sqrt{\frac{1}{2m\omega}}(\hat{a}+\hat{a}^{\dagger})\left(\sum_{i\in
I_{0}}c_{i0}e^{-i(E^{(S)}_{0}-\Delta_{i0})t}\ket{0^{(S)}}+\sum_{i\in
I_{1}}c_{i1}e^{-i(E^{(S)}_{1}-\Delta_{i1})t}\ket{1^{(S)}}\right)$ (57)
$\Rightarrow\langle\hat{X}\rangle_{t}=\sqrt{\frac{1}{2m\omega}}\left(\sum_{i\in
I_{0}}c^{*}_{i0}e^{i(E^{(S)}_{0}-\Delta_{i0})t}\bra{0^{(S)}}+\sum_{i\in
I_{1}}c^{*}_{i1}e^{i(E^{(S)}_{1}-\Delta_{i1})t}\bra{1^{(S)}}\right)\times\\\
\times\left(\sum_{i\in
I_{1}}c_{i1}e^{-i(E^{(S)}_{1}-\Delta_{i1})t}\ket{0^{(S)}}+\sum_{i\in
I_{0}}c_{i0}e^{-i(E^{(S)}_{0}-\Delta_{i0})t}\ket{1^{(S)}}+\sqrt{2}\sum_{i\in
I_{1}}c_{i1}e^{-i(E^{(S)}_{1}-\Delta_{i1})t}\ket{2^{(S)}}\right).$ (58)

From equation (58) we have:

$\langle\hat{X}\rangle_{t}=\sqrt{\frac{1}{2m\omega}}\sum_{i\in
I_{0}}\sum_{k\in
I_{1}}c^{*}_{i0}c_{k1}e^{i(E^{(S)}_{0}-E^{(S)}_{1}-\Delta_{i0}+\Delta_{k1})t}+\sqrt{\frac{1}{2m\omega}}\sum_{i\in
I_{0}}\sum_{k\in
I_{1}}c_{i0}c^{*}_{k1}e^{-i(E^{(S)}_{0}-E^{(S)}_{1}-\Delta_{i0}+\Delta_{k1})t}.$
(59)

Now, writing $c_{i0}=|c_{i0}|e^{i\varphi_{i0}}$,
$c_{k1}=|c_{k1}|e^{i\varphi_{k1}}$,
$\Delta\varphi^{(0,1)}_{ik}=\varphi_{k1}-\varphi_{i0}$ and considering that
$E^{(S)}_{1}-E^{(S)}_{0}=\omega$, we obtain:

$\langle\hat{X}\rangle_{t}=\sqrt{\frac{2}{m\omega}}\sum_{i\in I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|\cos((\omega+(\Delta_{i0}-\Delta_{k1}))t-\Delta\varphi^{(0,1)}_{ik})$
(60)

which is what we had to prove, since equation (60) is the same of the second
part of equation (27).

If we want instead to consider the expectation value
$\langle\hat{X}\rangle_{t}$ expressed in terms of the $\alpha_{j}(t)$, through
the definition (19) (that is $\alpha_{j}(t)=\sum_{i\in
I_{j}}c_{ij}e^{i\Delta_{ij}t}$), we can rewrite equation (59) as

$\langle\hat{X}\rangle_{t}=\sqrt{\frac{1}{2m\omega}}\left(\alpha^{*}_{0}(t)\alpha_{1}(t)e^{i(E^{(S)}_{0}-E^{(S)}_{1})t}+\alpha_{0}(t)\alpha^{*}_{1}(t)e^{-i(E^{(S)}_{0}-E^{(S)}_{1})t}\right).$
(61)

Writing $\alpha_{0}(t)=|\alpha_{0}(t)|e^{i\phi_{0}(t)}$,
$\alpha_{1}(t)=|\alpha_{1}(t)|e^{i\phi_{1}(t)}$,
$\Delta\phi(t)=\phi_{1}(t)-\phi_{0}(t)$ and considering again
$E^{(S)}_{1}-E^{(S)}_{0}=\omega$, (61) becomes:

$\langle\hat{X}\rangle_{t}=\sqrt{\frac{2}{m\omega}}|\alpha_{0}(t)||\alpha_{1}(t)|\cos(\omega
t-\Delta\phi(t)).$ (62)

Equation (62) is the same of the first part of equation (27) and shows the
expectation value $\langle\hat{X}\rangle_{t}$ expressed in terms of the time-
dependent coefficients $\alpha_{j}(t)$.

##  Appendix F Proof of equation (28)

We prove here that equation (27) reduces to (28) in the case of $t\ll
1/|\Delta_{i0}-\Delta_{k1}|$. We start considering the second part of equation
(27), that is:

$\begin{split}\langle\hat{X}\rangle_{t}=\sqrt{\frac{2}{m\omega}}\sum_{i\in
I_{0}}\sum_{k\in I_{1}}|c_{i0}||c_{k1}|\cos(\omega
t+(\Delta_{i0}-\Delta_{k1})t-\Delta\varphi^{(0,1)}_{ik}).\end{split}$ (63)

This can be rewritten as

$\begin{split}\langle\hat{X}\rangle_{t}&=\sqrt{\frac{2}{m\omega}}\left(\sum_{i\in
I_{0}}\sum_{k\in I_{1}}|c_{i0}||c_{k1}|\cos(\omega
t-\Delta\varphi^{(0,1)}_{ik})\cos((\Delta_{i0}-\Delta_{k1})t)\right)+\\\
&-\sqrt{\frac{2}{m\omega}}\left(\sum_{i\in I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|\sin(\omega
t-\Delta\varphi^{(0,1)}_{ik})\sin((\Delta_{i0}-\Delta_{k1})t)\right).\end{split}$
(64)

We impose now the condition $t\ll 1/|\Delta_{i0}-\Delta_{k1}|$, so we can
consider the Taylor expansions of $\cos((\Delta_{i0}-\Delta_{k1})t)$ and
$\sin((\Delta_{i0}-\Delta_{k1})t)$ thus obtaining

$\begin{split}\langle\hat{X}\rangle_{t}&=\sqrt{\frac{2}{m\omega}}\sum_{i\in
I_{0}}\sum_{k\in I_{1}}|c_{i0}||c_{k1}|\cos(\omega
t-\Delta\varphi^{(0,1)}_{ik})+\\\ &-\sqrt{\frac{2}{m\omega}}\sum_{i\in
I_{0}}\sum_{k\in I_{1}}|c_{i0}||c_{k1}|t(\Delta_{i0}-\Delta_{k1})\sin(\omega
t-\Delta\varphi^{(0,1)}_{ik})+\\\ &-\sqrt{\frac{1}{2m\omega}}\sum_{i\in
I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|t^{2}(\Delta_{i0}-\Delta_{k1})^{2}\cos(\omega
t-\Delta\varphi^{(0,1)}_{ik})+...\>\>.\end{split}$ (65)

With $\alpha_{0}(0)=|\alpha_{0}(0)|e^{i\phi_{0}(0)}$,
$\alpha_{1}(0)=|\alpha_{1}(0)|e^{i\phi_{1}(0)}$,
$\Delta\phi(0)=\phi_{1}(0)-\phi_{0}(0)$ and $\alpha_{j}(0)=\sum_{i\in
I_{j}}c_{ij}$, we have

$\sum_{i\in I_{0}}\sum_{k\in I_{1}}|c_{i0}||c_{k1}|\cos(\omega
t-\Delta\varphi^{(0,1)}_{ik})=|\alpha_{0}(0)||\alpha_{1}(0)|\cos(\omega
t-\Delta\phi(0)).$ (66)

So, considering up to the first order of approximation in
$t(\Delta_{i0}-\Delta_{k1})$, we have finally for the expectation value of the
position operator $\langle\hat{X}\rangle_{t}$:

$\langle\hat{X}\rangle_{t}\simeq\sqrt{\frac{2}{m\omega}}|\alpha_{0}(0)||\alpha_{1}(0)|\cos(\omega
t-\Delta\phi(0))+\\\ -\sqrt{\frac{2}{m\omega}}\sum_{i\in I_{0}}\sum_{k\in
I_{1}}|c_{i0}||c_{k1}|t(\Delta_{i0}-\Delta_{k1})\sin(\omega
t-\Delta\varphi^{(0,1)}_{ik})$ (67)

that is what we needed to show being (67) the same of equation (28).

## References

  * (1) S. Goldstein, J. L. Lebowitz, R. Tumulka and N. Zanghi, Phys. Rev. Lett. 96, 050403 (2006).
  * (2) S. Popescu, A. J. Short and A. Winter, Nat. Phys. 2, 754–758 (2006).
  * (3) S. Deffner and W. H. Zurek, New J. Phys. 18, 063013 (2016).
  * (4) D. N. Page and W. K. Wootters, Phys. Rev. D 27, 2885 (1983).
  * (5) W. K. Wootters, Int. J. Theor. Phys. 23, 701–711 (1984).
  * (6) V. Giovannetti, S. Lloyd and L. Maccone, Phys. Rev. D 92, 045033 (2015).
  * (7) E. Moreva, G. Brida, M. Gramegna, V. Giovannetti, L. Maccone and M. Genovese, Phys. Rev. A 89, 052122 (2014).
  * (8) E. Moreva, M. Gramegna, G. Brida, L. Maccone and M. Genovese, Phys. Rev. D 96, 102005 (2017).
  * (9) C. Marletto and V. Vedral, Phys. Rev. D 95, 043510 (2017).
  * (10) V. Vedral, Time, (Inverse) Temperature and Cosmological Inflation as Entanglement. In: R. Renner, S. Stupar, Time in Physics, 27-42, Springer (2017).
  * (11) L. Maccone and K. Sacha, Phys. Rev. Lett. 124, 110402 (2020).
  * (12) V. Baumann, F. Del Santo, A. R. H. Smith, F. Giacomini, E. Castro-Ruiz and C. Brukner, Quantum 5, 524 (2021).
  * (13) A. R. H. Smith and M. Ahmadi, Quantum 3, 160 (2019).
  * (14) A. Boette and R. Rossignoli, Phys. Rev. A 98, 032108 (2018).
  * (15) A. Boette, R. Rossignoli, N. Gigena and M. Cerezo, Phys. Rev. A 93, 062127 (2016).
  * (16) J. Leon and L. Maccone, Found. Phys. 47, 1597–1608 (2017).
  * (17) P. A. Hohn, A. R. H. Smith and M. P. E. Lock, Phys. Rev. D 104, 066001 (2021).
  * (18) T. Favalli and A. Smerzi, Quantum 4, 354 (2020).
  * (19) A. Vanrietvelde, P. A. Hoehn, F. Giacomini and E. Castro-Ruiz, Quantum 4, 225 (2020).
  * (20) A. R. H. Smith and M. Ahmadi, Nat. Commun. 11, 5360 (2020).
  * (21) N. L. Diaz, J. M. Matera and R. Rossignoli, Phys. Rev. D 100, 125020 (2019).
  * (22) N. L. Diaz and R. Rossignoli, Phys. Rev. D 99, 045008 (2019).
  * (23) C. Foti, A. Coppo, G. Barni, A. Cuccoli and P. Verrucchi, Nat. Commun. 12, 1787 (2021).
  * (24) AC. de la Hamette, S. L. Ludescher and M. P. Mueller, Entanglement/Asymmetry correspondence for internal quantum reference frames, arXiv:2112.00046.
  * (25) E. Castro-Ruiz, F. Giacomini, A. Belenchia and Č. Brukner, Nat. Commun. 11, 2672 (2020).
  * (26) F. Giacomini, Quantum 5, 508 (2021).
  * (27) H. Everett, The Theory of the Universal Wave Function. In: The Many Worlds Interpretation of Quantum Mechanics, Princeton University Press, Department of Physics, 1-140 (1957).
  * (28) P. Facchi and G. Garnero, Int. J. Geom. Methods Mod. Phys. 14, No. 08, 1740001 (2017).
  * (29) D. T. Pegg, Phys. Rev. A 58, 4307 (1998).
  * (30) K.V. Kuchar, Int. J. Mod. Phys. D 20, No. supp01, 3-86 (2011).
  * (31) R. Gambini, R. A. Porto, J. Pullin and S. Torterolo, Phys. Rev. D 79, 041501(R) (2009).
  * (32) Planck Collaboration, P. Ade et al., Planck 2015 results. XX. Constraints on inflation, Astron. Astrophys. 594, No. A20 (2016).
  * (33) E. Bianchi, L. Hackl and N. Yokomizo, Int. J. Mod. Phys. D 24, No. 12 (2015).
  * (34) S. Stupar and V. Vedral, Was inflation necessary for the existence of time?, arXiv:1710.04260.