# Multi-Relevance: Coexisting but Distinct Notions of Scale in Large Systems

Adam G. Kline  Department of Physics, The University of Chicago, Chicago IL
60637  Stephanie E. Palmer  Department of Organismal Biology and Anatomy and
Department of Physics, The University of Chicago, Chicago IL 60637

###### Abstract

Renormalization group (RG) methods have begun to see usage in biology and
computer science to support the search for simplifying structure in complex
models in high dimensional spaces. In this report, we show that mixture models
can be thought of as having multiple coexisting, exactly independent RG flows,
each with its own notion of scale, and define this property as multi-
relevance. To understand how this property affects RG analysis, we construct a
model that has two distinct notions of scale, each corresponding to the state
of an unobserved latent variable. In the regime where this latent variable can
be inferred using a linear classifier, the vertex expansion approach in non-
perturbative RG can be applied successfully but will give different answers
depending the choice of expansion point in state space. Next, in the regime
where linear estimation of the latent state fails, we show that an application
of the vertex expansion predicts a decrease in the total number of relevant
couplings from four to three. One consequence of this is that PCA-like coarse-
graining can eliminate important, resolvable parameters, and further, it may
only identify meaningless collective degrees of freedom. Taken together, these
results suggest that using RG in these new domains may lead to
oversimplification if multi-relevance is not accounted for.

Recently, renormalization group (RG) methods have been used in fields at the
boundaries of traditional physics and beyond Berman and Klinger (2022);
Goldman  _et al._ (2023); Erbin  _et al._ (2022); Lahoche  _et al._ (2021);
Mehta and Schwab (2014); Koch-Janusz and Ringel (2018); Di Sante  _et al._
(2022); Balog  _et al._ (2022); Bény and Osborne (2014); Cotler and Rezchikov
(2022); Strandkvist  _et al._ (2020); Pessoa and Caticha (2018); Jentsch and
Lee (2022). In the domain of theoretical biology, and in particular
neuroscience, there have been attempts to apply RG methods to data Bradde and
Bialek (2017); Meshulam  _et al._ (2019) and models Brinkman (2023); Tiberi
_et al._ (2022) in order to find simplifying structure through the language of
many-body physics. Yet, many biological systems are more complex than field
theories, and when directly adopting RG techniques, there is danger of
oversimplifying, or of failing to define important collective variables
correctly. This raises the question of how one should choose a cutoff scheme
when applying RG in biological contexts.

To frame the problem, let us consider the firing patterns of $N$ neurons in a
biological neural network. The joint distribution $P(\sigma)$ of spiking
events is a probability distribution over $\sigma\in\\{0,1\\}^{N}$, and $-\log
P(\sigma)$ is analogous to a Hamiltonian for classical Ising spins. Models for
this system have been fit to real data taken from vertebrate retina, and the
resulting energy landscapes have many local minima when constrained to a fixed
total firing rate Tkačik  _et al._ (2014); Prentice  _et al._ (2016). Because
each of these collective states is thought to encode a different set of
features in the stimulus, we expect the definitions of important collective
variables to change depending on the region of state space under
consideration. The problem is that most RG approaches implement linear coarse-
graining, meaning a single collective basis is applied globally. Using the
non-perturbative RG (NPRG), these statements can be made precise for
distributions over finitely many variables.

In Wetterich construction of NPRG, the cutoff is enforced by a adding a
regularization term to the Hamiltonian which limits fluctuations in IR
variables. Working with continuous degrees of freedom $\phi$, the Helmholtz
energy is given by

$W_{k}(J)=\log\int d\phi\,e^{-\mathcal{H}(\phi)-\frac{1}{2}\phi\cdot
R_{k}\cdot\phi+J\cdot\phi},$ (1)

where, $R_{k}$ adds mass to long wavelength collective variables ($|q|<k$),
suppressing their fluctuation contributions to the integral, but leaves short-
wavelength modes ($|q|>k$) unaffected. As $k\to 0$, it is required that
$R_{k}\to 0$, so that $W_{k}$ becomes the exact connected generating function.
In practice, the goal is to compute the effective average action

$\Gamma_{k}(\varphi)=\max_{J}\left\\{\varphi\cdot
J-W_{k}(J)\right\\}-\frac{1}{2}R_{k}\cdot\varphi^{2}$ (2)

which obeys the exact flow equation Wetterich (1993)

$\partial_{k}\Gamma_{k}(\varphi)=\frac{1}{2}\operatorname{Tr}\left\\{\partial_{k}R_{k}\cdot\left[R_{k}+\Gamma_{k}^{(2)}(\varphi)\right]^{-1}\right\\}\,.$
(3)

As $k\to 0$, $\Gamma_{k}$ approaches the effective potential, denoted
$\Gamma$, sometimes referred to as the Gibbs energy. NPRG methods are quite
diverse in implementation and motivation, and we refer the reader to a number
of excellent reviews for broader background Delamotte (2012); Kopietz  _et
al._ (2010); Bagnuls and Bervillier (2001); Dupuis  _et al._ (2021).

The fact that this RG procedure implements linear coarse-graining can be seen
in the quadratic dependence of the regulator term on $\phi$ in eq. (1) On the
other hand, a coarse-graining scheme which is defined to capture biologically
relevant features might use nonlinear coarse-graining Kline and Palmer (2022),
but this would invalidate the Wetterich equation (3). What is needed is some
way to add variability to the cutoff scheme without making the RG calculation
intractable.

Our first contribution is to provide a simple mechanism that allows for this
variability while leaving the essential flow equations unchanged. We show by
construction that some models are most naturally described by a collection of
coexisting, exactly independent RG flows, where each describes the collective
physics according to a different notion of scale. We term systems with this
property multi-relevant.

A system described by a state $\phi$ and Hamiltonian $\mathcal{H}(\phi)$ is
multi-relevant when there exists a finite set of Hamiltonians
$\\{\mathcal{H}_{z}(\phi)\\}$ which are polynomial in $\phi$ and which have
finite partition functions such that

$\exp\left(-\mathcal{H}(\phi)\right)=\sum_{z}\exp\left(-\mathcal{H}_{z}(\phi)\right)\,.$
(4)

In other words, $\exp(-\mathcal{H})$ is described by a mixture of models which
are normalizable and well-described by power-series expansions. Note that $z$
can be interpreted as a latent state label which is not explicitly part of the
system state. An immediate consequence of this condition is that the total
partition function is simply a sum of the component partition functions. This
can be expressed in terms of the effective actions as
$\exp\mathcal{L}[\Gamma](J)=\sum_{z}\exp\mathcal{L}[\Gamma_{z}](J)$,

where $\mathcal{L}[F](J)=\max_{\varphi}\left\\{\varphi\cdot
J-F(\varphi)\right\\}$ (recall that $\mathcal{L}$ is involutive), and by
$\Gamma_{z}$ we mean the effective potential obtained from $\mathcal{H}_{z}$
by running the RG flow (3) down to $k\to 0$. Because of this, each component
model $\mathcal{H}_{z}$ can be renormalized independently subject to its own
cutoff scheme $R_{z,k}$, then the resulting effective potentials can be
combined to give the full effective potential. In this sense, a multi-relevant
system has multiple exactly independent RG flows.

A similar mixture construction has been used to perform accurate multi-scale
modeling of complex systems in which metastable states have a significant
impact on the overall thermodynamics Liu  _et al._ (2022, 2019). In order to
calculate the full effective potential using our method, i.e. by decomposing
the flow into independent parts and re-combining, one must keep track of the
constant parts of the effective actions during the flow.

We now explore the consequences of this definition of multi-relevance in a
specific model. First, we will see that a standard NPRG approach applied to a
multi-relevant system can either succeed or fail to correctly describe the
low-energy physics. Using this, we will argue that principal component
analysis (PCA) fails to identify meaningful collective degrees of freedom in
data taken from multi-relevant systems, and leads to oversimplification in a
certain precise sense.

Our model has two mixture components, $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$.
Each is a finite, nonlocal analogue of scalar $\phi^{4}$ theory, centered
about a point $\varphi_{z}$, with $\phi\in\mathbb{R}^{N}$ and $N$ large.
Explicitly,

$\displaystyle\mathcal{H}_{z}(\phi+\varphi_{z})$
$\displaystyle=\sum_{ab}\frac{1}{2}(K_{ab}^{z}+u_{2}^{z}\delta_{ab})\phi_{a}\phi_{b}+\frac{u_{4}^{z}}{4!}\phi_{a}^{2}\phi_{b}^{2}$
(5) $\displaystyle\mathcal{H}(\phi)$
$\displaystyle=-\log[\exp(-\mathcal{H}_{A}(\phi))+\exp(-\mathcal{H}_{B}(\phi))]$
(6)

We will refer to (6) as the “$\phi^{4}$ mixture model”. Here, the usual
kinetic term has been replaced by a generic positive semi-definite matrix
$K_{ab}^{z}$ plus a “mass” part $u_{2}^{z}\delta_{ab}$, which is defined by
requiring that the smallest eigenvalue of $K^{z}$ is unity. The eigenvector
matrix $V^{z}$ of $K^{z}$ is chosen randomly from the Haar measure on $O(N)$,
and $\varphi_{z}$ has similarly a random orientation.

Although no notion of locality is present, we can still use RG. Following
Bradde and Bialek (2017), when analyzing model $z$, we take the collective
degrees of freedom to be the eigenvectors of $K^{z}$. Assume that the $K^{z}$
eigenvalues are well-described by a density $\rho_{z}(\lambda)$. In the
collective basis, the scale-free part of the kinetic term looks like

$\sum_{a}\lambda_{z,a}\phi_{a}^{2}=\int
d\lambda\rho_{z}(\lambda)\lambda\phi(\lambda)^{2}\sim\int\frac{d^{D_{z}}q}{(2\pi)^{D_{z}/2}}q^{2}|\phi(q)|^{2}$

By assumption, $\rho_{z}$ depends on $\lambda$ as a power law, and we can
identify the exponent with a spatial dimensionality:

$\rho_{z}(\lambda)\sim\lambda^{\alpha_{z}-1}\Rightarrow\alpha_{z}=D_{z}/2$

To be clear, we have not introduced anything like a spatial manifold in which
these degrees of freedom are localized, and this relation to dimension is only
an analogy. The spectral density exponent $\alpha$ controls the scaling
properties of couplings in the same way that dimensionality does in field
theory. The UV scale $\Lambda_{z}$ is the largest eigenvalue of $K^{z}$ and
scales as $N^{1/\alpha_{z}}$.

For our analysis we use the vertex expansion method with expansion points
$\varphi_{z}$. Although the NPRG formalism allows for all entries of the
effective vertices to be tracked, we perform a low-dimensional
parameterization in line with the parameters of the original Hamiltonian. Our
cutoff scheme is the Litim regulator Litim (2001), without field-strength
renormalization. For component $z$, this this takes the form

$(V_{z}^{T}R_{z,k}V_{z})_{ab}=\max\left\\{k-\lambda_{z,a},0\right\\}\delta_{ab}$
(7)

Our ansatz for the flowing effective action contains two couplings
$u_{2,k}^{z}$ and $u_{4,k}^{z}$ for each component $z$:

$\displaystyle\Gamma^{(2)}_{z,k\,ab}$
$\displaystyle=K_{ab}^{z}+u_{2,k}^{z}\delta_{ab}$
$\displaystyle\Gamma^{(4)}_{z,k\,abcd}$
$\displaystyle=\frac{u_{4,k}^{z}}{3}(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})$

The full RG analysis is given in supplemental information. We summarize the
essential results here. For $\alpha_{z}\in(3/2,2)$, (correspondingly
$D_{z}\in(3,4)$) the mass and interaction couplings $u_{z,2,k}$ and
$u_{z,4,k}$ are relevant and describe the asymptotic properties. They obey
standard flow equations, up to differences in numerical factors due to
neglecting sub-extensive diagrams. There are two fixed points, one at the
origin and one at the WF fixed point. There is a $\mathbb{Z}_{2}$ symmetry-
broken phase with two degenerate minima of the effective average action and a
symmetric phase with a single minimum at $\varphi_{z}$.

Because the $\phi^{4}$ mixture model is multi-relevant, $\Gamma$ can be
obtained from the $k\to 0$ limits of $\Gamma_{A,k}$ and $\Gamma_{B,k}$. These
individual analyses of $A$ and $B$ are simple, including only two relevant
parameters each. When we combine them to compute the whole effective potential
at the end, there are still only four relevant parameters which constrain
measurements, though their contributions take on a fairly complicated form. By
comparison, it is clear that if we had not used the mixture construction and
instead computed the vertex expansion flow for the whole $\Gamma_{k}$, our
final answer would look a lot simpler, namely it would be some polynomial in
the fields. This simpler representation is wrong in some regimes, since it
cannot capture the essentially non-polynomial dependence of $\Gamma(\phi)$ on
$\phi$, due to the mixture construction. Yet, there are regimes where the
standard approach is actually justified.

Whether standard vertex expansion succeeds or fails is essentially determined
by the parameters $\varphi_{A}$ and $\varphi_{B}$. These act as displacements
in state space, centering one $\phi^{4}$ model at $\varphi_{A}$ and the other
at $\varphi_{B}$. Let $\varphi=\varphi_{A}-\varphi_{B}$ be the separation
vector. In Fig. 1, we show how the probability density of $\phi$ projected
along the inter-basin axis $\hat{\varphi}$ depends on the separation. There is
a transition as $\varphi$ crosses some critical threshold $\varphi_{c}$, which
we calculate in the SI. To be concrete, we define this transition as the
largest separation at which $P(\hat{\varphi}\cdot\phi)$ has a local maximum at
$\hat{\varphi}\cdot\phi=0$. This approximates the separation beyond which the
$A$ and $B$ component densities can be separated by a hyperplane.

Figure 1: Demonstration of the coalescence transition. Top: spectrum of the
energy Hessian at the saddle point separating basins $A$ and $B$. Model
parameters are $N=50$, $\alpha_{A}=3/2$, $\Lambda_{A}=N^{1/\alpha_{A}}$,
$u_{2}^{A}=-\Lambda_{A}/10$, $u_{4}^{A}=0.58\Lambda_{A}^{s_{4}}$. All $A$ and
$B$ parameters are equal except basis. The grey lines depict the true
eigenvalues (neglecting $O(\varphi^{4})$ effects) while the black line
represents
$\lambda_{\varphi}=\hat{\varphi}\cdot\mathcal{H}^{(2)}\cdot\hat{\varphi}$. The
separation $\varphi$ at which $\lambda_{\varphi}$ crosses 0 is the mean-field
estimate. Bottom: probability density of the $\phi\cdot\hat{\varphi}$ as a
function of separation. Picking a value on the $x$-axis selects a single
distribution with a specified separation. The red dashed line indicates the
estimate (21), while the blue dashed line gives (22).

Consider first the well-separated phase, that is $\varphi>\varphi_{c}$. The
total energy has a saddle point near $\varphi_{A}$, and largely the local
structure looks like $\mathcal{H}_{A}$. This is because the coefficient data
due to the presence of $\mathcal{H}_{B}$ are suppressed exponentially with
respect to a power of the separation $\varphi$. For states $\phi$ near
$\varphi_{A}$,

$\displaystyle\mathcal{H}(\phi)$
$\displaystyle=\mathcal{H}_{A}(\phi)-e^{-\Delta\mathcal{H}_{BA}(\phi)}+O\left(e^{-2\Delta\mathcal{H}_{BA}(\phi)}\right),$

where $\Delta\mathcal{H}_{AB}=\mathcal{H}_{B}-\mathcal{H}_{A}$. Then, using
the regulator $R_{A,k}$ in (7), one can compute $\Gamma_{k}$ in the
neighborhood of $\varphi_{A}$, using the above as the initial conditions of
the flow. For sufficiently large separations $\varphi>\varphi_{c}$, the
perturbations to $\mathcal{H}_{A}$ due to $\mathcal{H}_{B}$ do not change the
fact that the calculation of $\Gamma_{k}$ recovers a good approximation of
$\Gamma_{A,k}$. This follows from the fact that when the basins are well-
separated, a series expansion of $\Gamma(\phi)$ will approximate that of
$\Gamma_{z}(\phi)$ for $\phi$ in a suitable neighborhood of $\varphi_{z}$. We
discuss this further in the SI.

The vertex expansion approach can therefore, in the well-separated case, be
sufficient to analyze a multi-relevant energy function $\mathcal{H}$. However,
to get a full description of $\Gamma$, one needs to perform several RG flows,
one at each $\varphi_{z}$. In general, each requires a different regulator
scheme, and will give different asymptotic properties than the others. In
practice, this is equivalent to using the mixture construction. It is
conceptually interesting to note that this analysis demonstrates how multi-
relevant models can have differing low-energy physics near different points in
state space. If there were additionally some notion of dynamics that led to
effective ergodicity breaking, one could consider the system getting trapped
in basin $A$ or $B$ for long periods of time, and in each case the other basin
should have little effect.

Next, let us examine where the vertex expansion calculation of the full
effective potential fails. In the coalesced phase, the separation $\varphi$ is
smaller than $\varphi_{c}$. For simplicity, we consider the case $\varphi=0$.
Expanding (6) about the point $\varphi_{A}=\varphi_{B}=0$ yields at the
quadratic level,

$\mathcal{H}^{(2)}_{ab}=K_{ab}+\frac{1}{2}(u_{2}^{A}+u_{2}^{B})\delta_{ab}\,,$

where $K=(K^{A}+K^{B})/2$. Because $K^{A}$ and $K^{B}$ are diagonal in
different, randomly chosen bases, their sum $2K$ defines a new notion of scale
that does not completely agree with either system $A$ or $B$. As before, we
define scale using the eigenvalues of $K$, with the collective degrees of
freedom (analogous to Fourier modes) defined by the eigenvectors. In our
construction, $K$ acquires a mass gap $K_{0}$, and after subtracting this off,
the remaining eigenvalues will have some asymptotically scale-free behavior.
Therefore let $\\{\lambda_{a}+K_{0}\\}_{a=1}^{N}$ denote the eigenvalues of
$K$, and construct the density $\rho(\lambda)$ of these eigenvalues as we did
in the finite $\phi^{4}$ case. The units associated to these eigenvalues are
defined through the scaling properties of this density. For example, since
$\rho(\lambda)\sim\lambda^{\alpha-1}$, an extensive quantity has dimension
$\alpha$.

Expanding (6) to higher orders in $\phi$, all of the terms come in the form of
symmetrized outer products involving $K^{A}$, $K^{B}$, and $\delta$, with
$\delta$ the identity matrix. For example, the $\mathcal{H}^{(4)}$ has a term
like $(K^{A})_{(ab}(K^{B})_{cd)}$, which we denote $(K^{A}\circ
K^{B})_{abcd}$. In the eigenbasis of $K$, the matrices $K^{A}$ and $K^{B}$ are
in fact approximately diagonal, but the diagonal elements pick up a mass gap
just as $K$ itself does. By defining $\tilde{K}^{A}=K^{A}-K^{A}_{0}\delta$ as
the scale-free part of the $K^{A}$ diagonal elements (and repeating for
$K^{B}$), we can write the expansion in terms of approximately scale-free
objects. After computing scaling dimensions, the relevant and marginal
couplings are captured by the ansatz (see SI):

$\displaystyle\Gamma_{k\,ab}^{(2)}$
$\displaystyle=(v^{A}_{k}\tilde{K}^{A}+v^{B}_{k}\tilde{K}^{B}+v^{\delta}_{k}\delta)_{ab}$
(8) $\displaystyle\Gamma_{k\,abcd}^{(4)}$
$\displaystyle=v^{\delta\delta}_{k}(\delta\circ\delta)_{abcd}$ (9)
$\displaystyle\Gamma^{(6)}_{k\,a\dots f}$
$\displaystyle=v^{\delta\delta\delta}_{k}(\delta\circ\delta\circ\delta)_{a\dots
f}$ (10)

It is worth noting that the spectral density exponent $\alpha$ determined by
this new quadratic coupling $K$ is smaller than $\alpha_{A}$ and $\alpha_{B}$
when we take them to be equal. Because of this, the sixtic interaction becomes
relevant and must be included. Yet, for $\alpha_{A}=\alpha_{B}$ in the range
$(3/2,2)$, no other terms appearing in the initial conditions generated by the
expansion of $\mathcal{H}$ are relevant. The relevant couplings we include
provide three unique constraints on the original model parameters, and the
marginal coupling initial conditions provide no additional constraints.

To compute the flow equations, we employ the Litim regulator (without field-
strength renormalization) as before. Keeping only extensive terms yields

$\displaystyle\dot{g}^{\delta}_{l}$
$\displaystyle=g^{\delta}_{l}+\frac{1}{6}\frac{g^{\delta\delta}_{l}}{(1+g^{\delta}_{l})^{2}}$
$\displaystyle\dot{g}^{\delta\delta}_{l}$
$\displaystyle=(2-\alpha)g^{\delta\delta}_{l}-\frac{1}{3}\frac{(g^{\delta\delta}_{l})^{2}}{(1+g^{\delta}_{l})^{3}}+\frac{1}{10}\frac{g^{\delta\delta\delta}}{(1+g^{\delta})^{2}}$
$\displaystyle\dot{g}^{\delta\delta\delta}_{l}$
$\displaystyle=(3-2\alpha)g^{\delta\delta\delta}_{l}+\frac{5}{3}\frac{(g_{l}^{\delta\delta})^{3}}{(1+g_{l}^{\delta})^{4}}-\frac{g_{l}^{\delta\delta}\,g_{l}^{\delta\delta\delta}}{(1+g_{l}^{\delta})^{3}}$

where $l=-\log(k/\Lambda)$. These flow equations yield three finite fixed
points: The Gaussian point at the origin, the WF point, and one other point
that appears for $\alpha\leq 3/2$. In effect, a standard application of the
vertex expansion approach predicts that the IR properties of the mixture model
in the coalesced phase are the same as some kind of nonlocal $\phi^{6}$ model
near and below 3 dimensions (but still above $D=8/3$ where $\phi^{8}$ becomes
relevant).

The crucial result that we come to now is that the vertex expansion approach
at $\varphi_{A}=\varphi_{B}$ has led to oversimplification. The primary
indication of this is there are only three relevant parameters, while
according to the separate analysis there should be four. As we construct ever
larger systems by increasing $N$, the net effect of irrelevant couplings on
the moments (measurements) decreases, and our ability to constrain their
initial values using data disappears. By contrast, we know that all four
parameters, $u_{z,2}$ and $u_{z,4}$, may be inferred if the mixture
construction is used, regardless of the basin separation $\varphi$. A
secondary indication is that the separate treatment admits an obvious
truncation scheme, while the combined analysis is significantly more
complicated. In the latter, due to the non-polynomial structure of the
microscopic model, an infinite number of irrelevant couplings with large
initial values are thrown away. However, these have a significant effect on
the flow of relevant and marginal couplings in the UV and for accurate
predictions must be included.

An important consequence of this calculation is that when PCA is applied to
data drawn from a multi-relevant system, it can lead to oversimplification.
This follows from the fact that by performing a single RG flow starting at
$\varphi_{A}=\varphi_{B}$ (or at the saddle point between the basins), we are
choosing a collective basis and cutoff scheme which emulates PCA. While it
still might succeed in identifying a low-dimensional description that spans
the minima of the effective potential, PCA identifies a collective basis which
is otherwise meaningless, and cannot be used to find $K^{A}$ and $K^{B}$.

In this letter, we have introduced the concept of multi-relevance. This
property endows a system with multiple coexisting and exactly independent RG
flows, each with its own potentially distinct notion of scale. The existence
of latent categorical variables is a simple mechanism which leads directly to
multi-relevance, essentially by definition, and so we expect these results to
be an important step towards implementing RG techniques in theoretical biology
and complex systems more broadly. We further found that in analyzing multi-
relevant systems, the vertex expansion approach in NPRG can be sufficient, but
only when the state space can be decomposed into regions with approximate
support for only one mixture component in each. In this case, the RG flow
needs to be computed in each such region. This condition is related to the
linear inferability of the latent variables $z$. In our $\phi^{4}$ mixture
model, the breakdown of this condition leads to a failure of the vertex
expansion approach to capture the dependence of average quantities on all of
the model parameters. In essence, this oversimplification is due to the issues
of approximating the mixture model Hamiltonian, which is non-polynomial in the
state variables, with a Taylor expansion.

##  I Acknowledgements

We would like to thank Cheyne Weis and Patrick Jentsch for useful feedback and
discussions during the preparation of this manuscript.

## References

  * Berman and Klinger (2022) D. S. Berman and M. S. Klinger, “The Inverse of Exact Renormalization Group Flows as Statistical Inference,”  (2022), arXiv:2212.11379 [hep-th].
  * Goldman  _et al._ (2023) S. Goldman, N. Lashkari, R. G. Leigh, and M. Moosa, “Exact Renormalization of Wave Functionals yields continuous MERA,”  (2023), arXiv:2301.09669 [hep-th, physics:quant-ph].
  * Erbin  _et al._ (2022) H. Erbin, V. Lahoche, and D. O. Samary, Machine Learning: Science and Technology 3, 015027 (2022), publisher: IOP Publishing.
  * Lahoche  _et al._ (2021) V. Lahoche, D. O. Samary, and M. Tamaazousti, arXiv:2002.10574 [cond-mat, physics:hep-th] (2021), arXiv: 2002.10574.
  * Mehta and Schwab (2014) P. Mehta and D. J. Schwab, “An exact mapping between the Variational Renormalization Group and Deep Learning,”  (2014), arXiv:1410.3831 [cond-mat, stat].
  * Koch-Janusz and Ringel (2018) M. Koch-Janusz and Z. Ringel, Nature Physics 14, 578 (2018).
  * Di Sante  _et al._ (2022) D. Di Sante, M. Medvidović, A. Toschi, G. Sangiovanni, C. Franchini, A. M. Sengupta, and A. J. Millis, “Deep Learning the Functional Renormalization Group,”  (2022), arXiv:2202.13268 [cond-mat].
  * Balog  _et al._ (2022) I. Balog, A. Rançon, and B. Delamotte, “Renormalization group and generalized Central Limit Theorems: The critical probability distributions of the order parameter of the Ising model,”  (2022), arXiv:2206.03769 [cond-mat, physics:hep-th].
  * Bény and Osborne (2014) C. Bény and T. J. Osborne, “Renormalisation as an inference problem,”  (2014), arXiv:1310.3188 [cond-mat, physics:hep-th, physics:quant-ph].
  * Cotler and Rezchikov (2022) J. Cotler and S. Rezchikov, “Renormalization Group Flow as Optimal Transport,”  (2022), number: arXiv:2202.11737 arXiv:2202.11737 [hep-th, physics:math-ph, physics:quant-ph].
  * Strandkvist  _et al._ (2020) C. Strandkvist, P. Chvykov, and M. Tikhonov, “Beyond RG: from parameter flow to metric flow,”  (2020), arXiv:2011.12420 [cond-mat].
  * Pessoa and Caticha (2018) P. Pessoa and A. Caticha, Entropy 20, 25 (2018), number: 1 Publisher: Multidisciplinary Digital Publishing Institute.
  * Jentsch and Lee (2022) P. Jentsch and C. F. Lee, “Simultaneous Critical Order-Disorder and Critical Phase Separation in Compressible Polar Active Fluids: Dynamic and Functional Renormalization Group Analyses,”  (2022), arXiv:2210.03830 [cond-mat, physics:physics].
  * Bradde and Bialek (2017) S. Bradde and W. Bialek, Journal of Statistical Physics 167, 462 (2017).
  * Meshulam  _et al._ (2019) L. Meshulam, J. L. Gauthier, C. D. Brody, D. W. Tank, and W. Bialek, Physical Review Letters 123, 178103 (2019).
  * Brinkman (2023) B. A. W. Brinkman, “Non-perturbative renormalization group analysis of nonlinear spiking networks,”  (2023), arXiv:2301.09600 [cond-mat, q-bio].
  * Tiberi  _et al._ (2022) L. Tiberi, J. Stapmanns, T. Kühn, T. Luu, D. Dahmen, and M. Helias, Physical Review Letters 128, 168301 (2022).
  * Tkačik  _et al._ (2014) G. Tkačik, O. Marre, D. Amodei, E. Schneidman, W. Bialek, and M. J. Berry, PLoS Computational Biology 10, e1003408 (2014).
  * Prentice  _et al._ (2016) J. S. Prentice, O. Marre, M. L. Ioffe, A. R. Loback, G. Tkačik, and M. J. B. Ii, PLOS Computational Biology 12, e1005148 (2016), publisher: Public Library of Science.
  * Wetterich (1993) C. Wetterich, Physics Letters B 301, 90 (1993), arXiv: 1710.05815.
  * Delamotte (2012) B. Delamotte, arXiv:cond-mat/0702365 852, 49 (2012), arXiv: cond-mat/0702365.
  * Kopietz  _et al._ (2010) P. Kopietz, L. Bartosch, and F. Schütz,  _Introduction to the Functional Renormalization Group_, Lecture Notes in Physics, Vol. 798 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010).
  * Bagnuls and Bervillier (2001) C. Bagnuls and C. Bervillier, Physics Reports 348, 91 (2001), arXiv: hep-th/0002034.
  * Dupuis  _et al._ (2021) N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, Physics Reports 910, 1 (2021), arXiv: 2006.04853.
  * Kline and Palmer (2022) A. G. Kline and S. E. Palmer, New Journal of Physics 24, 033007 (2022).
  * Liu  _et al._ (2022) Z.-K. Liu, Y. Wang, and S.-L. Shang, Journal of Phase Equilibria and Diffusion 43, 598 (2022).
  * Liu  _et al._ (2019) Z.-K. Liu, B. Li, and H. Lin, Journal of Phase Equilibria and Diffusion 40, 508 (2019).
  * Litim (2001) D. F. Litim, Physical Review D 64, 105007 (2001).

##  Appendix A Appendixes

###  A.1 NPRG on the finite $\phi^{4}$ model

In this appendix, we work through the NPRG analysis of a finite $\phi^{4}$
model as defined in the main text, Equation 5:

$\mathcal{H}(\phi)=\frac{1}{2}\sum_{ab}(K_{ab}+u_{2}\delta_{ab})\phi_{a}\phi_{b}+\frac{u_{4}}{4!}\phi_{a}^{2}\phi_{b}^{2}$

In Wetterich’s formulation of the NPRG, the flow is defined in terms a
slightly modified version of the 1PI generating functional, or Gibbs energy,
known as the effective average action:

$\Gamma_{k}(\varphi)=\varphi\cdot J-\log\int
d\phi\exp\left[-\mathcal{H}(\phi)-\frac{1}{2}R_{k\,ab}\,\phi_{a}\phi_{b}+J\cdot\phi\right]-\frac{1}{2}R_{k\,ab}\,\varphi_{a}\varphi_{b}$
(11)

The matrix $R$ is known as the regulator. It ensures that the the largest
eigenvalue of the propagator is no larger than $1/k$ (or $1/k^{2}$,
traditionally, though we choose an alternative parameterization for
convenience). The full propagator at scale $k$ can be found using the
effective action and the regulator:

$\left[G^{-1}_{k}\right]_{ab}=\frac{\partial}{\partial\varphi_{a}}\frac{\partial}{\partial\varphi_{b}}\Gamma_{k}(\varphi)+R_{k\,ab}$

It is important to note that unlike a hard cutoff scheme, the whole
configurational integral over $\phi$ in (11) is carried out. The regulator
takes the role as the cutoff, and given that it satisfies the necessary
boundary conditions, the resulting RG flow is totally well-defined for a
finite number of degrees of freedom. This is not the case with a hard cutoff
scheme, wherein one needs to rescale the total number of degrees of freedom in
order to map the parameters of the coarse-grained Hamiltonian back in to the
original space of couplings. This does not mean that there is not a rescaling
step eventually, rather this new rescaling step happens after the flow has
already been defined.

We use the vertex expansion, meaning the flow equations for vertices
$\Gamma^{(n)}_{a_{1},\dots,a_{n}}$ are found by differentiating the Wetterich
flow equation directly and keeping track of coefficient data:

$\displaystyle\dot{\Gamma}_{k}(\varphi)$
$\displaystyle=\frac{1}{2}\operatorname{Tr}\left[\dot{R}_{k}G_{k}(\varphi)\right]$
$\displaystyle\dot{\Gamma}_{ab}^{(2)}$
$\displaystyle=-\frac{1}{2}\operatorname{Tr}\left[\Gamma^{(4)}_{ab}G_{k}\dot{R}_{k}G_{k}\right]+\operatorname{Tr}\left[\Gamma^{(3)}_{a}G_{k}\Gamma^{(3)}_{b}G_{k}\dot{R}_{k}G_{k}\right]$
$\displaystyle\dot{\Gamma}_{abc}^{(3)}$ $\displaystyle=\dots$
$\displaystyle\dot{\Gamma}_{abcd}^{(4)}$
$\displaystyle=-\frac{1}{2}\operatorname{Tr}\left[\Gamma^{(6)}_{abcd}G_{k}\dot{R}_{k}G_{k}\right]+3\operatorname{Tr}\left[\Gamma^{(4)}_{ab}G_{k}\Gamma^{(4)}_{cd}G_{k}\dot{R}_{k}G_{k}\right]_{\text{(symmetrized)}}$
$\displaystyle\qquad+\text{terms involving }\Gamma^{(3)}$

From here on, we implicitly symmetrize where necessary. To evaluate these flow
equations we need to truncate the hierarchy at some highest power $n$ and
parameterize the vertices with an ansatz. It is easiest to start with a good
parameterization and verify that higher-order vertices will not be relevant.
Working from the high-temperature phase, all odd vertices such as
$\Gamma^{(3)}$ are set to zero. The initial condition on $\Gamma_{k}$ is
$\Gamma_{k\to\infty}=\mathcal{H}$. We approximate this by setting
$\Gamma_{\Lambda}=\mathcal{H}$. Truncating at order $n=4$, the effective
action can be written in the same form as $\mathcal{H}$:

$\displaystyle\Gamma^{(2)}_{k\,ab}$
$\displaystyle=(\lambda_{a}+u_{2\,k})\delta_{ab}$
$\displaystyle\Gamma^{(4)}_{k\,abcd}$
$\displaystyle=\frac{u_{4\,k}}{3}\left(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\right)$
$\displaystyle u_{2\,\Lambda}$ $\displaystyle=u_{2}$ $\displaystyle
u_{4\,\Lambda}$ $\displaystyle=u_{4}$

From the truncation and $\phi\to-\phi$ symmetry we therefore have the
following closed set of finitely many flow equations:

$\displaystyle\dot{\Gamma}_{ab}^{(2)}$
$\displaystyle=-\frac{1}{2}\operatorname{Tr}\left[\Gamma^{(4)}_{a,b}G_{k}\dot{R}_{k}G_{k}\right]$
$\displaystyle\dot{\Gamma}_{abcd}^{(4)}$
$\displaystyle=3\operatorname{Tr}\left[\Gamma^{(4)}_{ab}G_{k}\Gamma^{(4)}_{cd}G_{k}\dot{R}_{k}G_{k}\right]$

First consider inserting the ansatz for $\Gamma^{(4)}$ into the RHS of the
$\Gamma^{(2)}$ flow equation:

$\displaystyle\partial_{k}\Gamma^{(2)}_{k\,ab}$
$\displaystyle=-\frac{1}{2}\frac{u_{4\,k}}{3}\sum_{cd}{(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})\left(G_{k}\dot{R}_{k}G_{k}\right)_{cd}}$
$\displaystyle=-\frac{u_{4\,k}}{6}\left(\operatorname{Tr}\left(G_{k}\dot{R}_{k}G_{k}\right)\delta_{ab}+2\left(G_{k}\dot{R}_{k}G_{k}\right)_{ab}\right)$

The first and second terms on the second line differ by a factor proportional
to the number of degrees of freedom in the system, as we will show shortly.
For this reason, we refer to the first term as extensive, while the second is
sub-extensive. To actually calculate these terms, we need to introduce a
regulator and get an expression for the propagator. In this context, the Litim
regulator without field-strength renormalization can be written in the
$\Gamma^{(2)}$ eigenbasis as

$R_{k\,ab}=\max\left\\{k-\lambda_{a},0\right\\}\delta_{ab}$

This gives us

$\displaystyle\dot{R}_{k\,ab}$ $\displaystyle=0\,;\quad a>n_{k}$
$\displaystyle=\delta_{ab}\,;\quad a<n_{k}$ $\displaystyle G_{k\,ab}$
$\displaystyle=\left(\Gamma^{(2)}_{k}\right)^{-1}_{ab}\,;\quad a>n_{k}$
$\displaystyle=\frac{1}{k+u_{2\,k}}\delta_{ab}\,;\quad a>n_{k}$

And so

$\displaystyle\dot{\Gamma}^{(2)}_{k\,ab}$
$\displaystyle=-\frac{1}{2}\frac{u_{4\,k}}{3}\left(\frac{n_{k}}{(k+u_{2\,k})^{2}}\delta_{ab}+2\frac{1}{(k+u_{2\,k})^{2}}\dot{R}_{k\,ab}\right)$

By $n_{k}$ we mean the number of eigenvalues $\lambda_{a}$ below the cutoff
$k$. We will mainly be interested in the intermediate part of the flow, for
which $n_{k}$ is a number much larger than unity, even though it may be much
less than $N$. This is the regime in which the scale-free parts of these flow
equations are accurate and look similar to traditional results in infinite
systems. Dropping the sub-extensive term and inserting our ansatz for
$\Gamma^{(2)}$ into the RHS of the above equation, we find:

$\dot{u}_{2\,k}=-\frac{1}{6}\frac{n_{k}u_{4\,k}}{(k+u_{2\,k})^{2}}$

The computation for $\partial_{k}u_{4\,k}$ proceeds similarly. Inserting the
effective action ansatz into the RHS of the truncated flow equations yields

$\displaystyle\dot{\Gamma}^{(4)}_{k\,abcd}$
$\displaystyle=3\operatorname{Tr}\left[\Gamma^{(4)}_{ab}G_{k}\Gamma^{(4)}_{cd}G_{k}\dot{R}_{k}G_{k}\right]$
$\displaystyle=3u_{4\,k}^{2}\sum_{a^{\prime}b^{\prime}c^{\prime}d^{\prime}}\frac{1}{3}(\delta_{ab}\delta_{a^{\prime}b^{\prime}}+\delta_{aa^{\prime}}\delta_{bb^{\prime}}+\delta_{ab^{\prime}}\delta_{ba^{\prime}})\times$
$\displaystyle\qquad\qquad\qquad\>\>\frac{1}{3}(\delta_{cd}\delta_{c^{\prime}d^{\prime}}+\delta_{cc^{\prime}}\delta_{dd^{\prime}}+\delta_{cb^{\prime}}\delta_{dc^{\prime}})(G_{k}\dot{R}_{k}G_{k})_{a^{\prime}c^{\prime}}\,G_{k\,b^{\prime}d^{\prime}}$
$\displaystyle=\frac{1}{3}u_{4\,k}^{2}\left(\delta_{ab}\delta_{cd}\operatorname{Tr}\left(G^{3}_{k}\dot{R}_{k}\right)+\text{sub-
extensive terms}\right)$
$\displaystyle\to\frac{1}{3}\frac{n_{k}u_{4\,k}^{2}}{(k+u_{2\,k})^{3}}\delta_{ab}\delta_{cd}$

Where the RHS is symmetrized across $abcd$. Substituting on the LHS, we have
the scale-free part

$\dot{u}_{4\,k}=\frac{1}{3}\frac{n_{k}}{N}\frac{u_{4\,k}^{2}}{(k+u_{2\,k})^{3}}$

In this case, we can assign dimensions to the couplings simply by inspecting
the flow equations. The key is to note that, for $[k]=1$ by definition,

$\alpha=\partial_{\log\lambda}\log\rho(\lambda)+1\quad\Rightarrow\quad[n_{k}]=\alpha\,.$

A heuristic way to see this is to give $\lambda_{n}$ a simple expression
consistent with $\rho(\lambda)\sim\lambda^{\alpha-1}$. This is satisfied by

$\lambda_{n}\sim\Lambda\left(\frac{n}{N}\right)^{1/\alpha}\,.$

Proceeding, we can straightforwardly identify:

$s_{2}=[u_{2}]=1\quad s_{4}=[u_{4}]=2-\alpha$

Finally, define rescaled dimensionless couplings $g_{2}=k^{-s_{2}}u_{2}$ and
$g_{4}=k^{-s_{4}}u_{4}$, and let $l=-\log(k/\Lambda)$. The dimensionless flow
equations are given by:

$\displaystyle\partial_{l}g_{2\,l}$
$\displaystyle=g_{2\,l}+\frac{1}{6}\frac{g_{4\,l}}{(1+g_{2\,l})^{2}}$
$\displaystyle\partial_{l}g_{4\,l}$
$\displaystyle=(2-\alpha)g_{4\,l}-\frac{1}{3}\frac{g_{4\,l}^{2}}{(1+g_{2\,l})^{3}}$

These have the standard properties describing the $\phi^{4}$ system, up to
some differences in coefficients, and with anomalous scaling neglected. Save
for the infinite fixed points, when $\alpha>2$ the only physical ($g_{4}>0$)
fixed point is at $(0,0)$. For $\alpha<2$, the Wilson-Fisher fixed point
appears at

$g_{2}^{*}=\frac{2-\alpha}{4-\alpha}\qquad
g_{4}^{*}=\frac{24(2-\alpha)}{(4-\alpha)^{2}}\,.$

###  A.2 Approximate behavior of $\Gamma$ in well-separated phase

In the main text, we discussed that for a multi-relevant Hamiltonian
$\mathcal{H}$ with component Hamiltonians $\mathcal{H}_{z}$, the full
effective potential $\Gamma$ satisfies the exact relation

$\exp\mathcal{L}[\Gamma](J)=\exp W(J)=\sum_{z}\exp
W_{z}(J)=\sum_{z}\exp\mathcal{L}[\Gamma_{z}](J)$

In this section, we argue that in the well-separated phase, $\Gamma$ behaves
like $\Gamma_{z}$ for $\phi$ in the neighborhood of $\varphi_{z}$. Further, we
argue that although $\Gamma$ is convex in principle, it may still be useful to
consider it as the convex hull of a non-convex combination of the potentials
$\Gamma_{z}$.

To begin, we work out a toy model. Consider a state variable
$\phi\in\mathbb{R}$ which is distributed according to a mixture of Gaussians
with unit variance and with means at $\pm a$. That is,

$\mathcal{H}(x)=-\log\left[\exp(-\mathcal{H}_{+}(\phi))+\exp(-\mathcal{H}_{-}(\phi))\right]=-\log\left[\exp\left(-\frac{1}{2}(\phi-a)^{2}\right)+\exp\left(-\frac{1}{2}(\phi+a)^{2}\right)\right]$

The Helmholtz energy of $\mathcal{H}_{\pm}(\phi)$ is

$W_{\pm}(J)=\log\int
d\phi\exp(-\mathcal{H}_{\pm}(\phi)+J\phi)=\frac{1}{2}J^{2}\pm Ja$

By applying a Legendre transform to these functions we get the effective
potentials

$\Gamma_{z}(\phi)=\frac{1}{2}(\phi-\varphi_{z})^{2}=\mathcal{H}_{z}(\phi)$

Where $\varphi_{\pm}=\pm a$. Define $\varphi=\partial_{J}W(J)$, the average of
$\phi$ given source $J$. Then

$\varphi=J+a\tanh Ja$

So for $J>1/a$, $\varphi\to J+a$. As the separation $a$ is made larger, $1/a$
decreases, and this asymptotic behavior occurs for smaller values of $J$.
Consider the full effective potential in this case, where $a$ is large
(relative to unity, the standard deviation of a single mixture component).

When $J>1/a$, the sum over $z$ is dominated by the term $\exp W_{+}(J)$. This
yields:

$\Gamma(\varphi)=\mathcal{L}[W](\varphi)\approx\mathcal{L}[W_{+}](\varphi)=\Gamma_{+}(\varphi)$

In words, the full Helmholtz energy $W(J)$ is well-approximated by $W_{+}(J)$
for large positive sources $J$, i.e. $J>1/a$ (and $W_{-}(J)$ for $J<-1/a$).
Therefore, when we calculate $\Gamma$ from $W$ via Legendre transform, the
fact that $W$ is dominated individual components $W_{z}$ in various regions
causes $\Gamma$ to approximate $\Gamma_{z}$ in those regions.

How does this relate to our model and the definition of the critical
separation $\varphi_{c}$? First note that the toy model can be
straightforwardly extended to a larger state space and non-Gaussian component
energy functions. What is necessary is simply that the component densities are
well-separated, so that $\langle\phi\rangle_{J}$ can take on large values at
relatively small sources $J$. Next, our choice to define $\varphi_{c}$ in
terms of the maximum of the probability density of $\phi\cdot\hat{\varphi}$
was not crucial. The important point is that the separation between the basins
is greater than the variance (in the inter-basin direction) of each basin
individually. In fact, we compute $\varphi_{c}$ according to this condition
later in the SI. At $J=0$, the second derivative of $W(J)$ with respect to $J$
gives the total covariance $\langle\phi^{2}\rangle-\langle\phi\rangle^{2}$.
For large separations, this covariance value drops off quickly as $J$ is moved
away from zero, since a slight source “pushes” the system into one basin or
another. By assumption, these basins individually have much less variance than
the total distribution at $J=0$.

Typically, this structure is associated with symmetry breaking and ergodicity
breaking, wherein a small external source causes the system to find a state in
an ergodically broken region. In this paper, we have only been concerned with
distributions over finitely many variables, and adding dynamics would require
additional assumptions and structure. However, given the close analogies to
standard analysis in many-body systems, we believe that an effective
description of the structure of $\Gamma(\phi)$ in the well-separated phase is
provided by $\Gamma_{z}(\phi)$ when $\phi$ is near $\varphi_{z}$, i.e., the
minima of $\Gamma_{z}(\phi)$ can sometimes be treated as metastable states.
While this manifestly breaks convexity of $\Gamma$, the convex hull can always
be taken if needed, and otherwise $\Gamma$ can be interpreted as a constrained
free energy.

###  A.3 RG flow in the coalesced phase

####  A.3.1 Units and scaling dimensions

In the first appendix, the units of couplings $u_{2}$ and $u_{4}$ were not
determined until after the flow equations were found. In general, this is the
correct approach, namely that the units of different couplings should be
verified by showing that they remove explicit scale-dependence in the flow
equations. In the following discussion, the full (before irrelevant terms are
removed) flow equation hierarchy is very complicated and a similar, ‘by
inspection’ approach is difficult. In this section, we discuss heuristic
reasoning that provides scaling dimensions in terms of the original energy
function. However, we wish to emphasize that the true scaling dimensions may
not always appear when using this heuristic, since they are really defined by
a cutoff-scheme and are revealed in the flow equations.

The basic strategy follows Bradde and Bialek (2017) and Lahoche  _et al._
(2021), though we do not take into account non-power law dependence of $\rho$
on $\lambda$. As the cutoff scale $k$ is changed, the number $n_{k}$ of
remaining modes below this cutoff scales like $n_{k}\sim k^{\alpha}$, which we
denote by $[n_{k}]=\alpha$. What we would like to compute are the so-called
‘engineering dimensions’ of our couplings, which can be defined as the
critical exponents at the non-interacting (Gaussian) fixed point. First,
consider the trace of the propagator under a change of $k$, where $k$ is the
UV cutoff:

$\operatorname{Tr}_{k}[G]=\sum_{a<n_{k}}\langle\phi_{a}^{2}\rangle=\int_{0}^{k}d\lambda\,\rho(\lambda)\langle\phi(\lambda)^{2}\rangle=\int_{0}^{k}d\lambda\,\rho(\lambda)\lambda^{-1}=\frac{\alpha
N}{\Lambda^{\alpha}}\int_{0}^{k}d\lambda\,\lambda^{\alpha-2}\sim k^{\alpha-1}$

We conclude that

$[\phi]=\frac{1}{2}(\alpha-1)$

Note that under the identification $\lambda\sim q^{2}$ and $\alpha\sim D/2$ as
pointed out in the main text, this agrees with the standard result $[\phi]\sim
D/2-2$. Now consider a positive-definite matrix $C$ which is diagonal in the
same basis as $\langle\phi_{a}\phi_{b}\rangle$, and with eigenvalues
$\tilde{\lambda}_{a}$, where $a$ is ordered so that $\lambda_{a}$ are
monotonically increasing. Let the density of $\tilde{\lambda}_{a}$ be denoted
$\tilde{\rho}$, which integrates to $N$. Then define $\tilde{\alpha}(k)$ as

$\tilde{\alpha}=\partial_{\log
k}\log\tilde{n}_{k}\,;\quad\tilde{n}_{k}=\int_{0}^{k}d\lambda\,\tilde{\rho}(\lambda)$

This gives a mapping between the eigenvalues $\tilde{\lambda}_{a}$ of $C$ and
the eigenvalues $\lambda_{a}$ of $K$ through
$n_{\lambda}=\tilde{n}_{\tilde{\lambda}(\lambda)}$. When $\tilde{\rho}$ has
power-law dependence on $\tilde{\lambda}$, $\tilde{n}_{k}\sim
k^{\tilde{\alpha}}$, which yields:

$\tilde{\lambda}(\lambda)\sim\lambda^{\alpha/\tilde{\alpha}}\equiv\lambda^{\beta}$

Now the scaling dimension of $C$ can be computed in the same way:

$\operatorname{Tr}_{k}[CG]=\sum_{a<n_{k}}\langle
C_{aa}\phi_{a}^{2}\rangle=\sum_{a}\lambda_{a}^{-1}\tilde{\lambda}_{a}\sim\int_{0}^{k}d\lambda\,\rho(\lambda)\lambda^{\beta-1}\sim
k^{\alpha+\beta-1}$

Hence

$[C]=\frac{\alpha}{\tilde{\alpha}}=\beta$

Finally, let us consider a quartic coupling comprised of the symmetrized outer
product of two matrices $C_{1}$ and $C_{2}$ which satisfy the same properties
as $C$ in the previous calculation and have dimensions $\beta_{1}$ and
$\beta_{2}$ respectively. Suppose only the $n_{k}$ modes below the cutoff are
involved. Further, recall that these dimensions are defined as the exponents
of the non-interacting theory, so the four-point function can be broken into
two-point functions.

$\sum_{a,b<n_{k}}(C_{1}\circ
C_{2})_{aabb}\langle\phi_{a}^{2}\phi_{b}^{2}\rangle\sim\left(\int^{k}_{0}d\lambda\,\rho(\lambda)\lambda^{1-\beta_{1}}\right)\left(\int_{0}^{k}d\lambda\,\rho(\lambda)\lambda^{1-\beta_{2}}\right)+\text{sub-
leading terms}\sim k^{2\alpha+\beta_{1}+\beta_{2}-2}$

And so the engineering dimensions of couplings like $C_{1}\circ C_{2}$ are
just the sum of the dimensions of $C_{1}$ and $C_{2}$, that is:

$[C_{1}\circ C_{2}]=\beta_{1}+\beta_{2}\,.$

The dimensions of couplings in the energy can be found straightforwardly,
using these rules. As an example, consider the energy of the finite $\phi^{4}$
model, analyzed in the first appendix:

$\mathcal{H}(\phi)=\frac{1}{2}\sum_{ab}(K_{ab}+u_{2}\delta_{ab})\phi_{a}\phi_{b}+\frac{u_{4}}{4!}\phi_{a}^{2}\phi_{b}^{2}$

The dimension of $\mathcal{H}$ must be $\alpha$:

$[\mathcal{H}]=[K\cdot\phi^{2}]=1+2[\phi]=\alpha$

Following the rules above, we further find:

$[u_{2}]=1,\quad[u_{4}]=2-\alpha$

We reiterate that this method is heuristic, and rescaling by these dimensions
will not always give scale-free flows. The primary reasons why this
calculation may fail are anomalous scaling, the presence of operators like
$\tilde{K}^{z}$ (which will be defined in the next section) which do not
commute with the propagator $G$. When this latter situation is the case,
$\tilde{K}^{z}$ is not diagonal in the collective variable basis, and so the
$k$-scaling dimensions of loop sums like
$\operatorname{Tr}[\dot{R}_{k}G_{k}\tilde{K}^{A}G_{k}\tilde{K}^{A}G_{k}]$
aren’t simply given by the sum of dimensions of operators under the trace. We
are fortunate in this work that although this non-commutativity is present,
the $\tilde{K}^{z}$ matrices appear to be approximately or effectively
diagonal in the collective basis, and the loop sums have scaling dimensions
which can be naïvely approximated.

A third mechanism which can break this heuristic is the presence of operators
which are not scale-free. One example could be if $K$ were not scale free,
causing the eigenvalue density $\rho(\lambda)$ to not have simple power-law
dependence on $\lambda$. This can cause the scaling dimensions to change along
the flow Lahoche  _et al._ (2021) and demonstrates why it is necessary to
define the true scaling dimensions in terms of the flow equations and not a
priori. This contingency also demonstrates the power of NPRG formalism, since
the flow equations are defined whether or not a scale-free description is
available.

####  A.3.2 RG calculation

In the coalesced phase, the separation $\varphi$ between $\varphi_{A}$ and
$\varphi_{B}$ is smaller than a critical value $\varphi_{c}$, which is
calculated in the next appendix. Here we set $\varphi=0$ which significantly
simplifies the analysis. Expanding (6) with $\varphi_{A}=\varphi_{B}$ about
the point $\phi=0$ yields at the quadratic level,

$\displaystyle\mathcal{H}^{(2)}$
$\displaystyle=\frac{1}{2}\left(K^{A}+K^{B}+(u^{A}_{2}+u^{B}_{2})\delta\right)$
$\displaystyle=K+u_{2}\delta\,.$

In the above expressions, we use $\delta$ to denote the identity matrix. Let
us first consider the quadratic part, which we use to define scale. The
density of eigenvalues of $K$ is given by:

$\rho(\lambda)=\sum_{a=1}^{N}\delta(\lambda-\lambda_{a})\,,$

Which for large $N$ can be treated as an effectively smooth distribution, with
proper care. Approximating these eigenvalues to be power-law distributed,
define the scaling dimension $\alpha$ as:

$\alpha(k)=\partial_{\log
k}\log\int_{0}^{k}d\lambda\,\rho(\lambda)\approx\alpha$

Next, We can calculate the scaling dimensions of $K^{A}$ and $K^{B}$ as
discussed in the previous section. When written in the $K$-eigenbasis, $K^{A}$
and $K^{B}$ appear as approximately diagonal. The precise reasons for this are
beyond the scope of this work, so we merely take it as an experimental fact.
We must be careful to note that these matrices are not truly diagonal in the
$K$-eigenbasis. Additionally, $K^{A}$ and $K^{B}$ display mass gaps in their
diagonal elements just as $K$ does. We define their scale-free parts by
subtracting off these mass gaps:

$\displaystyle\tilde{K}^{A}$ $\displaystyle=K^{A}-K^{A}_{0}\delta$
$\displaystyle\tilde{K}^{B}$ $\displaystyle=K^{B}-K^{B}_{0}\delta$

This allows us to compute the $k$-scaling dimension of $\tilde{K}^{A}$, which
we also approximate as constant with respect to $k$:

$[\tilde{K}^{A}]=\beta_{A}(k)=\partial_{\log
k}\log\int_{0}^{k}d\lambda\,\sum_{a}\delta\left(\lambda-\operatorname{diag}_{a}\\{\tilde{K}^{A}\\}\right)\approx\beta_{A}$

And similarly for $\tilde{K}^{B}$, where by
$\operatorname{diag}\\{\tilde{K}^{A}\\}$ we mean the diagonal element in the
$K$-eigenbasis corresponding to the $\lambda+K_{0}$ eigenvalue of $K$.

In the topmost plot of Fig. 2, we show that by restricting to
$\alpha_{A}=\alpha_{B}$ and sweeping, $\beta_{A}$ and $\beta_{A}$ are
essentially unity, and that $\alpha$ has a predictable functional dependence.
While small deviations from power-law behavior may be present, we approximate
$\beta_{A}=\beta_{B}=1$ for the rest of the analysis, unless otherwise
specified. This approximation breaks down when $\alpha_{A}$ and $\alpha_{B}$
differ, or when $\Lambda_{A}$ and $\Lambda_{B}$ differ.

Figure 2: Numerical examination of the exponent $\alpha$ describing the
scaling properties of eigenvalues of $K$. Areas without grey shading denote
$\alpha_{A}=\alpha_{B}\in(1.5,2)$, which corresponds to $D\in(3,4)$. Top:
Estimates of the exponent $\alpha$ were obtained by explicitly constructing
$K$ and measuring power laws in the scale-free part of their spectra. Error
bars represent the root mean square error over 1 000 trials at each value of
$\alpha_{A}$. Matrix size is $N=75$. Black line is a polynomial fit which we
use to find scaling dimensions. Middle: Estimates of $\beta_{A}$ and
$\beta_{B}$ show no significant deviation from unity. Bottom: scaling
dimensions of various couplings. The only couplings which are relevant in the
range $\alpha_{A}=\alpha_{B}\in(1.5,2)$ are $v^{\delta\delta}$ and
$v^{\delta\delta\delta}$

We now turn our attention to the quartic couplings, which will contain most of
the rest of the building blocks of our ansatz.

$\displaystyle\mathcal{H}^{(4)}$ $\displaystyle=-\frac{3}{4}K^{A}\circ
K^{A}-\frac{3}{4}K^{B}\circ K^{B}+\frac{3}{2}K^{A}\circ K^{B}$
$\displaystyle-\frac{3(u_{2}^{A}-u_{2}^{B})}{2}K^{A}\circ\delta-\frac{3(u_{2}^{B}-u_{2}^{A})}{2}K^{B}\circ\delta$
$\displaystyle+\left(\frac{u_{4}^{A}}{2}+\frac{u_{4}^{B}}{2}-\frac{3((u_{2}^{A})^{2}+(u_{2}^{B})^{2})}{4}\right)\delta\circ\delta\,.$

The use of the symbol $\circ$ here denotes a symmetrized outer product:

$\displaystyle(K^{A}\circ K^{A})_{abcd}$
$\displaystyle=K^{A}_{(ab}K^{A}_{cd)}$
$\displaystyle=\frac{1}{4!}\sum_{\text{perms}\,p}K^{A}_{p_{1}p_{2}}K^{A}_{p_{3}p_{4}}\,.$

Again, it is useful to separate out the couplings into scale-free parts. This
yields:

$\displaystyle\mathcal{H}^{(4)}$
$\displaystyle=-\frac{3}{4}\tilde{K}^{A}\circ\tilde{K}^{A}-\frac{3}{4}\tilde{K}^{B}\circ\tilde{K}^{B}+\frac{3}{2}\tilde{K}^{A}\circ\tilde{K}^{B}$
$\displaystyle+\frac{3}{2}(u_{2}^{B}+K^{B}_{0}-u_{2}^{A}-K^{A}_{0})\tilde{K}^{A}\circ\delta$
$\displaystyle+\frac{3}{2}(u_{2}^{A}+K^{A}_{0}-u_{2}^{B}-K^{B}_{0})\tilde{K}^{B}\circ\delta$
$\displaystyle+\left(\frac{u_{4}^{A}}{2}+\frac{u_{4}^{B}}{2}-\frac{3}{4}\left(u_{2}^{A}+K_{0}^{A}-u_{2}^{B}-K_{0}^{B}\right)^{2}\right)\delta\circ\delta\,.$

Finally, for the sixth-order part, we discard all but the
$\delta\circ\delta\circ\delta$ term, which we shall justify shortly.

$\displaystyle\mathcal{H}^{(6)}_{a\dots f}$
$\displaystyle=\frac{15}{4}(u_{4}^{A}-u_{4}^{B})\times$
$\displaystyle\qquad(u_{2}^{B}-K_{0}^{B}-u_{2}^{A}+K_{0}^{A})(\delta\circ\delta\circ\delta)_{a\dots
f}$

For the moment, we explicitly include all of these operators in our ansatz for
the effective action $\Gamma_{k}$. This gives at least ten couplings in total;
three quadratic operators, six quartic, and one sixtic:

$\displaystyle\Gamma_{k\,ab}^{(2)}$
$\displaystyle=(v^{A}_{k}\tilde{K}^{A}+v^{B}_{k}\tilde{K}^{B}+v^{\delta}_{k}\delta)_{ab}$
(12) $\displaystyle\Gamma_{k\,abcd}^{(4)}$
$\displaystyle=(v^{AA}_{k}\tilde{K}^{A}\circ\tilde{K}^{A}+v^{AB}_{k}\tilde{K}^{A}\circ\tilde{K}^{B}+$
$\displaystyle\qquad
v^{BB}_{k}\tilde{K}^{B}\circ\tilde{K}^{B}+v^{A\delta}_{k}\tilde{K}^{A}\circ\delta\,+$
$\displaystyle\qquad
v^{B\delta}_{k}\tilde{K}^{B}\circ\delta\,+v^{\delta\delta}_{k}\delta\circ\delta)_{abcd}$
(13) $\displaystyle\Gamma^{(6)}_{k\,a\dots f}$
$\displaystyle=(\dots+v^{\delta\delta\delta}_{k}\delta\circ\delta\circ\delta)_{a\dots
f}$ (14)

Many of these couplings are irrelevant for $\alpha_{A}=\alpha_{B}\in(1.5,2)$,
including all sixtic couplings except $v^{\delta\delta\delta}$. This follows
from the dimensional analysis rules we discussed, together with our definition
of $\alpha$.

$\displaystyle[\phi]$
$\displaystyle=\frac{1}{2}(\alpha-1),\quad[\mathcal{H}^{(n)}]=\alpha\left(1-\frac{n}{2}\right)+\frac{n}{2}$
$\displaystyle[v^{A}]$
$\displaystyle=1-\beta_{A}=0,\quad[v^{B}]=1-\beta_{B}=0,\quad[v^{\delta}]=1$
$\displaystyle[v^{z_{1}z_{2}}]$
$\displaystyle=2-\alpha-\beta_{z_{1}}-\beta_{z_{2}}\,;\quad
z_{i}\in\\{A,B,\delta\\},\quad\beta_{\delta}=0$
$\displaystyle[v^{\delta\delta\delta}]$ $\displaystyle=3-2\alpha$

A few of these scaling dimensions are represented as a function of
$\alpha_{A}=\alpha_{B}$ in Fig. 2. Terms like $\tilde{K}^{A}\circ\delta$ and
$\tilde{K}^{A}\circ\tilde{K}^{B}$ have negative dimensions for the range of
spectral exponents we are considering, so these are termed irrelevant. Such
terms can be important for obtaining accurate estimates of non-universal
properties in terms of the original couplings, but deep in the IR, they only
have the effect of shifting values of relevant couplings. As discussed, these
dimensions are derived heuristically, and must be verified.

The ansatz (12-14) generates an approximately closed set of flow equations
given the truncation, since the new couplings generated in the flow are driven
only by sub-extensive contributions. As an example, we verify the predicted
value of $s_{z_{1}z_{2}}=[v^{z_{1}z_{2}}_{k}]$ above. Looking only at the
diagram with two 4-point vertices and keeping only extensive terms (see first
appendix for explanation),

$\partial_{k}v_{k}^{z_{1}z_{2}}=\frac{1}{3}\sum_{\begin{subarray}{c}{x\in\\{A,B,\delta\\}}\\\
y\in\\{A,B,\delta\\}\end{subarray}}v^{z_{1}x}_{k}v^{z_{2}y}_{k}\operatorname{Tr}\left[\dot{R}_{k}G_{k}C_{x}G_{k}C_{y}G_{k}\right]\,;\quad
C_{A}=\tilde{K}^{A},\,C_{B}=\tilde{K}^{B},\,C_{\delta}=\delta$

The LHS must have dimension $s_{z_{1}z_{2}}-1$, while the RHS is a little
tricker. Consider a single term:

$s_{z_{1}z_{2}}=1+s_{z_{1}x}+s_{z_{2}y}+[\Phi_{k}]\,;\quad\Phi_{k}^{xy}=\operatorname{Tr}\left[\dot{R}_{k}G_{k}C_{x}G_{k}C_{y}G_{k}\right]$
(15)

In reality, these $C$ matrices do not generally commute with $G_{k}$. Yet,
given the setup of our problem, the $\tilde{K}^{z}$ matrices effectively act
as if they are diagonal when in the $K$ eigenbasis. The catch is that their
diagonals are not eigenvalues and the scale free parts are no longer described
by the exponents $\alpha_{A}$ and $\alpha_{B}$, but instead with
$\alpha(\alpha_{A},\alpha_{B})$, as can be easily numerically verified.
Therefore,

$[\Phi^{xy}_{k}]=-3+\alpha+\beta_{x}+\beta_{y}$

Which shows, after direct substitution into (15) that the free-field guess at
$s_{z_{1}z_{2}}$ is correct.

The ansatz (12-14) can be used to find flow equations for all couplings, and
these can be integrated numerically. At present we are only interested in
asymptotic properties, so we can thin out the effective action ansatz by
removing irrelevant terms. This takes us from ten down to five couplings:

$\displaystyle\Gamma_{k\,ab}^{(2)}$
$\displaystyle=(v^{A}_{k}\tilde{K}^{A}+v^{B}_{k}\tilde{K}^{B}+v^{\delta}_{k}\delta)_{ab}$
(16) $\displaystyle\Gamma_{k\,abcd}^{(4)}$
$\displaystyle=v^{\delta\delta}_{k}(\delta\circ\delta)_{abcd}$ (17)
$\displaystyle\Gamma^{(6)}_{k\,a\dots f}$
$\displaystyle=v^{\delta\delta\delta}_{k}(\delta\circ\delta\circ\delta)_{a\dots
f}$ (18)

Because $v^{A\delta}$ and $v^{B\delta}$ are irrelevant and have been dropped,
$v^{A}$ and $v^{B}$ are not significantly renormalized by any couplings in
this ansatz. The dimensionless flow equations for relevant parameters,
computed using the Litim regulator and not accounting for anomalous scaling,
are given by:

$\displaystyle\dot{g}^{\delta}_{l}$
$\displaystyle=g^{\delta}_{l}+\frac{1}{6}\frac{g^{\delta\delta}_{l}}{(1+g^{\delta}_{l})^{2}}$
$\displaystyle\dot{g}^{\delta\delta}_{l}$
$\displaystyle=(2-\alpha)g^{\delta\delta}_{l}-\frac{1}{3}\frac{(g^{\delta\delta}_{l})^{2}}{(1+g^{\delta}_{l})^{3}}+\frac{1}{10}\frac{g^{\delta\delta\delta}}{(1+g^{\delta})^{2}}$
$\displaystyle\dot{g}^{\delta\delta\delta}_{l}$
$\displaystyle=(3-2\alpha)g^{\delta\delta\delta}_{l}+\frac{5}{3}\frac{(g_{l}^{\delta\delta})^{3}}{(1+g_{l}^{\delta})^{4}}-\frac{g_{l}^{\delta\delta}\,g_{l}^{\delta\delta\delta}}{(1+g_{l}^{\delta})^{3}}$

###  A.4 Calculation of multi-relevance breakdown separation $\varphi_{c}$

####  A.4.1 Mean-field saddle point method

To begin, assume the $A$ and $B$ models have all the same parameters, e.g.
$u_{2}^{A}=u_{2}^{B}$, etc., but the eigenbases of $K^{A}$ and $K^{B}$ are
random with respect to each other. (Each drawn from the Haar measure on
$O(N)$.) The saddle point state is given by:

$\phi_{\text{sp}}=\frac{1}{2}M^{-1}\,(M^{B}-M^{A})\cdot\varphi+O(\varphi^{3})$

Where $M^{A}=K^{A}+u_{2}^{A}$ and $M=M^{A}+M^{B}$. The truncation in powers of
$\varphi$ is due to the fact that $\varphi_{c}$ should go as $\Lambda^{-1/2}$,
which we take to be small. Set $\varphi_{A}=\varphi/2$ and
$\varphi_{B}=-\varphi_{A}$ without loss of generality. The saddle point
condition is approximately satisfied at the origin:

$\frac{1}{2}M^{-1}\,(M^{B}-M^{A})\cdot\varphi\approx\frac{1}{2}M^{-1}\left(\frac{1}{N}\operatorname{Tr}M^{B}-\frac{1}{N}\operatorname{Tr}M^{A}\right)\varphi\approx
0$

An important approximation we made here is that since $\varphi$ is randomly
oriented with respect to the eigenbases of $K^{A}$ and $K^{B}$, it is
approximately an eigenvector for the ranges $\alpha$ we are dealing with, as
can be verified numerically. Explicitly,

$M^{A}\cdot\varphi\approx\frac{1}{N}\left(\operatorname{Tr}M^{A}\right)\varphi\,.$

At the origin, the energy Hessian up to $O(\varphi^{4})$ is

$\displaystyle\mathcal{H}^{(2)}_{ab}$
$\displaystyle=M_{ab}-\frac{1}{4}(M\cdot\varphi)_{a}(M\cdot\varphi)_{b}+\frac{1}{2^{4}3}(u_{4}^{A}+u_{4}^{B})(\varphi^{2}\delta_{ab}+2\varphi_{a}\varphi_{b})+O(\varphi^{4})$
(19)
$\displaystyle=M_{ab}+\frac{1}{2^{3}3}u_{4}\varphi^{2}\delta_{ab}-\frac{1}{4}\left(\left(\frac{1}{N}\operatorname{Tr}M\right)^{2}-\frac{1}{3}u_{4}\right)\varphi_{a}\varphi_{b}$
(20)

Though $\varphi$ is not an exact eigenvector of this Hessian, we approximate
it as one. Its ‘eigenvalue’ can be found from the above:

$(\mathcal{H}^{(2)}\cdot\varphi)_{a}=\left(\frac{1}{N}\operatorname{Tr}{M}+\frac{1}{2^{3}}u_{4}\varphi^{2}-\frac{1}{2^{2}}\left(\frac{1}{N}\operatorname{Tr}M\right)^{2}\varphi^{2}\right)\varphi_{a}$

Setting this to zero, we can solve for the point at which curvature along the
separation vector between the basins goes from positive to negative:

$\displaystyle\varphi_{c}^{2}=4\left(\left(\frac{1}{N}\operatorname{Tr}M\right)^{2}-\frac{1}{2}u_{4}\right)^{-1}$
(21)

####  A.4.2 Mixture component variance calculation

The alternative method for calculating the critical separation is to leverage
the fact that our model is a mixture of densities:

$P_{z}(\phi)\propto\exp(-\mathcal{H}_{z}(\phi))\,,\quad
P(\phi)\propto\exp(-\mathcal{H}_{A}(\phi))+\exp(-\mathcal{H}_{B}(\varphi))$

Because we have assumed all of the same parameters for the $A$ and $B$ states
(except basis),

$P(\phi)=\frac{1}{2}P_{A}(\phi)+\frac{1}{2}P_{B}(\phi)$

To calculate the probability density of $\phi$ along $\hat{\varphi}$, the unit
vector along the basin separation, we evaluate

$\displaystyle P(\hat{\varphi}\cdot\phi=s)$ $\displaystyle=\int
d\phi\,\delta(s-\hat{\varphi}\cdot\phi)P(\phi)$
$\displaystyle=\frac{1}{2}P_{A}(\hat{\varphi}\cdot\phi=s)+\frac{1}{2}P_{B}(\hat{\varphi}\cdot\phi=s)$

In the independent phase, $P_{A}(\hat{\varphi}\cdot\phi)$ does not
significantly overlap $P_{B}(\hat{\varphi}\cdot\phi)$. If $A$ and $B$ are in
their disordered phases, these distributions will be Gaussian, while in their
low temp phases they will each look like a mixture two Gaussians placed at
$\pm\hat{\varphi}_{0}|\langle\phi_{0}\rangle_{z}|$, where by
$|\langle\phi_{0}\rangle_{z}|$ we simply mean the magnetization in the
standard sense. In both cases the coalescence transition occurs roughly when
the distance between the centers of these distributions is twice the standard
deviation of one of them, causing overlap:

$\left(\frac{\varphi_{c}}{2}\right)^{2}=\langle(\hat{\varphi}\cdot\phi)^{2}\rangle_{z}=\sum_{ab}\hat{\varphi}_{a}\hat{\varphi}_{b}\left(\Gamma_{z,k\to
0}^{(2)}\right)^{-1}_{ab}$ (22)

From the ansatz for $\Gamma_{k}$,

$\left(\Gamma_{z,0}^{(2)}\right)_{ab}=(\lambda^{z}_{a}+u_{2,0}^{z})\delta_{ab}$

In the $K^{z}$ eigenbasis. Using the fact that $\varphi$ has a random
orientation in this basis,

$\sum_{ab}\hat{\varphi}_{a}\hat{\varphi}_{b}\left(\Gamma_{z,0}^{(2)}\right)^{-1}_{ab}\approx\frac{1}{N}\operatorname{Tr}\left(\Gamma_{z,0}^{(2)}\right)^{-1}$

In the disordered phase, $u^{z}_{2,0}>0$, so

$\frac{1}{N}\operatorname{Tr}\left(\Gamma_{z,0}^{(2)}\right)^{-1}=\frac{1}{N}\int_{0}^{\Lambda_{z}}d\lambda\,\rho_{z}(\lambda)\frac{1}{\lambda+u_{2,0}^{z}}$

In the disordered phase, $u_{2,0}^{z}<0$, so expand the potential instead
around one of the minima at
$\langle\phi_{0}\rangle=\pm(-6u_{2,0}^{z}/u_{4,0}^{z})^{1/2}$. The total
variance is $\langle\phi_{0}\rangle^{2}$ plus trace the inverse of
$\Gamma^{(2)}$ evaluated at the symmetry-broken expansion point:

$\sum_{ab}\hat{\varphi}_{a}\hat{\varphi}_{b}\left(\Gamma_{z,0}^{(2)}\right)^{-1}_{ab}\approx\frac{1}{N}\left(-\frac{6u_{2,0}^{z}}{u_{4,0}^{z}}+\int_{0}^{\Lambda_{z}}d\lambda\,\rho_{z}(\lambda)\frac{1}{\lambda-2u_{2,0}^{z}}\right)$

There is a normalization factor on $\rho_{z}$ so that it integrates to $N$.
For large $N$, our assumption is that

$\rho_{z}(\lambda)\to\frac{N\alpha_{z}}{\Lambda_{z}^{\alpha_{z}}}\lambda^{\alpha_{z}-1}$

Hence the disordered phase gives

$\left(\frac{\varphi_{c}}{2}\right)^{2}\approx\frac{\alpha_{z}}{\Lambda_{z}^{\alpha_{z}}}\int_{0}^{\Lambda_{z}}d\lambda\,\frac{\lambda^{\alpha_{z}-1}}{\lambda+u_{2,0}^{z}}$

While the ordered phase gives

$\left(\frac{\varphi_{c}}{2}\right)^{2}\approx-\frac{6}{N}\frac{u_{2,0}^{z}}{u_{4,0}^{z}}+\frac{\alpha_{z}}{\Lambda_{z}^{\alpha_{z}}}\int_{0}^{\Lambda_{z}}d\lambda\,\frac{\lambda^{\alpha_{z}-1}}{\lambda+u_{2,0}^{z}}$


# Testing Quantum Coherence in Stochastic Electrodynamics with Squeezed
Schrödinger Cat States

Wayne Cheng-Wei Huang  Center for Fundamental Physics, Northwestern
University, Evanston, Illinois 60208, USA  Herman Batelaan  email:
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Nebraska-Lincoln, Lincoln, Nebraska 68588, USA

###### Abstract

The interference pattern in electron double-slit diffraction is a hallmark of
quantum mechanics. A long standing question for stochastic electrodynamics
(SED) is whether or not it is capable of reproducing such effects, as
interference is a manifestation of quantum coherence. In this study, we use
excited harmonic oscillators to directly test this quantum feature in SED. We
use two counter-propagating dichromatic laser pulses to promote a ground-state
harmonic oscillator to a squeezed Schrödinger cat state. Upon recombination of
the two well-separated wavepackets, an interference pattern emerges in the
quantum probability distribution but is absent in the SED probability
distribution. We thus give a counterexample that rejects SED as a valid
alternative to quantum mechanics.

##  I Introduction

Over the past decades, there has been sustained interest in developing
classical alternatives to quantum mechanics (QM) with the goal of solving the
quantum-classical boundary problem. Despite the proposed classical
alternatives Pena1996 ; Pena2015 ; Cavalleri , there is a lack of quantitative
tests of such theories against QM, mostly because analytic solutions to
concrete physical systems such as two-level atoms have not be found. Arguably
one of the most developed classical alternatives is stochastic electrodynamics
(SED) Milonni ; Boyer1975-1 . Studies of SED harmonic systems have found many
examples that are in exact agreement with QM. These include the retarded van
der Waals force Boyer1973 , ground state distribution of harmonic oscillators
Boyer1975-2 ; Huang2013 , Landau diamagnetism Boyer1980 ; Boyer2016 , Planck
spectrum of blackbody radiation Boyer1984 ; Boyer2017 , and Debye specific-
heat law for solids Santos . Numerical studies of hydrogen have given some
qualitative features Puthoff ; Cole , but have not led to a clear success
Nieuwenhuizen2015 ; Nieuwenhuizen2016 . Recently, it was further shown that
parametric interaction can give rise to discrete SED excitation spectra that
is in excellent agreement with QM predictions Huang2015 .

Figure 1: Two counter-propagating laser fields with frequencies $\omega_{1}$
and $\omega_{2}$, collaboratively drive the harmonic oscillator with a
spatially modulated Kapitza-Dirac force in the direction of wave propagation.
The force has a modulation periodicity of $2\pi c/(\omega_{1}+\omega_{2})$ and
it oscillates at the difference frequency $\omega_{1}-\omega_{2}$. A particle
subject to the perturbation of the classical zero-point electromagnetic field
can get pushed to either directions.

Figure 2: Time evolution of QM and SED probability distributions after laser
excitation. (a) The Wigner function of the QM oscillator is plotted after the
laser excitation. Positive values are color coded in red and negative values
are in blue. Quantum coherence between the two well-separated squeezed states
is manifested by the fringe structure in the center. The Wigner function
rotates counter-clockwise at the oscillator’s resonant frequency
$\omega_{0}=10^{16}$ (rad/s). At the moment when the distribution is depicted,
the position quadrature uncertainty $\sigma_{x}$ of the squeezed state is
smaller than that of the ground state $\Delta_{x}=\sqrt{\hbar/2m\omega_{0}}$.
Meanwhile, the quadrature uncertainty product satisfies the Heisenberg
relation $\sigma_{x}\sigma_{p}=\hbar/2$ at all times. (b) The two wavepackets
oscillate back and forth in the harmonic potential giving rise to a double
sinusoidal trajectory of the QM probability distribution. An interference
pattern appears when the two wavepackets merge. (c) The phase space
distribution of the SED oscillator shows two well-separated sub-ensembles.
Each sub-ensemble has a squeezed structure that mimics the QM squeezed state
shown in (a). (d) As the SED phase space distribution rotates counter-
clockwise, the probability distribution bundles into two ”macroscopic”
trajectories. In each trial of the SED simulation, there is no knowledge which
macroscopic trajectory a particle will follow unless the initial phase of the
background zero-point field is known. No interference-like patterns are found
in the SED probability distribution when the two macroscopic trajectories of
distributions cross.

However, a major drawback to the generality of the SED approach is that none
of the investigated effects involves quantum coherence. In light of this, some
have proposed to study electron double-slit diffraction within the framework
of SED as interference is a manifestation of quantum coherence Boyer1975-1 ;
Huang2013 ; Bach . Within the SED community the proposed view of electron
diffraction is that the double-slit poses boundary conditions that modify the
classical zero-point electromagnetic field and in turn it acts as a guiding
wave for free electrons Cavalleri ; Boyer1975-1 ; Avendano ; Kracklauer1992 .
The appeal of this idea is that the guiding field can be affected by both
slits, while the particle passes through only one slit, similar to the idea
that has pushed oil droplet analogues Wolchover ; Couder ; Andersen ; Bohr ;
Batelaan2016 ; Pucci . As appealing as the idea may sound, so far there has
been no concrete calculation or simulation demonstrating such an effect
because of two major theoretical obstacles: (1) the effective spectrum of the
zero-point field is unbounded for free electrons, (2) The radiation damping of
free electrons gives rise to runaway solutions.

Rather than focusing on the specific theoretical difficulties that are
relevant to free electrons, we develop a paradigm that can be used as a direct
test for quantum coherence in stochastic electrodynamics. Building on our
previous results Huang2013 ; Huang2015 , we devise a laser excitation scheme
to promote a ground-state quantum harmonic oscillator to a squeezed
Schrödinger cat state Kienzler ; Etesse ; Huang . The Schrödinger cat state of
a harmonic oscillator is an analogy to the electron double-slit state,
$|L\rangle+|R\rangle$, where $|L\rangle$ and $|R\rangle$ indicate the left and
right electron slit states in the position space Bach ; Tonomura . Comparing
the QM probability distribution with that of SED harmonic oscillators, we
observe some interesting similarities. Nevertheless, the interference pattern
is missing in the SED probability distribution.

##  II Kapitza-Dirac force on harmonic oscillators

Let us consider the setup in FIG. 1. Two counter-propagating laser fields
propagate along the $x$-axis, and the electric fields are linearly polarized
along the $z$-axis,

$\begin{split}\mathbf{E}_{1}=&A_{1}\omega_{1}\cos{(k_{1}x-\omega_{1}t)}\hat{\epsilon}_{z}\\\
\mathbf{E}_{2}=&-A_{2}\omega_{2}\cos{(k_{2}x+\omega_{2}t)}\hat{\epsilon}_{z},\end{split}$
(1)

where $\hat{\epsilon}_{z}$ is the unit vector along the $z$-axis,
$k_{1,2}=\omega_{1,2}/c$ are the wave numbers, and $A_{1,2}$ are amplitudes of
the corresponding vector potentials
$\mathbf{A}_{1}=A_{1}\sin{(k_{1}x-\omega_{1}t)}\hat{\epsilon}_{z}$ and
$\mathbf{A}_{2}=A_{2}\sin{(k_{2}x+\omega_{2}t)}\hat{\epsilon}_{z}$. The
electric and magnetic components of the combined laser field are

$\begin{split}E_{z}=A_{1}\omega_{1}\cos{(k_{1}x-\omega_{1}t)}-A_{2}\omega_{2}\cos{(k_{2}x+\omega_{2}t)}\\\
B_{y}=-A_{1}k_{1}\cos{(k_{1}x-\omega_{1}t)}-A_{2}k_{2}\cos{(k_{2}x+\omega_{2}t)}.\end{split}$
(2)

Assuming that the particle is free in the direction of the electric field
(i.e. the $z$-axis in FIG. 1), the cross terms between the electric and
magnetic components can give rise to a spatially modulated force in the
direction of wave propagation (i.e. the $x$-axis in FIG. 1). Herein we term
this force the Kapitza-Dirac (KD) force. The KD force can be derived using the
following equations of motion,

$\left\\{\begin{array}[]{l}\displaystyle m\frac{dv_{z}}{dt}=qE_{z}\\\ \\\
\displaystyle F^{(\textrm{KD})}_{x}=qv_{z}B_{y}\end{array}\right.,$ (3)

where $m$ and $q$ are mass and charge of the particle. Now, there are two
scenarios: (1) if frequencies of the two laser fields are identical
($\omega_{1}=\omega_{2}$), the KD force will be constant in time, which is
also known as the pondermotive force Batelaan2007 , (2) if the laser
frequencies are different ($\omega_{1}\neq\omega_{2}$), the KD force will
oscillate in time with the sum and difference frequencies,
$\omega_{1}+\omega_{2}$ and $\omega_{1}-\omega_{2}$. If the charge particle is
confined by a harmonic potential $U(x)=m\omega_{0}^{2}x^{2}/2$ with a resonant
frequency $\omega_{0}=(\omega_{1}-\omega_{2})/2$, the KD force that can
resonantly drive the harmonic oscillator will be (see derivation in Appendix)

$F_{\textrm{KD}}=\frac{q^{2}A_{1}A_{2}}{m}\left(\frac{k_{1}+k_{2}}{2}\right)\sin{\left((k_{1}+k_{2})x\right)}\cos{\left((\omega_{1}-\omega_{2})t\right)}.$
(4)

Accordingly, the corresponding time-varying KD potential is

$U_{\textrm{KD}}=\frac{q^{2}A_{1}A_{2}}{2m}\cos{\left((k_{1}+k_{2})x\right)}\cos{\left((\omega_{1}-\omega_{2})t\right)}.$
(5)

We note that at any given time $t=t_{0}$ a trapping site in the KD potential
$U_{\textrm{KD}}(x,t_{0})$ (i.e. a minimum in the potential) will turn into an
unstable point after a quarter of the natural period
$T_{0}/4=\pi/|\omega_{1}-\omega_{2}|$, where $T_{0}=2\pi/\omega_{0}$, since
the potential polarity is reversed,

$U_{\textrm{KD}}(x,t_{0}+T_{0}/4)=-U_{\textrm{KD}}(x,t_{0}).$ (6)

This feature will later be used to coherently split the ground-state
wavepacket of a quantum oscillator.

##  III Generation of squeezed Schrödinger cat states

The KD effect for quantum harmonic oscillators can be modeled by adding the KD
potential in Eq. (5) to the unperturbed oscillator Hamiltonian. We replace the
continuous-wave laser fields in Eq. (1) with pulsed fields,

$\begin{split}\mathbf{E}_{1}(x,t)&=A_{1}\omega_{1}\cos{(k_{1}x-\omega_{1}t)}e^{-(t/\tau)^{2}}{\hat{\epsilon}_{z}}\\\
\mathbf{E}_{2}(x,t)&=-A_{2}\omega_{2}\cos{(k_{2}x+\omega_{2}t)}e^{-(t/\tau)^{2}}{\hat{\epsilon}_{z}},\end{split}$
(7)

where $\tau$ is the pulse duration, in order to avoid indefinite sequential
excitation of the oscillator’s ladder levels. Given the appropriate pulse
amplitudes and durations, the final population distribution can have a peaked
structure. The quantum Hamiltonian is thus

$\hat{H}=\frac{m\omega_{0}}{2}\hat{x}^{2}+\frac{\hat{p}^{2}}{2m}+U_{KD}(\hat{x},t)e^{-2(t/\tau)^{2}}.$
(8)

The difference frequency of the laser fields is twice the oscillator’s
resonant frequency $\omega_{1}-\omega_{2}=2\omega_{0}$, so the ground-state
oscillator will be parametrically excited to the even-symmetry states
$|n=2k\rangle$, which is a prerequisite for cat state generation because a cat
state has an even symmetry. We obtained the QM result by numerically solving
the Schödinger equation as in Huang2015 . The oscillator’s parameters are
$m=9.11\times 10^{-35}$ (kg), $q=1.60\times 10^{-19}$ (C), and
$\omega_{0}=10^{16}$ (rad/s). The laser parameters are chosen to be
$\omega_{1}=2.3\,\omega_{0}$, $\omega_{2}=0.3\,\omega_{0}$, $\tau=5\times
10^{-15}$ (s), and $A_{1}=A_{2}=4.5\times 10^{-8}$ (V$\cdot$s/m). The mass is
chosen at this unusual value in order to make the computational time for SED
simulation more manageable Huang2013 . In FIG. 2(b) the probability trajectory
of the excited state is shown. Upon excitation the ground-state wavepacket is
coherently split into two wavepackets, so the oscillator is in a superposition
of two macroscopically distinct states. The two wavepackets oscillate back and
forth in the harmonic potential. As they recombine, interference fringes
emerge in the probability distribution due to the quantum coherence between
the two wavepackets. The quantum coherence is readily illustrated by the
fringe structure in the oscillator’s Wigner function as shown in FIG. 2(a). We
note that one of the two quadrature uncertainties of each wavepacket,
$\sigma_{x}\equiv\sqrt{\langle x^{2}\rangle-\langle x\rangle^{2}}$ (or
$\sigma_{p}\equiv\sqrt{\langle p^{2}\rangle-\langle p\rangle^{2}}$), is
smaller than the ground state uncertainty
$\Delta_{x}=\sqrt{\hbar/2m\omega_{0}}$ (or $\Delta_{p}=\sqrt{\hbar
m\omega_{0}/2}$), while their product remains the same,
$\sigma_{x}\sigma_{p}=\hbar/2$, at all times (see FIG. 2(a)). This implies
that the state generated by the laser excitation is a squeezed Schrödinger cat
state (that is, a superposition of two displaced squeezed states with opposite
phases) Etesse ; Huang .

For SED oscillators, we first prepare the classical ensemble in a ground state
with the $x$ and $p$ probability distributions identical to those of a quantum
oscillator (see details in Huang2013 ). Under excitation of the same laser
pulses, the equation of motion for the SED harmonic oscillator is

$m\frac{d^{2}x}{dt^{2}}=-m\omega_{0}^{2}x-m\Gamma\omega_{0}^{2}\frac{dx}{dt}+qE^{(x)}_{vac}(x,t)+F_{\textrm{KD}}(x,t)e^{-2(t/\tau)^{2}},$
(9)

where
$\displaystyle\Gamma\equiv\frac{2q^{2}}{3mc^{3}}\frac{1}{4\pi\epsilon_{0}}$ is
the radiation damping coefficient. The zero-point field $E^{(x)}_{vac}(x,t)$
is configured according to Huang2013 . Apart from the damping term
$m\Gamma\omega_{0}^{2}dx/dt$ and the coupling to the zero-point field
$E^{(x)}_{vac}(x,t)$, Eq. (9) is formally equivalent to a quantum Heisenberg
equation derived from the Hamiltonian in Eq. (8). This suggests that the
excitation dynamics in SED should be identical to QM assuming (1) the pulse
duration $\tau$ is much shorter than the damping time
$\tau_{d}=2/\Gamma\omega_{0}^{2}$, and (2) the KD force is much stronger than
the fluctuating force from the zero-point field Huang2013 ,

$\frac{q^{2}A_{1}A_{2}}{m}\left(\frac{k_{1}+k_{2}}{2}\right)\gg\frac{q}{2\pi}\sqrt{\frac{\hbar\Gamma\omega_{0}^{5}}{\epsilon_{0}c^{3}}}.$
(10)

In our simulation these two conditions are satisfied, and the excitation
dynamics in SED and QM are the same as shown in FIG. 3(a), where the time
evolutions of the expectation value of the QM energy and the ensemble average
of the SED energy are compared. The two energy trajectories stay overlapped
through most of the pulse and only deviate at the end of the excitation.
Furthermore, the SED and QM energy distributions have similar shapes (see FIG.
3(b)), despite that the QM distribution is discrete and the SED distribution
is continuous.

The probability trajectory of the SED oscillator ensemble along with its phase
space distribution are shown in FIG. 2(c)(d). The ensemble particle number is
$N_{p}=3\times 10^{4}$. The parametric interaction between the SED oscillator
and the laser fields were simulated using the same method as in Huang2015 .
Like the QM oscillator, upon excitation the ground-state SED probability
distribution also splits into two sub-ensembles that follow two distinct
sinusoidal trajectories. The initial phase spectrum of the zero-point field,
which is random and considered as the“hidden variable”, determines which
trajectory a particle will follow in each trial. A detailed comparison between
SED and QM probability distributions is given in FIG. 4 for two time points:
(1) when the two QM wavepackets are separated, (2) when they recombine.
Although there are some overall similarities between the SED and QM
distributions, there are no interference fringes in the SED probability
distribution.

Figure 3: The QM and SED energy distributions after laser excitation. (a) The
ensemble average of the SED oscillator energy (red line) is compared with the
expectation value of the QM oscillator energy (blue line) during the course of
laser excitation. The shaded area (red) represents the excitation laser pulse
$e^{-(t/\tau)^{2}}$. The pulse duration is $\tau=5\times 10^{-15}$(s) which is
much smaller than the SED oscillator damping time $\tau_{d}\approx 3.2\times
10^{-13}$ (s). Therefore, damping has no significant effect on the
oscillator’s dynamics during the excitation process. (a) Energy distributions
of QM and SED oscillators are compared at $t=t_{m}$ when the SED energy
reaches its maximum and the radiation damping starts to dominate. The QM
energy distribution (blue bar) is discrete and occupies only even energy
levels $E_{2k}=E_{0}(2k+1/2)$, where $E_{0}=\hbar\omega_{0}$. The SED energy
distribution (red bar) is continuous but has a similar width and average value
as the QM distribution. This indicates that the excitation process is
identical for QM and SED oscillators.

Figure 4: Comparison between the QM and SED probability distributions at $t=0$
and $t=T_{0}/4$ in FIG. 2. (a) The agreement between QM and SED probability
distributions is good when the two macroscopically distinct QM wavepackets are
well-separated by a peak-to-peak distance $a$ beyond the ground state width
$\Delta_{x}$, $a\gg\Delta_{x}$. (b) After a quarter of the natural period
$T_{0}/4$, the two QM wavepackets recombine in the harmonic potential.
Interference fringes appear in the QM distribution (blue line) but not in the
SED distribution (red dots). The fringe periodicity $\lambda$ is determined by
the wavepacket separation $a$ in (a) through the relation
$\lambda=2\pi\hbar/m\omega_{0}a$, which resembles the well-known double-slit
diffraction formula. The SED probability distribution captures the outline of
the QM distribution as if the quantum coherence between the two QM wavepackets
is lost.

##  IV Discussion and Conclusion

While the zero-point electromagnetic field only introduces small radiative
corrections to nonrelativistic QM, such as Lamb shifts, it drastically changes
the particle dynamics in classical mechanics and leads to the reproduction of
QM effects in some classical systems Boyer1973 ; Boyer1975-2 ; Huang2013 ;
Boyer1980 ; Boyer1984 ; Boyer2017 ; Santos ; Huang2015 . Our work aims to
investigate to what extent such a classical theory can reproduce QM features
by comparing results obtained from SED, Eq. (9), with those obtained from the
QM Hamiltonian, Eq. (8). The qualitative difference between the probability
distributions of QM and SED in FIG. 2 and FIG. 4(b) establishes that SED in
its traditional form does not support physical effects that involve quantum
coherence Milonni ; Boyer1975-1 . The squeezed Schrödinger cat state used in
this work is an analogy to the electron double-slit experiment Bach . The
peak-to-peak separation $a$ between the two wavepackets in FIG. 4(a)
determines the fringe periodicity $\lambda$ in FIG. 4(b) through the relation
$\lambda=2\pi\hbar/m\omega_{0}a$, which mimics the double-slit diffraction
formula. Our analysis provides evidence that coherence-like behavior is absent
in SED, and thus we predict that SED electron double-slit diffraction, if ever
calculated, will not show fringes.

On the other hand, our result helps to establish the validity range of SED. We
note that the partial agreement between the SED and QM results stems from the
formal resemblance between the SED equation of motion and the QM Heisenberg
equation assuming that the laser excitation pulses satisfy certain criteria.
While $x$ and $p$ are independent dynamic variables in SED, the canonical
commutation relation $[\hat{x},\hat{p}]=i\hbar$ makes them a Fourier pair in
QM. This difference makes the distinction between QM and SED in terms of
quantum coherence. Therefore, SED may be seen as the decoherence limit of QM
Sonnentag . Although the zero-point field bestows a special phase relation
between $x$ and $p$ of a SED harmonic oscillator, which leads to the quantum
ground-state distributions Boyer1975-2 ; Huang2013 , the phase relation serves
only as initial conditions and does not affect the dynamical evolution of $x$
and $p$ during laser excitation. Therefore, we speculate that any mechanism or
theoretical operation that restores (or deteriorates) the Fourier relation
between $x$ and $p$ for SED (or QM), will make the proper decoherence theory
that bridges the gap between SED and QM.

##  V Acknowledgement

The authors thank A. M. Steinberg, C. Monroe, and P. W. Milonni, for advice.
W. C. Huang wishes to give special thanks to Yanshuo Li for helpful
discussions. This work utilized high-performance computing resources from the
Holland Computing Center of the University of Nebraska. We gratefully
acknowledge funding support from NSF PHY-1602755.

*

##  Appendix A Derivation of the resonant Kapitza-Dirac force

In this appendix we derive the KD force in Eq. (4) from Eq. (3) using the
electric and magnetic components of the combined laser field given in Eq. (2).
First, we solve the velocity $v_{z}$ by integrating the equation of motion
$mdv_{z}/dt=qE_{z}$,

$v_{z}(t)=-\frac{q}{m}\left(A_{1}\sin{(k_{1}x-\omega_{1}t)}+A_{2}\sin{(k_{2}x+\omega_{2}t)}\right).$
(11)

Substituting $v_{z}(t)$ to $F_{x}=qv_{z}B_{y}$, we obtain the Lorentz force

$\begin{split}F_{x}=&\frac{q^{2}}{m}\left(A_{1}\sin{(k_{1}x-\omega_{1}t)}+A_{2}\sin{(k_{2}x+\omega_{2}t)}\right)\\\
&\times\left(A_{1}k_{1}\cos{(k_{1}x-\omega_{1}t)}+A_{2}k_{2}\cos{(k_{2}x+\omega_{2}t)}\right).\end{split}$
(12)

We can see four frequency components in the Lorentz force by expanding Eq.
(12),

$\begin{split}F_{x}=&\frac{q^{2}}{2m}\left[A_{1}^{2}k_{1}\sin{(2k_{1}x-2\omega_{1}t)}\right.\\\
&+A_{1}A_{2}k_{1}\left(\sin{((k_{1}+k_{2})x-(\omega_{1}-\omega_{2})t)}-\sin{((k_{1}-k_{2})x-(\omega_{1}+\omega_{2})t)}\right)\\\
&+A_{1}A_{2}k_{2}\left(\sin{((k_{1}+k_{2})x-(\omega_{1}-\omega_{2})t)}+\sin{((k_{1}-k_{2})x-(\omega_{1}+\omega_{2})t)}\right)\\\
&+\left.A_{1}^{2}k_{2}\sin{(2k_{2}x+2\omega_{2}t)}\right].\end{split}$ (13)

Because the frequency components $2\omega_{1}$, $2\omega_{2}$, and
$\omega_{1}+\omega_{2}$ are not integer multipliers of $\omega_{0}$, we can
drop these terms and keep only the parametrically resonant term
$\omega_{1}-\omega_{2}=2\omega_{0}$,

$F_{x}\approx\frac{q^{2}}{2m}A_{1}A_{2}(k_{1}+k_{2})\sin{((k_{1}+k_{2})x-(\omega_{1}-\omega_{2})t)}.$
(14)

The force has a traveling wave profile which can be decomposed into two
standing-wave components with even and odd symmetries,

$\sin{((k_{1}+k_{2})x-(\omega_{1}-\omega_{2})t)}=\sin{((k_{1}+k_{2})x)}\cos{((\omega_{1}-\omega_{2})t)}-\cos{((k_{1}+k_{2})x)}\sin{((\omega_{1}-\omega_{2})t)}.$
(15)

The force needs to have a potential profile with even symmetry in order to
resonantly drive the oscillator from the ground state with an even frequency
($\omega_{1}-\omega_{2}=2\omega_{0}$). This implies that the resonant KD force
should an odd symmetry, thus it takes the form

$F_{KD}=\frac{q^{2}}{2m}A_{1}A_{2}(k_{1}+k_{2})\sin{((k_{1}+k_{2})x)}\cos{((\omega_{1}-\omega_{2})t)}.$
(16)

## References

  * (1) L. de la Peña and A. M. Cetto, The Quantum Dice: an Introduction to Stochastic Electrodynamics, Springer (1996).
  * (2) L. de la Peña, A. M. Cetto, A. Valdes-Hernandez, Emerging Quantum: the physics behind quantum mechanics, Springer (2015).
  * (3) G. Cavalleri et al., Front. Phys. China 5, 107 (2010).
  * (4) P. W. Milonni, The Quantum Vacuum: An Introduction to Quantum Electrodynamics, Academic Press (1994).
  * (5) T. H. Boyer, Phys. Rev. D 11, 790 (1975).
  * (6) T. H. Boyer, Phys. Rev. A 7, 1832 (1973).
  * (7) T. H. Boyer, Phys. Rev. D 11, 809 (1975).
  * (8) W. C. Huang and H. Batelaan, J. Comput. Methods Phys. 2013, 308538 (2013).
  * (9) T. H. Boyer, Phys. Rev. A 21, 66 (1980).
  * (10) T. H. Boyer, Eur. J. Phys. 37, 065102 (2016).
  * (11) T. H. Boyer, Phys. Rev. D 29, 1096 (1984).
  * (12) T. H. Boyer, Eur. J. Phys. 38, 045101 (2017).
  * (13) R. Blanco, H. M. Franca, and E. Santos, Phys. Rev. A 43, 693 (1991).
  * (14) H. E. Puthoff, Phys. Rev. D 35, 3266 (1987).
  * (15) D. C. Cole and Y. Zou, Phys. Lett. A 317, 14 (2013).
  * (16) T. M. Nieuwenhuizen and M. T. P. Liska, Phys. Scr. T165, 014006 (2015).
  * (17) T. M. Nieuwenhuizen, Entropy 18, 135 (2016).
  * (18) W. C. Huang and H. Batelaan, Found. Phys. 45, 333 (2015).
  * (19) R. Bach, D. Pope, S.-H. Liou, and H. Batelaan, New J. Phys. 15, 033018 (2013).
  * (20) J. Avendaño and L. de la Peña, Phys. Rev. E 72, 066605 (2004).
  * (21) A. F. Kracklauer, Phys. Essays 5, 226 (1992).
  * (22) N. Wolchover, Famous Experiment Dooms Alternative to Quantum Weirdness, Quanta Magazine Oct. 11 (2018).
  * (23) Y. Couder and E. Fort, Phys. Rev. Lett. 97, 154101 (2006).
  * (24) A. Andersen, J. Madsen, C. Reichelt, S. R. Ahl, B. Lautrup, C. Ellegaard, M. T. Levinsen, and T. Bohr, Phys. Rev. E 92, 013006 (2015).
  * (25) T. Bohr, A. Andersen, and B. Lautrup, Recent Advances in Fluid Dynamics with Environmental Applications, pp. 335-349, Springer (2016).
  * (26) H. Batelaan, E. Jones, W. C. Huang, R. Bach, J. Phys.: Conf. Ser. 701, 012007 (2016).
  * (27) G. Pucci, D. M. Harris, L. M. Faria, and J. W. M. Bush, J. Fluid Mech. 835, 1136 (2018).
  * (28) D. Kienzler, C. Flühmann, V. Negnevitsky, H.-Y. Lo, M. Marinelli, D. Nadlinger, and J. P. Home, Phys. Rev. Lett. 116, 140402 (2016).
  * (29) J. Etesse, M Bouillard, B. Kanseri, and R. Tualle-Brouri, Phys. Rev. Lett. 114, 193602 (2015).
  * (30) K. Huang, H. Le Jeannic, J. Ruaudel, V. B. Verma, M. D. Shaw, F. Marsili, S. W. Nam, E. Wu, H. Zeng, Y.-C. Jeong, R. Filip, O. Morin, and J. Laurat, Phys. Rev. Lett. 115, 023602 (2015).
  * (31) A. Tonomura, J. Endo, T. Matsuda, T. Kawasaki, and H. Ezawa, Am. J. Phys, 57, 117 (1989).
  * (32) H. Batelaan, Rev. Mod. Phys. 79, 929 (2007).
  * (33) P. Sonnentag and F. Hasselbach, Phys. Rev. Lett. 98, 200402 (2007)