$\sqrt{N}\cdot\operatorname{\mathds{E}}\exp\left(\left\|\Psi^{y_{N}}\right\|_{2}^{2}-N\theta\right)\cdot\mathds{1}_{\left\|\Psi^{y_{N}}\right\|_{2}^{2}\leqslant
N\theta}$

converges to a constant. We have that

$\exp\left(\left\|\Psi^{y_{N}}\right\|_{2}^{2}-N\theta\right)=\int^{\infty}_{0}e^{-s}\cdot\mathds{1}_{s\geqslant
N\theta-\left\|\psi^{y_{N}}\right\|_{2}^{2}}\cdot ds,$

and thus

$\operatorname{\mathds{E}}\exp\left(\left\|\Psi^{y_{N}}\right\|_{2}^{2}-N\theta\right)\cdot\mathds{1}_{\left\|\Psi^{y_{N}}\right\|_{2}^{2}\leqslant
N\theta}=\int_{0}^{\infty}e^{-s}\cdot\mathds{P}\left(s\leqslant
N\theta-\left\|\Psi\right\|_{2}^{2}\leqslant 0\right)ds.$

The problem now comes down to analyzing

$\sqrt{N}\cdot\mathds{P}\left(s\leqslant
N\theta-\left\|\Psi^{y_{N}}\right\|_{2}^{2}\leqslant 0\right).$

Note that we may equivalently consider

$\displaystyle\sqrt{N}\cdot\mathds{P}\left(\frac{s}{\sqrt{N}}\leqslant\frac{N\theta-\left\|\Psi^{y_{N}}\right\|_{2}^{2}}{\sqrt{N}}\leqslant
0\right).$ (68)

Recall that $\left\|\Psi^{y_{N}}\right\|$ is expressible as a sum of
independent random variables, we have that

$\left\|\Psi^{y_{N}}\right\|_{2}^{2}\stackrel{{\scriptstyle\mathrm{d}}}{{=}}\sum_{\boldsymbol{k}\in[n]^{d}}\frac{1}{y_{N}+\lambda_{\boldsymbol{k}}}X_{\boldsymbol{k}}.$

Further, $N\theta-\left\|\Psi^{y_{N}}\right\|_{2}^{2}$ is a centered sum by
definition, and each random variable in the sum has a uniformly bounded third
moment. The Berry-Esseen Theorem then tells us that (68) is bounded above by a
constant. The local limit theorem yields convergence. $\blacksquare$

Theorem 8.1 gives us a means of evaluating probabilities of events with
respect to the reference measure, in terms of the MGFF, where of course we pay
a penalty due to the mass restriction. In the next sequence, we will examine
the effect of exponentially tilting the massive free field with respect to the
nonlinearity.

######  Theorem 8.2.

Let $\Psi^{y}$ be distributed according to the massive Gaussian Free Field on
$\mathds{T}^{d}$ with mass parameter $y$. Then we have a constant $C$
depending on $\nu$, $y$ and $d$ such that

$\operatorname{\mathds{E}}\exp\left(\frac{2\nu
N^{\frac{1-p}{2}}}{p+1}\left\|\Psi^{y}\right\|^{p+1}_{p+1}\right)\cdot\mathds{1}_{\left\|\Psi\right\|_{2}^{2}\leqslant
N}\leqslant e^{C}$

The mass constraint is a significant obstacle, as it is not expressible as a
product, nor is it smooth for us to apply some standard techniques to deal
with Gaussian integration. We resolve both these issues with the following
upper bound. We introduce the following auxilary function in order to deal
with this issue:

$\displaystyle h(x):=\frac{2x^{p+1}}{1+x^{p-1}}.$ (69)

Observe that when $x<1$, we have that $x^{p+1}<h(x)$, and when $x>1$,
$h(x)<2x^{2}$. Crucially, $h$ is smoothly differentiable.

######  Lemma 8.3.

For $p<5+\sqrt{32}$ we have that

$h^{\prime\prime}(x)+x^{-1}\cdot h^{\prime}(x)>0$

###### Proof.

It is equivalent to verify that $(xh^{\prime}(x))^{\prime}>0$. We merely carry
out the calculation. We have

$(xh^{\prime}(x))=\frac{2(p+1)x^{p}}{1+x^{p-1}}-\frac{2(p-1)x^{2p-1}}{(1+x^{p-1})^{2}}+\frac{(p-1)(p+1)x^{p}}{(1+x^{p-1})^{2}}-\frac{2(p-1)^{2}x^{2p-1}}{(1+x^{p-1})^{3}}.$

Adding and simplifying, we find that the numerator is given by

$x^{p}\left(4x^{2p-2}+(3+6p-p^{2})x^{p-1}+(p+1)^{2}\right),$

This equation has no real roots when $1<p<5+\sqrt{32}.$ $\blacksquare$

Merely by the definition of $h$, it follows that

$\displaystyle\operatorname{\mathds{E}}\exp\left(\frac{2\nu
N^{\frac{1-p}{2}}}{p+1}\left\|\Psi\right\|^{p+1}_{p+1}\right)\cdot\mathds{1}_{\left\|\Psi\right\|_{2}^{2}\leqslant
N}\leqslant\operatorname{\mathds{E}}\prod_{{\boldsymbol{x}}\in\mathds{T}_{n}^{d}}\exp\left(\frac{2\nu
N}{p+1}h\left(\frac{|\Psi_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right)$ (70)

The intuition now is to work with the right side of (70) and attempt to use
the decay of correlations to separate as a product of expectations,
essentially comparing to the i.i.d. situation. We will first establish a
counterpart of Theorem 8.2 for an i.i.d. standard complex Gaussian vector
$\\{\Phi_{{\boldsymbol{x}}}\\}_{{\boldsymbol{x}}\in\mathds{T}^{d}}$, and then
use Gaussian interpolation to compare the expectations under the MGFF and
i.i.d. cases.

######  Lemma 8.4.

Let $\Phi$ be an i.i.d standard complex Gaussian vector. Then for $p>3$ and

$\nu<(p+1)/8,$

we have a constant $C$ depending on $\nu$ and $p$ such that

$\operatorname{\mathds{E}}\exp\left(\frac{2\nu
N}{p+1}h\left(\frac{|\Phi_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right)\leqslant
e^{C}.$

###### Proof.

We may immediately express the expectation as a product using the
independence, and noting that $|\Phi_{x}|^{2}$ is distributed according to
exponential $1/2$, we have

$\displaystyle\operatorname{\mathds{E}}\exp\left(\frac{2\nu
N}{p+1}h\left(\frac{|\Phi_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right)=\left(\frac{1}{2}\int_{0}^{\infty}\exp\left(\frac{2\nu
N}{p+1}h\left(\sqrt{N^{-1}t}\right)\right)\cdot e^{-t/2}dt\right)^{N}.$ (71)

We split the integral into separate regions and bound their contributions
individually. Let $\alpha=1-4/(p+1)$. We will consider the intervals
$[0,N^{\alpha})$ and $[N^{\alpha},\infty)$. This is where $p>3$ becomes
important, as we want $N^{\alpha}\to\infty$ as $N\to\infty$. For the first
interval,

$\displaystyle\frac{1}{2}\int_{0}^{N^{\alpha}}\exp\left(\frac{2\nu
N}{p+1}h\left(\sqrt{\frac{t}{N}}\right)\right)e^{-\frac{t}{2}}dt$
$\displaystyle\leqslant\exp\frac{4\nu}{(p+1)N}\cdot\frac{1}{2}\int_{0}^{N^{\alpha}}e^{-t/2}dt$
$\displaystyle\leqslant\exp\frac{2\nu(a+1)}{a(p+1)N}$

In this bound, we have used the fact that on the interval $[0,N^{\alpha-1})$,
$h(x)\leqslant 2x^{p+1}$. Now take a fixed $\varepsilon>0$, we know that for
$N$ sufficiently large,

$\exp\frac{2\nu(a+1)}{a(p+1)N}\leqslant
1+(1+\varepsilon)\frac{2\nu(a+1)}{a(p+1)N}.$

The case to consider is $t\in[N^{\alpha},\infty)$. On this interval, we may
use the bound $h(x)\leqslant(a+1)x^{2}$ to obtain

$\int_{N^{\alpha}}^{N}\exp\left(\frac{2\nu
N}{p+1}h\left(\sqrt{\frac{t}{N}}\right)\right)e^{-\frac{t}{2}}dt\leqslant
N\exp\biggl{(}\bigl{(}-\frac{1}{2}+\frac{4\nu}{p+1}\bigr{)}N^{\alpha-1}\biggr{)}.$

This decays to zero rapidly as $N\to\infty$, so long as
$\nu<\frac{1}{8}(p+1).$ Thus, combining we obtain that the following is an
upper bound for (71)

$\left(1+\frac{6\nu(1+\varepsilon)}{(p+1)N}\right)^{N}\leqslant\exp\left(\frac{6\nu(1+\varepsilon)}{p+1}\right).$

$\blacksquare$

To use Lemma 8.4 to prove Theorem 8.2, we need one more ingredient to compare
the expectation w.r.t. the massive Gaussian Free Field and the i.i.d. case. We
may use Gaussian interpolation and correlation decay to accomplish this.

######  Lemma 8.5.

Let $\Psi^{y}$ denote the MGFF on $\mathds{T}^{d}_{n}$, let $\Phi$ denote a
standard complex Gaussian vector, and let $C_{y}$ be a constant such that

$\displaystyle C_{y}\geqslant
2\sum_{{\boldsymbol{x}}_{2}}(y-\Delta)^{-1}_{{\boldsymbol{x}}_{1}{\boldsymbol{x}}_{2}}.$
(72)

Then we have that

$\operatorname{\mathds{E}}\prod_{x\in\mathds{T}_{d}}\exp\left(\frac{2\nu
N}{p+1}h\left(\frac{|\Psi^{y}_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right)\leqslant\left(\operatorname{\mathds{E}}\exp\left(\frac{2\nu
N}{p+1}h\left(\frac{\sqrt{C_{y}}|\Phi_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right)\right)^{N}.$

Lemma 8.5 hinges on the application of Gaussian interpolation. We state the
version used below, easily adapted from the real-valued case.

######  Lemma 8.6.

Let $\Psi^{0}$ and $\Psi^{1}$ be independent, centered, complex-valued
Gaussian processes on $\mathds{T}^{d}_{n}$, with covariance matrices $G_{0}$
and $G_{1}$ respectively. Let $r:\mathds{C}^{N}\to\mathds{R}$ be integrable
w.r.t. the laws of both $\Psi^{0}$ and $\Psi^{1}$. Let
$\Psi^{t}:=\sqrt{1-t}\cdot\Psi^{0}+\sqrt{t}\cdot\Psi^{1}$, and define
$R(t):=\operatorname{\mathds{E}}r(\Psi^{t})$. Then we have that

$R^{\prime}(t)=\sum_{{\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2}\in\mathds{T}^{d}_{n}}\bigl{(}G_{1}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})-G_{0}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})\bigr{)}\cdot\operatorname{\mathds{E}}\frac{\partial^{2}}{\partial\psi_{{\boldsymbol{x}}_{1}}\partial\psi_{{\boldsymbol{x}}_{2}}}r(\Psi^{t}).$

It is clear that $R(0)$ denotes the expectation of $r$ under the law of
$\Psi^{0}$, and $R(1)$ the expectation under the law of $\Psi^{1}$. Thus, if
$R^{\prime}$ can be shown to be nonnegative, it is clear that $R(1)\geqslant
R(0)$. We are now ready to prove Lemma 8.5, and then obtain Theorem 8.2 for
free.

###### Proof of Lemma 8.5.

We apply Lemma 8.6 with $\Psi^{0}=\Psi^{y}$, $\Psi^{1}=\sqrt{C_{y}}\ \Phi$ as
defined in (72), and

$r(\psi)=\exp\left(\frac{2\nu
N}{p+1}\sum_{{\boldsymbol{x}}\in\mathds{T}^{d}}h\left(\frac{|\psi_{{\boldsymbol{x}}}|}{\sqrt{N}}\right)\right).$

Observe that for any smooth function $f:\mathds{R}\to\mathds{R}$, we have

$(\exp(f))^{\prime}=\exp(f)\cdot f^{\prime}\text{ and
}(\exp(f))^{\prime\prime}=\exp(f)\cdot(f^{\prime\prime}+(f^{\prime})^{2}).$

In our cases, $\Re{\Psi}$ and $\Im{\Psi}$ are independent copies of one
another, and therefore any calculations done for the real part are exactly
replicated for the imaginary part. For convenience, we introduce the notation
$\mathfrak{J}$ which denotes either $\Re$ or $\Im$. With our choice of $r$, we
have

$\frac{\partial^{2}}{\partial\mathfrak{J}_{1}\psi_{{\boldsymbol{x}}_{1}}\partial\mathfrak{J}_{2}\psi_{{\boldsymbol{x}}_{2}}}r(\psi)=N\cdot\left(\frac{2\nu}{p+1}\right)^{2}\cdot
h^{\prime}\left(\frac{|\psi_{{\boldsymbol{x}}_{1}}|}{\sqrt{N}}\right)\cdot\frac{\partial|\psi_{\boldsymbol{x}_{1}}|}{\partial\mathfrak{J}_{1}\psi_{\boldsymbol{x}_{1}}}\cdot
h^{\prime}\left(\frac{|\psi_{{\boldsymbol{x}}_{2}}|}{\sqrt{N}}\right)\cdot\frac{\partial|\psi_{\boldsymbol{x}_{2}}|}{\partial\mathfrak{J}_{2}\psi_{\boldsymbol{x}_{2}}}\cdot
r(\psi).$

For the non mixed second derivative, we have

$\displaystyle\frac{\partial^{2}}{\partial\mathfrak{J}\psi_{{\boldsymbol{x}}}^{2}}r(\psi)$
$\displaystyle=N\cdot r(\psi)\cdot\left(\frac{2\nu}{p+1}\cdot
h^{\prime}\left(\frac{|\psi_{\boldsymbol{x}}|}{N}\right)\cdot\frac{\partial|\psi_{\boldsymbol{x}}|}{\partial\mathfrak{J}\psi_{\boldsymbol{x}}}\right)^{2}+r(\psi)\cdot\frac{2\nu}{p+1}\cdot\left(\frac{\partial|\psi_{\boldsymbol{x}}|}{\partial\mathfrak{J}\psi_{\boldsymbol{x}}}\right)^{2}h^{\prime\prime}\left(\frac{|\psi_{{\boldsymbol{x}}}|}{N}\right)$
$\displaystyle\qquad\qquad+N\cdot
r(\psi)\cdot\left(\frac{2\nu}{p+1}\right)h^{\prime}\left(\frac{|\psi_{{\boldsymbol{x}}}|}{N}\right)\cdot\frac{\partial^{2}|\psi_{\boldsymbol{x}}|}{\partial\mathfrak{J}\psi_{\boldsymbol{x}}^{2}}.$

We focus our attention on the latter two terms, that is

$\displaystyle
r(\psi)\cdot\frac{2\nu}{p+1}\cdot\left(\frac{\partial|\psi_{\boldsymbol{x}}|}{\partial\mathfrak{J}\psi_{\boldsymbol{x}}}\right)^{2}h^{\prime\prime}\left(\frac{|\psi_{{\boldsymbol{x}}}|}{N}\right)+\sqrt{N}\cdot
r(\psi)\cdot\left(\frac{2\nu}{p+1}\right)h^{\prime}\left(\frac{|\psi_{{\boldsymbol{x}}}|}{N}\right)\cdot\frac{\partial^{2}|\psi_{\boldsymbol{x}}|}{\partial\mathfrak{J}\psi_{\boldsymbol{x}}^{2}}$
(73)

Evaluating the sum of (73) over both cases of $\mathfrak{J}$, that is $\Re$ or
$\Im$, we obtain

$r(\psi)\cdot\left(\frac{\partial^{2}}{\partial|\psi_{\boldsymbol{x}}|^{2}}+\frac{1}{|\psi_{\boldsymbol{x}}|}\frac{\partial}{\partial|\psi_{\boldsymbol{x}}|}\right)h\left(\frac{|\psi_{\boldsymbol{x}}|}{\sqrt{N}}\right).$

Using Lemma 8.3, we may conclude that

$\frac{\partial^{2}}{\partial\Re\psi_{\boldsymbol{x}}^{2}}r(\psi)+\frac{\partial^{2}}{\partial\Im\psi_{\boldsymbol{x}}^{2}}r(\psi)-\left(\frac{\partial}{\partial\Re\psi_{\boldsymbol{x}}}r(\psi)\right)^{2}-\left(\frac{\partial}{\partial\Im\psi_{\boldsymbol{x}}}r(\psi)\right)^{2}\geqslant
0$

Combining the inequality $a^{2}+b^{2}>2ab$ with (72), we have that

$R^{\prime}(t)\geqslant\sum_{x_{1}}\left(C_{y}-\sum_{{\boldsymbol{x}}_{2}}(y-\Delta)^{-1}_{{\boldsymbol{x}}_{1}{\boldsymbol{x}}_{2}}\right)\operatorname{\mathds{E}}\left(\frac{\partial^{2}}{\partial\Re\psi_{{\boldsymbol{x}}}^{2}}+\frac{\partial^{2}}{\partial\Im\psi_{{\boldsymbol{x}}}^{2}}\right)r(\psi)\geqslant
0.$

This completes the proof. $\blacksquare$

All parts are established to prove Theorem 8.2.

###### Proof.

A detailed proof is omitted; all that is required is to combine Lemma 8.5 with
a rescaled version of Lemma 8.4. $\blacksquare$

##  9\. Proofs of Main Theorems

In the sections prior, we have all the pieces required to prove our main
results. In this section, we tie them all together. We begin with the proof of
free energy convergence. All that needs to be done is to combine the upper and
lower bounds found in Section 6

###  9.1. Proof of Theorem 2.3

Essentially, the proof is an immediate corollary of Lemmas 6.1 and 6.5. It
suffices to specify $\kappa_{N}$ and $s_{N}$, and provide a rate of
convergence for $\gamma_{N}$. The sequence $\kappa_{N}$ is chosen such that
the mass shells are of adequate thickness such that the required concentration
of mass of the Gaussian free field holds. A larger $\kappa_{N}$ yields a
better concentration bound but a worse rate of convergence. It is necessary
that

$\lim_{N\to\infty}\frac{\kappa_{N}}{\sqrt{m_{N}(2)N}}=\infty\text{ and
}\lim_{N\to\infty}\frac{\kappa_{N}}{N}=0.$

As for $s_{N}$, which governs the size of the concentrated region in the upper
bound, observe that $s_{N}$ needs to be chosen such that the corresponding $U$
will never contain nontrivial cycles, $s_{N}\to 0$ as $N\to\infty$, and the
corresponding $U$ are sufficiently small so as to allow concentration. It
suffices to consider

$\displaystyle s_{N}=\frac{1}{\log N}.$ (74)

Finally, as for $\gamma_{N}$, this depends on the optimizing values $a$. We
note that if $0$ is optimizing, then by Lemma 5.8, we have that

$\gamma_{N}\leqslant\frac{2d\kappa_{N}}{N\log N}.$

On the contrary, if we have a non trivial minimezer $a_{\star}$, then we also
know that $\nu a_{\star}>R_{p}$, and thus by Lemma 5.7 we have that

$\gamma_{N}\leqslant C_{1}(\nu)(\nu a_{\star})^{p^{2}-1}\exp(-C_{2}(\nu)\cdot
N^{1/d}).$

Thus, the product of all the error terms arising in Lemmas 6.1 and 6.5 can be
bounded above by

$e^{C_{1}{\kappa_{N}}}\cdot\left(1-C\cdot\frac{Nm_{N}(2)}{\kappa_{N}^{2}}\right)^{2}.$

Using the fact that $m_{N}(2)=N^{-1+4/d}$, we find that the best choice of
$\alpha$ for $\kappa_{N}=N^{\alpha}$ is given by

$\alpha=\frac{1}{3}+\frac{4}{3d}$

$\blacksquare$

Next, the characterization of the phase transition curve. These are
essentially the finishing touches required to extract Theorem 2.4 from Section
7.

###  9.2. Proof of Theorem 2.4

We address the four parts Theorem 2.4 individually.

  1. a)

Immediate corollary of Lemma 7.2.

  2. b)

This is exactly established in Lemma 7.4.

  3. c)

We begin with the bound and asymptotic for $\nu\downarrow R_{p}$. Indeed,
Lemma 7.3 may be rephrased as

$\theta_{c}(\nu)\leqslant C_{d}\cdot\frac{\nu}{\nu-R_{p}}.$

Further, note that we have a non trivial minimizer $a_{\star}$ iff we have an
$a$ such that $G(a)<0$. In particular, thus, if $(\theta,\nu)$ corresponds to
the dispersive phase, we must have that $G(a)\geqslant 0$ for $a\in[0,1)$. Now
consider the curve given by

$\frac{R_{p}}{\nu}+\frac{C_{d}}{\theta}=(1+\varepsilon)$

We know that the function $J$ has a derivative, moreover $J^{\prime}=0$ for
all $a<R_{p}$. As $\theta\uparrow\infty$,
$\nu\downarrow(1+\varepsilon)^{-1}\cdot R_{p}$. Thus we must have that along
this curve

$\lim_{\theta\to\infty}\theta\cdot J(a\nu)=0$

for all values of $a\in[0,1)$. Thus, along this curve as $\theta$ increases,
$(\nu_{\varepsilon},\theta_{\varepsilon})$ will eventually correspond to the
dispersive phase. Since choice of $\varepsilon$ is arbitrary, this verifies
that

$\lim_{\nu\downarrow R_{p}}(\nu-R_{p})\theta_{c}(\nu)=C_{d}R_{p}.$

Finally, we address the bound and the asymptotic as $\nu\uparrow\infty$. To do
this, recall the functions $\widehat{W}$ and $\widehat{J}$ introduced in (44)
and (29) respectively. We may then write

$\displaystyle\log{\theta}$
$\displaystyle+W(\theta(1-a))+\frac{\theta}{\nu}I(a\nu)$
$\displaystyle=\widehat{W}(\theta(1-a))+\theta
a\widehat{J}(a\nu)-\frac{2\theta}{p+1}a(a\nu\vee R_{p})^{(p-1)/2}-\log(1-a).$

We know that the functions $\widehat{W}$ and $\widehat{J}$ are bounded,
moreover, we may discard the $\log\theta$ as it is irrelevant to the phase
behavior (independent of $a$). Thus, it suffices to characterize the behavior
of

$G_{\theta,\nu}(a)=\widehat{W}(\theta(1-a))+\theta
a\widehat{J}(a\nu)-\frac{2\theta}{p+1}a^{(p+1/2)}\cdot\nu^{(p-1)/2}-\log(1-a).$

If we take $C\cdot\theta=\nu^{-(p-1)/2}$ for a constant $C$ then observe that
then the limit as $\nu\to\infty$ is non negative iff

$C\leqslant\xi_{p}(0)$

for $\xi_{p}$ defined in (18). This completes the proof.

  4. d)

It suffices to consider $\nu>R_{p}$ and $\theta>\theta_{c}$. This implies that
there is a $a_{\star}>0$ which is the smallest optimizer. Now, suppose that as
$\theta\downarrow\theta_{c}$, and we have that
$a_{\star}(\theta_{c})\downarrow 0$. Since $a_{\star}(\theta)$ continue to be
non trivial minimizers we must still have $I(a_{\star}\nu)<0$, a contradiction
as eventually $a_{\star}\nu<R_{p}$. $\blacksquare$

We conclude this section by finally proving Theorem 2.8. Essentially, we need
to combine Theorems 8.1 and 8.2.

###  9.3. Proof of Theorem 2.8

Let $\mathcal{A}\subset\ell^{2}(\mathds{T}^{d}_{n})$. Using the trivial fact
that $\tilde{Z}_{N}\geqslant Z_{\text{ref}}$, we may write

$\displaystyle\tilde{\mu}_{N}(A)\leqslant\frac{Z^{y_{N}}e^{N\theta
y_{N}}}{\sqrt{N}\cdot
Z_{\text{ref}}(\theta)}\operatorname{\mathds{E}}\left(\mathds{1}_{\mathcal{A}}\cdot\frac{\exp(\left\|\Psi^{y}\right\|_{2}^{2})}{e^{N\theta
y_{N}}}\cdot\exp\left(\frac{\nu(\theta
N)^{(1-p)/2}}{p+1}\left\|\Psi^{y_{N}}\right\|_{p+1}^{p+1}\right)\cdot\mathds{1}_{\left\|\Psi^{y}\right\|_{2}^{2}\leqslant
N\theta}\right).$ (75)

Two applications of Hölder’s inequality are all that are required now. Indeed,
let $\epsilon_{1}$ and $\epsilon_{2}$ be positive real numbers, such that
$(1+\epsilon_{1})\nu\theta^{(p+1)/2}$ satisfies the hypothesis of Theorem 8.2.
By the first application, we have that (75) can be further bounded above by
the product of the following

$\displaystyle\operatorname{\mathds{E}}\left(\exp((1+\epsilon_{1})\cdot\frac{\nu
N^{-(p-1)/2}}{p+1}\cdot\left\|\Psi^{y_{N}}\right\|_{p+1}^{p+1})\cdot\mathds{1}_{\left\|\Psi^{y}\right\|_{2}^{2}\leqslant
N\theta}\right)^{1/(1+\epsilon_{1})}$ (76)

and

$\displaystyle\operatorname{\mathds{E}}\left(\exp\left(\frac{1+\epsilon_{1}}{\epsilon_{1}}y_{N}\left\|\psi^{y_{N}}\right\|_{2}^{2}-\frac{1+\epsilon_{1}}{\epsilon_{1}}y_{N}N\theta\right)\cdot\mathds{1}_{\left\|\Psi^{y}\right\|_{2}^{2}\leqslant
N\theta}\right)^{\epsilon_{1}/(1+\epsilon_{1})}$ (77)

We, for convenience, denote

$\epsilon_{1,2}=\epsilon_{1}+\epsilon_{2}+\epsilon_{1}\epsilon_{2}.$

Applying Hölder’s inequality again to (77), we obtain the upper bound

$\displaystyle\operatorname{\mathds{E}}\left(\exp\left(\frac{1+\epsilon_{1,2}}{\epsilon_{1}\epsilon_{2}}y_{N}\left\|\psi^{y_{N}}\right\|_{2}^{2}-\frac{1+\epsilon_{1,2}}{\epsilon_{1}\epsilon_{2}}y_{N}N\theta\right)\cdot\mathds{1}_{\left\|\Psi^{y}\right\|_{2}^{2}\leqslant
N\theta}\right)^{\epsilon_{1}\epsilon_{2}/(1+\epsilon_{1,2})}\operatorname{\mathds{P}}\left(\Psi^{y}\in\mathcal{A}\right)^{\epsilon_{1}/(1+\epsilon_{1,2})}.$
(78)

Thus, by Theorem 8.2, we have that (76) is bounded above by a constant. By
Theorem 8.1, we know that (78) is bounded above by
$N^{-\epsilon_{1}\epsilon_{2}/(2+2\epsilon_{1,2})}$. Combining, we obtain that

$\tilde{\mu}_{N}(A)\leqslant C\cdot
N^{\frac{1+\epsilon_{1}+\epsilon_{2}}{2+2\epsilon_{1,2}}}\cdot\mathds{P}\left(\Psi^{y}\in\mathcal{A}\right)^{\epsilon_{2}/(1+\epsilon_{1,2})}$

Indeed, if $\mathds{P}(\Psi^{y_{N}}\in\mathcal{A})\leqslant N^{-\alpha}$ for
some $\alpha>0$, we have that

$\tilde{\mu}_{N}(\mathcal{A})\leqslant C\cdot
N^{\frac{1+\epsilon_{1}+(1-2\alpha)\epsilon_{2}}{2+2\epsilon_{1,2}}}$

We have flexibility in choosing $\epsilon_{2}$, and so long as $\alpha>1/2$,
we may select $\epsilon_{2}$ such that the numerator of the exponent is
negative. $\blacksquare$

##  10\. Discussion and Further Questions

###  10.1. On the Question of Multi-Soliton Solutions

The question of multi-soliton phases is closely related to having multiple
minimizers for the variational formula in the soliton phase; _i.e.,_ $0$ is
not a minimizer. We conjecture that this is impossible. Indeed, we know that
$I(a)$ is concave for large values of $a$. However, we need detailed behavior
of the $I$ function near $R_{p}$ to prove this.

###  10.2. Ergodicity

The behavior of the invariant measure under the dynamics is yet to be
explored. In [CK12], the corresponding analysis was possible due to the mass
of typical functions being concentrated at a single lattice site,
significantly simplifying computations. This is no longer possible here; for
us, the corresponding concentration occurs on a region of size $O(1)$. It
should be possible to consider the dynamics introduced in [OL] with the regime
of scaling considered in this article. In [WEI3], a hierarchy of local minima
for the Hamiltonian (5) was provided, which would be the collection of
metastable states. In addition, the Witten Laplacian approach in [LEB] can be
considered with restriction to the sphere instead of the ball, so as to remove
technicalities arising from carrying out Morse theoretic calculations on a
manifold with boundary.

###  10.3. Dimension two analysis

We have explicitly used the finiteness of $C_{d}$ for our maximum bounds. One
of the obstacles in extending our results to the two-dimensional case is that
$C_{2}=\infty$. We work with a massive field for most parts, so this issue
does not arise often. We anticipate that one can work around it. The more
serious obstacle is that the techniques used here do not yield adequate mass
concentration in two dimensions. This being said, the asymptotics is more
interesting in two dimensions; by choosing the correct scaling, it might be
possible to use the 2-D Gaussian Free Field characterization in [RAY] to
evaluate our scaling limit explicitly in the dispersive phase.

Acknowledgments. We would like to thank Gayana Jayasinghe, Gourab Ray, and
Arnab Sen for many useful discussions.

## References


$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{|}\tilde{c}_{10}^{\tau}\langle\Sigma^{\prime\prime}_{\tau}\rangle+{\langle
f\rangle\over\langle
g\rangle}i\tilde{b}_{8}^{\tau}\langle\Sigma^{\prime\prime(0)}_{\tau}\rangle+{\langle
f\rangle\over\langle
g\rangle}i\tilde{b}_{16}^{\tau}\langle\Sigma^{\prime\prime(2)}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}~{}\Big{]}~{}\bigg{\\}}.$
(120)

If we specialize to the case of interactions involving scalar or vector
mediators, then we can use Tables 7 and 8 to re-express Eq. (119) in terms of
the 20 relativistic LECs $\tilde{d}_{i}$. One finds

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\omega={G_{F}^{2}\over\pi}~{}{q_{\mathrm{eff}}^{2}\over
1+{q\over M_{T}}}~{}{4\pi\over
2j_{N}+1}\sum_{\tau=0,1}\sum_{\tau^{\prime}=0,1}$
$\displaystyle\bigg{\\{}~{}\sum_{J=0,2,...}^{\infty}\Big{[}~{}\Big{|}(\tilde{d}_{1}^{\tau}-{q\over
m_{L}}\tilde{d}_{9}^{\tau})\langle
M^{g}_{\tau}-M^{(2)f}_{\tau}\rangle+\tilde{d}_{5}^{\tau}\langle
M_{\tau}^{g}+M^{(2)f}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{|}(\tilde{d}_{3}^{\tau}+{q\over
m_{L}}\tilde{d}_{17}^{\tau})\langle
M^{g}_{\tau}-M^{(2)f}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+\Big{|}\tilde{d}_{13}^{\tau}\langle
M^{g}_{\tau}+M^{(2)f}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}~{}\Big{]}$
$\displaystyle+\sum_{J=1,3,...}^{\infty}\Big{[}~{}\Big{|}\tilde{d}_{15}^{\tau}\langle\Sigma^{\prime
g}_{\tau}+\Sigma^{\prime(0)f}_{\tau}\rangle-{q_{\mathrm{eff}}\over
m_{N}}(\tilde{d}_{5}^{\tau}+{q\over
m_{L}}\tilde{d}_{9}^{\tau})\langle\Delta_{\tau}^{g}\rangle+{q\over
2m_{N}}(\tilde{d}_{5}^{\tau}-2\tilde{d}_{6}^{\tau}+{q\over
m_{L}}(\tilde{d}_{9}^{\tau}-2\tilde{d}_{10}^{\tau}))\langle\Sigma^{\prime
g}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{|}{q\over
m_{L}}\tilde{d}_{19}^{\tau}\langle\Sigma^{\prime
g}_{\tau}+\Sigma^{\prime(2)f}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Big{|}\tilde{d}_{7}^{\tau}\langle\Sigma^{\prime
g}_{\tau}+\Sigma^{\prime(0)f}_{\tau}\rangle+{q\over
m_{L}}\tilde{d}_{11}^{\tau}\langle\Sigma^{\prime
g}_{\tau}+\Sigma_{\tau}^{\prime(2)f}\rangle-{q_{\mathrm{eff}}\over
m_{N}}\tilde{d}_{13}^{\tau}\langle\Delta_{\tau}^{g}\rangle+{q\over
2m_{N}}(\tilde{d}_{13}^{\tau}-2\tilde{d}_{14}^{\tau})\langle\Sigma_{\tau}^{\prime
g}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{|}{q\over
m_{L}}{q_{\mathrm{eff}}\over
m_{N}}\tilde{d}^{\tau}_{17}\langle\Delta_{\tau}^{g}\rangle+{q\over
m_{L}}{q\over
2m_{N}}(2\tilde{d}^{\tau}_{18}-\tilde{d}_{17}^{\tau})\langle\Sigma^{\prime
g}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}\Big{]}$
$\displaystyle+\sum_{J=2,4,...}^{\infty}\Big{[}~{}\Big{|}(\tilde{d}_{3}^{\tau}+{q\over
m_{L}}\tilde{d}_{17}^{\tau})\langle
M_{\tau}^{(1)f}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+\Big{|}\tilde{d}_{13}^{\tau}\langle
M_{\tau}^{(1)f}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+\Big{|}(\tilde{d}_{1}^{\tau}-\tilde{d}_{5}^{\tau}-{q\over
m_{L}}\tilde{d}_{9}^{\tau})\langle
M_{\tau}^{(1)f}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}\Big{]}$
$\displaystyle+\sum_{J=1,3,...}^{\infty}\Big{[}~{}\Big{|}\tilde{d}_{15}^{\tau}\langle\Sigma^{\prime\prime
g}_{\tau}-\Sigma^{\prime\prime(0)f}_{\tau}\rangle-{q\over
2m_{N}}(\tilde{d}_{4}^{\tau}-{2q\over
m_{L}}\tilde{d}^{\tau}_{20})\langle\Sigma^{\prime\prime
g}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+\Big{|}{q\over
m_{L}}\tilde{d}_{19}^{\tau}\langle\Sigma^{\prime\prime(0)f}_{\tau}+\Sigma^{\prime\prime(2)f}_{\tau}\rangle+{q\over
m_{N}}\tilde{d}_{16}^{\tau}\langle\Sigma^{\prime\prime
g}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{|}{q\over
2m_{N}}(\tilde{d}_{2}^{\tau}-2\tilde{d}_{8}^{\tau}+{2q\over
m_{L}}\tilde{d}_{12}^{\tau})\langle\Sigma^{\prime\prime
g}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}+\Big{|}\tilde{d}_{7}^{\tau}\langle\Sigma^{\prime\prime
g}_{\tau}-\Sigma^{\prime\prime(0)f}_{\tau}\rangle+{q\over
m_{L}}\tilde{d}_{11}^{\tau}\langle\Sigma^{\prime\prime(0)f}_{\tau}+\Sigma^{\prime\prime(2)f}_{\tau}\rangle\Big{|}^{2}_{\tau\tau^{\prime}}\Big{]}\bigg{\\}}.$
(121)

These results are used Sec. IV.2.

##  Appendix C NRET Mathematica and Python Scripts

In this appendix, we describe the open-source Mathematica and Python scripts
that were employed for the calculations reported in this paper, which as noted
below are available to others. Possible uses include explorations of specific
UV theories of CLFV that can be reduced to a nucleon-level operator basis or
experimental investigations of alternative targets to explore relative LEC
sensitivities. The Mathematica script interrogates users about desired input,
and generally includes defaults as well as menus allowing parameter
investigation.

Table 11:  Input parameters and output quantities for the muon and electron Dirac solutions discussed in the text. Nucleus |  $c$ (fm) |  $\beta$ (fm) |  $\sqrt{\langle r^{2}\rangle}$ (fm) |  $E_{\mu}^{\mathrm{bind}}$ (MeV) |  $\int_{0}^{\infty}|F_{1s}|^{2}dr$ |  Z |  $Z_{\mathrm{eff}}$ |  R |  $q$ (MeV) |  $q_{\mathrm{eff}}$ (MeV)
---|---|---|---|---|---|---|---|---|---|---
12C | 2.215 | 0.491 | 2.505 | 0.1000 | 0.00047 | 6 | 5.7030 | 0.8587 | 105.07 | 108.40
16O | 2.534 | 0.514 | 2.739 | 0.1775 | 0.00083 | 8 | 7.4210 | 0.7982 | 105.11 | 109.16
19F | 2.580 | 0.567 | 2.904 | 0.2242 | 0.00104 | 9 | 8.2298 | 0.7646 | 105.12 | 109.44
23Na | 2.760 | 0.543 | 2.940 | 0.3337 | 0.00154 | 11 | 9.8547 | 0.7190 | 105.07 | 110.25
27Al | 3.070 | 0.519 | 3.062 | 0.4630 | 0.00211 | 13 | 11.3086 | 0.6583 | 104.98 | 110.81
28Si | 3.140 | 0.537 | 3.146 | 0.5346 | 0.00241 | 14 | 12.0009 | 0.6299 | 104.91 | 111.03
32S | 3.161 | 0.578 | 3.239 | 0.6924 | 0.00308 | 16 | 13.1839 | 0.5595 | 104.78 | 111.56
40Ca | 3.621 | 0.563 | 3.499 | 1.0585 | 0.00453 | 20 | 15.6916 | 0.4830 | 104.45 | 112.28
48Ti | 3.843 | 0.588 | 3.693 | 1.2615 | 0.00527 | 22 | 16.6562 | 0.4340 | 104.28 | 112.43
56Fe | 4.111 | 0.558 | 3.800 | 1.7182 | 0.00690 | 26 | 18.6028 | 0.3663 | 103.84 | 113.16
63Cu | 4.218 | 0.596 | 3.947 | 2.0884 | 0.00811 | 29 | 19.8563 | 0.3210 | 103.48 | 113.50
184W | 6.510 | 0.535 | 5.421 | 9.0851 | 0.01169 | 74 | 32.2914 | 0.0831 | 96.54 | 114.95
Table 12: Shell-model values of the six nuclear response functions $W_{O}^{\tau\tau^{\prime}}=W_{OO}^{\tau\tau^{\prime}}(q_{\mathrm{eff}})$. Response |  C |  O |  F |  Na |  Al |  Si |  S |  Ca |  Ti |  Fe |  Cu
---|---|---|---|---|---|---|---|---|---|---|---
$W^{00}_{M}$ | 72 | 1.3$\cdot 10^{2}$ | 1.5$\cdot 10^{2}$ | 2.1$\cdot 10^{2}$ | 2.6$\cdot 10^{2}$ | 2.8$\cdot 10^{2}$ | 3.3$\cdot 10^{2}$ | 4.4$\cdot 10^{2}$ | 5.2$\cdot 10^{2}$ | 5.9$\cdot 10^{2}$ | 6.5$\cdot 10^{2}$
$W^{11}_{M}$ | 4.7$\cdot 10^{-3}$ | 2.6$\cdot 10^{-3}$ | 0.28 | 0.28 | 0.28 | 4.5$\cdot 10^{-2}$ | 4.2$\cdot 10^{-2}$ | 4.2$\cdot 10^{-2}$ | 1.6 | 1.5 | 2.7
$W^{00}_{\Sigma^{\prime}}$ | 5.4$\cdot 10^{-4}$ | 0 | 0.51 | 7.9$\cdot 10^{-2}$ | 9.1$\cdot 10^{-2}$ | 1.9$\cdot 10^{-3}$ | 1.2$\cdot 10^{-4}$ | 0 | 3.7$\cdot 10^{-3}$ | 1.9$\cdot 10^{-4}$ | 0.15
$W^{11}_{\Sigma^{\prime}}$ | 3.8$\cdot 10^{-4}$ | 0 | 4.7$\cdot 10^{-2}$ | 5.7$\cdot 10^{-2}$ | 5.4$\cdot 10^{-2}$ | 1.1$\cdot 10^{-3}$ | 2.5$\cdot 10^{-5}$ | 0 | 2.5$\cdot 10^{-3}$ | 1.2$\cdot 10^{-4}$ | 0.13
$W^{00}_{\Sigma^{\prime\prime}}$ | 1.6$\cdot 10^{-3}$ | 0 | 0.25 | 6.1$\cdot 10^{-2}$ | 0.11 | 2.6$\cdot 10^{-3}$ | 5.8$\cdot 10^{-4}$ | 0 | 8.8$\cdot 10^{-3}$ | 3.9$\cdot 10^{-5}$ | 0.12
$W^{11}_{\Sigma^{\prime\prime}}$ | 1.3$\cdot 10^{-3}$ | 0 | 0.27 | 4.8$\cdot 10^{-2}$ | 9.2$\cdot 10^{-2}$ | 2.2$\cdot 10^{-3}$ | 3.2$\cdot 10^{-4}$ | 0 | 7.1$\cdot 10^{-3}$ | 7.4$\cdot 10^{-7}$ | 0.12
$W^{00}_{\Delta}$ | 1.5$\cdot 10^{-3}$ | 0 | 1.6$\cdot 10^{-4}$ | 0.21 | 0.50 | 1.2$\cdot 10^{-3}$ | 2.7$\cdot 10^{-3}$ | 0 | 5.9$\cdot 10^{-2}$ | 1.0$\cdot 10^{-3}$ | 8.0$\cdot 10^{-2}$
$W^{11}_{\Delta}$ | 5.8$\cdot 10^{-4}$ | 0 | 4.3$\cdot 10^{-2}$ | 4.6$\cdot 10^{-2}$ | 0.21 | 1.2$\cdot 10^{-3}$ | 1.2$\cdot 10^{-3}$ | 0 | 4.7$\cdot 10^{-2}$ | 1.3$\cdot 10^{-3}$ | 9.0$\cdot 10^{-3}$
$W^{00}_{\tilde{\Phi}^{\prime}}$ | 0 | 0 | 0 | 5.7$\cdot 10^{-3}$ | 8.8$\cdot 10^{-4}$ | 0 | 2.9$\cdot 10^{-4}$ | 0 | 2.8$\cdot 10^{-2}$ | 0 | 2.3$\cdot 10^{-3}$
$W^{11}_{\tilde{\Phi}^{\prime}}$ | 0 | 0 | 0 | 5.0$\cdot 10^{-4}$ | 6.6$\cdot 10^{-2}$ | 0 | 8.8$\cdot 10^{-5}$ | 0 | 1.4$\cdot 10^{-3}$ | 0 | 4.7$\cdot 10^{-2}$
$W^{00}_{\Phi^{\prime\prime}}$ | 0.90 | 1.0$\cdot 10^{-3}$ | 0.28 | 3.9 | 11 | 10 | 10 | 0.10 | 14 | 38 | 53
$W^{11}_{\Phi^{\prime\prime}}$ | 8.1$\cdot 10^{-4}$ | 1.0$\cdot 10^{-3}$ | 6.3$\cdot 10^{-2}$ | 5.6$\cdot 10^{-2}$ | 8.3$\cdot 10^{-2}$ | 5.0$\cdot 10^{-3}$ | 4.9$\cdot 10^{-2}$ | 0.10 | 3.7 | 0.56 | 1.4

Atomic and nuclear physics input: The atomic input is derived from solving the
Dirac equation in the nuclear Coulomb field, using the nuclear charge density
profile described in the text; the parameters $c$ and $\beta$ were either
taken directly from Ref. [43] or fit to the alternative density profiles
provided in that reference. The density profiles were determined from fits to
elastic electron scattering and thus are consistent with the known charge
radii of the targets. Table 11 gives for each target the values of $c$ and
$\beta$ we employed, as well as the resulting rms charge radius, the muon
binding energy, the muon lower-component probability, the computed
$Z_{\mathrm{eff}}$ that relates the muon wave function (averaged over the
nuclear density, weighted by $j_{0}(qr)$) to the point Schrödinger density at
the origin, the equivalent reduction factor $R$, the three-momentum transfer
$q$, and the computed $q_{\mathrm{eff}}$. The values for
$E_{\mu}^{\mathrm{binding}}$, $Z_{\mathrm{eff}}$, and $q_{\mathrm{eff}}$ are
embedded in the script. The atomic physics is evaluated for the principal
isotope indicated in the table, then used for all contributing isotopes.

The script includes options for 11 targets, ranging from C to Cu. The
selection includes the light- and medium-mass nuclei that were employed in
past $\mu\rightarrow e$ conversion experiments or are under consideration for
future experiments, as well as a few other common targets. As described in the
text, response function evaluation requires the one-body density matrix:
knowledge of the full one-body density matrix allows one to calculate the
exact many-body matrix element of any one-body operator. The shell model (SM)
approximation to the density matrix truncates the many-body Hilbert space to
some valence space and adopts a single-particle basis for representing the
Slater determinants in that space. We use the harmonic oscillator, a choice
that allows us to preserve translational invariance (as we use separable SM
spaces). With this choice, all multipole operator single-particle matrix
elements can be evaluated analytically. Table XII gives the resulting nuclear
response functions for the 11 targets we considered.

Typically, natural targets are used in experiments, so the script accesses
density matrices for every isotope with an abundance $\gtrsim$ 0.2%, computing
the corresponding responses by weighting each density matrix by the
appropriate abundance. There is an option in the script, however, to compute
the $\mu\rightarrow e$ conversion rate for a selected single isotope, in which
case the abundance for that isotope is set to unity.

The nuclear model dependence arises from the density matrix truncation
(inherent in the choice of the finite SM space of interacting Slater
determinants) and from the choice of the single-particle basis employed in the
determinants. The truncation generates the dependence on the choice of single-
particle basis; absent truncation, all complete bases are equivalent. The SM
spaces used are shown in Table 13 along with a list of effective interactions
that were used to generate the one-body density matrices in those spaces.
Effective interactions are typically constructed from iterated two-nucleon
potentials derived from $NN$ phase shift analyses, augmented by
phenomenological terms whose strengths are carefully tuned to reproduce
nuclear spectra and other properties in the mass range being addressed. In
most cases, density matrices are provided for more than one interaction, so
that users can explore uncertainties associated with alternative schemes for
fitting the effective interaction. The uncertainties associated with the SM
itself – the basis truncation and the choice of single particle basis – are
more difficult to assess. Shell model calculations were done with the BIGSTICK
code [70, 71].

Table 13:  Shell-model spaces employed and options available for the corresponding effective interactions. See text for descriptions of the two procedures used to determine the oscillator parameter $b$. Target | Isotopes | SM Space | $\begin{array}[]{c}\mathrm{Interaction}\\\ \mathrm{options}\end{array}$ | b (fm)
---|---|---|---|---
C | 12,13 | $1p$ | [72] | 1.67/1.70
O | 16,18 |  $2s$-$1d$ | [73, 97] | 1.73/1.83
|  | 4$\hbar\omega$ | [98] | 1.73/1.80
F | 19 |  $2s$-$1d$ | [73, 97] | 1.76/1.88
Na | 23 |  $2s$-$1d$ | as above | 1.80/1.83
Al | 27 |  $2s$-$1d$ | as above | 1.84/1.85
Si | 28-30 |  $2s$-$1d$ | as above | 1.85/1.89
S | 32-34 |  $2s$-$1d$ | as above | 1.88/1.91
Ca | 40,42,44 |  $2p$-$1f$ | [74, 99, 100] | 1.94/2.02
Ti | 46-50 |  $2p$-$1f$ | as above | 1.99/2.09
Fe | 54,56-58 |  $2p$-$1f$ | as above | 2.03/2.08
Cu | 63,65 |  $1f_{5\over 2}$-$2p$-$1g_{9\over 2}$ | [58, 101, 102] | 2.07/2.12

The script provides a default value for the harmonic oscillator size parameter
$b$, but allows the user to replace that with any desired alternative value.

User options: The script gives the user various choices:

  1. 1.

Target: The 11 options listed in Table 13 are available.

  2. 2.

Isotopic composition: The default option is a natural target, in which case
all isotopes with abundance $\gtrsim$ 0.2% are included, weighted by their
abundances, as described above. The user has the option to instead select a
single isotope.

  3. 3.

Density matrix: The user is prompted with a set of shell-model effective
interactions. Once a choice is made, the script extracts the corresponding
one-body density from a library, for use in evaluating the nuclear response
functions.

  4. 4.

Oscillator parameter: The first entries for $b$ in the last column of Table 13
are the default values obtained from the formula

$b=\sqrt{41.467\over 45\bar{A}^{-1/3}-25\bar{A}^{-2/3}}~{}\mathrm{fm}$

but the user can specify another value. (For example, the $b$ values listed
second in Table 13 reproduce the experimental rms charge radius when point
protons are assumed. This choice could be for consistency with the
distributions used to evaluate lepton Coulomb wave functions, as these were
determined from fits to measured electron scattering form factors.) Here
$\bar{A}$ is the atomic mass computed for the natural isotope, and thus
weighted over the abundances.

  5. 5.

Response functions: The user has the option to request that the calculated
nuclear response functions be displayed either as analytic functions in $y$ or
as a graph. If this option is exercised, the user is queried about each of the
six response functions and can chose to display any subset of them. The user
is also asked to specify the isospin, with options including isoscalar,
isovector, proton-only, and neutron-only couplings.

  6. 6.

Relativistic or nonrelativistic LECs: The user has the option of using the 20
relativistic operators of Table 7 that arise for scalar and vector exchanges,
rather than the basis of nonrelativistic operators. If this option is
exercised, it asks the user to specify the leptonic scale $m_{L}$ (in MeV)
that is used in Table 7.

  7. 7.

LEC input: The user either inserts the relativistic LECs
$\\{i,\tilde{d}_{i}(0),\tilde{d}_{i}(1)\\}$ or nonrelativistic LECs
$\\{i,\tilde{c}_{i}(0),\tilde{c}_{i}(1)\\}$. Here $i$ is the operator index
corresponding to the lists within this paper, and the two couplings are the
isoscalar and isovector components. Relativistic LECs are immediately
converted into their nonrelativistic equivalents, as in Table 7. The LECs are
dimensionless, given in units of $1/v^{2}$. Only those LECs of interest have
to be inserted: entry of a triad with $i=0$ terminates the input. LECs can be
inserted as numbers or as symbols.

The output is the $\mu\rightarrow e$ conversion rate in 1/sec.

Obtaining the Script and Accessing the Density matrix Library: The Mathematica
and Python scripts may be accessed via the public GitHub repository [103]. The
included README file describes the setup of the scripts. The associated
library of elastic one-body density matrices for relevant isotopes is
available from a second public repository [104]. The branch label for both
repositories corresponding to this paper is v1.0.

## References

  * [1] A. Baldini et al., arXiv:1812.06540.
  * [2] Y. Kuno and Y. Okada, Rev. Mod. Phys. 73, 151 (2001).
  * [3] W. J. Marciano, T. Mori, and JH. M. Roney, Annu. Rev. Nucl. Part. Sci. 58, 315 (2008).
  * [4] A. de Gouvêa and P. Vogel, Prog. Part. Nucl. Phys. 71, 75 (2013).
  * [5] R. Barbieri and L. J. Hall, Phys. Lett. B 338, 218 (1994).
  * [6] V. Cirigliano, R. Kitano, Y. Okada, and P. Tuzon, Phys. Rev. D 80, 013002 (2009).
  * [7] L. Calibbi and G. Signorelli, La Rivista del Nuovo Cimento 41, 71 (2018).
  * [8] R. H. Bernstein and P. S. Cooper, Phys. Rept. 532, 27 (2013).
  * [9] L. Bartoszek et al. (The Mu2e Collaboration), Fermilab Report No. Fermilab-TM-2594 and Fermilab-DESIGN-2014-1 (unpublished) arXiv:1501.05241.
  * [10] F. Abusalma et al., Fermilab report FN-1052 (unpublished), arXiv:1802.02599.
  * [11] R. Abramishvili et al. (COMET Collaboration), Prog. Theor. Exp. Phys. 2020, 033C01.
  * [12] J.-C. Angelique et al. (COMET Collaboration), arXiv:1812.07824.
  * [13] The DeeMe Collaboration, Proposal of an Experiment Searching for m-e Conversion in Nuclear Field at J-PARC Muon Facility, J-PARC proposal, 2010 (unpublished).
  * [14] A. L. Fitzpatrick, W. C. Haxton, E. Katz, N. Lubbers, and Y. Xu, J. Cosmol. Astropart. Phys. 02 (2013) 004.
  * [15] N. Anand, A. L. Fitzpatrick, and W. C. Haxton, Phys. Rev. C 89, 065501 (2014).
  * [16] Evan Rule, W. C. Haxton, and Ken McElvain, Phys. Rev. Lett 130, 131901 (2023).
  * [17] F. Bishara, J. Brod, B. Grinstein and J. Zupan, arXiv:1708.02678.
  * [18] K. Schneck et al. (SuperCDMS Collaboration), Phys. Rev. D 91, 092004 (2015).
  * [19] M. F. Albakry et al. (SuperCDMS Collaboration), arXiv:2205.11683.
  * [20] N. A. Larsen, Ph.D. thesis, Yale University, 2016 (unpublished).
  * [21] D. S. Akerib et al. (LUX Collaboration), Phys. Rev. D 103, 122005 (2021).
  * [22] J. Xia et al. (PandaX-II Collaboration), Phys. Lett. B 792, 193 (2019).
  * [23] S. Weinberg and G. Feinberg, Phys. Rev. Lett. 3, 111 and 245(E) (1959).
  * [24] W. J. Marciano and A. I. Sanda, Phys. Rev. Lett. 38, 1512 (1977).
  * [25] O. Shanker, Phys. Rev. D 20, 1608 (1979).
  * [26] T. S. Kosmas and J. D. Vergados, Nucl. Phys. A 510, 641 (1990).
  * [27] H. C. Chiang, E. Oset, T. S. Kosmas, A. Faessler, and J. D. Vergados, Nucl. Phys. A 559, 526 (1993).
  * [28] T. S. Kosmas, J. D. Vergados, O. Civitarese, and A. Faessler, Nucl. Phys. A 570, 637 (1994).
  * [29] A. Czarnecki, W. J. Marciano, and K. Melnikov, in The Workshop on Physics at the First Muon Collider and at the Front End of a Muon Collider, 6-9 November 1997, Batavia, IL., AIP Conf. Proc. No. 435 (AIP, New York, 1998), p. 409.
  * [30] T. Siiskonen, J. Suhonen, and T. S. Kosmas, Phys. Rev. C 62, 035502 (2000).
  * [31] T. S. Kosmas, Nucl. Phys. A 683, 443 (2001).
  * [32] R. Kitano, M. Koike, and Y. Okada, Phys. Rev. D 66, 096002 (2002).
  * [33] T. S. Kosmas, Nucl. Instum. Methods in Phys, Res., Sect. A 503, 247 (2003).
  * [34] A. Crivellin, S. Davidson, G. M. Pruna, and A. Signer, J. High Energy Phys. 2017, 117 (2017).
  * [35] A. Bartolotta and M. J. Ramsey-Musolf, Phys. Rev. C 98, 015208 (2018).
  * [36] V. Cirigliano, S. Davidson and Y. Kuno, Phys. Lett. B 771, 242 (2017).
  * [37] S. Davidson, Y. Kuno, and A. Saporta, Eur. Phys. J. C 78, 109 (2018).
  * [38] O. Civitarese and T. Tarutina, Phys. Rev. C 99, 065504 (2019).
  * [39] J. Heeck, R. Szafron, and Y. Uesaka, Nucl. Phys. B 980 115833 (2022).
  * [40] V. Cirigliano, K. Fuyoto, M. J. Ramsey-Musolf, and E. Rule, Phys. Rev. C 105, 055504 (2022).
  * [41] M. Hoferichter, J. Menéndez and F. Noël, Phys. Rev. Lett. 130, 1231902 (2023).
  * [42] J. D. Walecka, in Muon Physics, ed. V. W. Hughes and C. S. Wu (Academic, New York, 1975), Vol. II.
  * [43] H. De Vries, C. W. De Jager, and C. De Vries, Atomic Data and Nucl. Data Tables 36, 495 (1987).
  * [44] F. Lenz, Ph.D. thesis, University of Freiburg, Germany, 1971 (unpublished).
  * [45] J. Knoll, Nucl. Phys. A 223, 462 (1974).
  * [46] F. Lenz and R. Rosenfelder, Nucl. Phys. A 176, 513 (1971).
  * [47] A. Aste, C. von Arx, and D. Trautmann, Eur. Phys. J. A 26, 167 (2005).
  * [48] J. D. Bjorken and S. D. Drell, Relativistic Quantum Mechanics (McGraw-Hill, New York, 1964).
  * [49] T. W. Donnelly and W. C. Haxton, Atomic Data and Nucl. Data Tables 23, 103 (1979); W. C. Haxton and C. Lunardini, Computer Phys. Comm. 179, 345 (2008).
  * [50] T. DeForest and J. D. Walecka, Adv. Phys. 15, 1 (1966).
  * [51] J. Engel, Phys. Lett. B 264, 114 (1991).
  * [52] B. D. Serot, Nucl. Phys. A 308, 457 (1978).
  * [53] M. Baumgart et al., arXiv:2203.08204; also see arXiv:2210.03199.
  * [54] F. Bishara, J. Broad, B. Grinstein, and J. Zupan, J. Cosmol. Astropart. Phys. 02 (2017) 009.
  * [55] J. Brod, A. Gooties-Dreesbach, M. Tammaro, and J. Zupan, J. High Energy Phys. 10 (2018) 065.
  * [56] M. Hoferichter, P. Klos, J. Menendez, and A. Schwenk, Phys. Rev. D 99, 055031 (2019).
  * [57] T. Menzo and J. Zupan (private communication).
  * [58] A. Gniady, E. Caurier, and F. Nowacki (unpublished).
  * [59] B. A. Brown, W. Chung, and B. H. Wildenthal, Phys. Rev. Lett. 40, 1631 (1978).
  * [60] T. Park, D. Min, and M. Rho, Phys. Rep. 233, 341 (1993); A. Baroni, L. Girlanda, S. Pastore, R. Schiavilla, and M. Viviani, Phys. Rev. C 93, 015501 and 049902(E) (2016).
  * [61] W. C. Haxton, Phys. Rev. Lett 46, 698 (1981).
  * [62] J. Fan, M. Reece, and L.-T. Wang, J. Cosmol. Astropart. Phys. 11 (2010) 042.
  * [63] M. Cirelli, E. Del Nobile, and P. Panci, J. Cosmol. Astropart. Phys. 10 (2013) 019.
  * [64] G. Barello, S. Chang, and C. A. Newby, Phys. Rev. D 90, 094027 (2014).
  * [65] R. J. Hill and M. P. Solon,. Phys. Rev. D 91, 043505 (2015).
  * [66] C. Drischler, W. C. Haxton, K. McElvain, E. Mereghetti, A. Nicholson, P. Vranas, and A. Walker-Loud, Prog. Part. Nucl. Phys. 121, 103888 (2021).
  * [67] H. Hergert, Phys. Scr. 92, 023002 (2017).
  * [68] G. Gamow and E. Teller, Phys. Rev. 49, 895 (1936).
  * [69] W. C. Haxton, Evan Rule, and Ken McElvain, in preparation.
  * [70] C. W. Johnson, W. E. Ormand and P. G. Krastev, Comput. Phys. Commun. 184, 2761 (2013).
  * [71] C. W. Johnson, W. E. Ormand, K. S. McElvain, and H. Shan, arXiv:1801.08432.
  * [72] S. Cohen and D. Kurath, Nucl. Phys. 73, 1 (1965).
  * [73] B. A. Brown and W. A. Richter, Phys. Rev. C 74, 034315 (2006).
  * [74] A. Poves, J. Sánchez-Solano, E. Caurier, and F. Nowacki, Nucl. Phys. A 694,157 (2001) .
  * [75] Y. Kuno (COMET Collaboration), Prog. Theor. Phys. 2013, 022C01 (2013).
  * [76] A. Badertscher et al., Nucl. Phys. A 377, 406 (1982).
  * [77] D.A. Bryman, E. T. H. Clifford, M. J. Leitch, I. Navon, T. Numao, P. Schlatter et al., Phys. Rev. Lett. 55, 465 (1985).
  * [78] S. Ahmad, G. Azuelos, M. Blecher, D. A. Bryman, R. A. Burnham, E. T. H. Clifford et al., Phys. Rev. D 38, 2102 (1988).
  * [79] C. Dohmen, et al. (SINDRUM II Collaboration), Phys. Lett. B 317, 631 (1993).
  * [80] P. Wintz, in Proc. of the First International Symposium on Lepton and Baryon Number Violation, edited by H. V. Klapdor-Kleingrothaus and I. V. Krivosheina (IOP, Bristol, 1998), p. 534.
  * [81] D.A. Bryman, M. Blecher, K. Gotow, and R.J. Powers, Phys. Rev. Lett. 28, 1469 (1972).
  * [82] W. Burl et al. (SINDRUM II Collaboration), Eur. Phys. J. C 47, 337 (2006).
  * [83] W. Honecker, C. Dohmen, H. Haan, D. Junker, G. Otter, M. Starlinger et al. (SINDRUM II Collaboration), Phys. Rev. Lett. 76, 200 (1996).
  * [84] T. Suzuki, D. F. Measday, and J. P. Roalsvig, Phys. Rev. C 35, 2212 (1987).
  * [85] N. Teshima, PoS (NuFact2019), 082 (2020).
  * [86] E. Fermi and E. Teller, Phys. Rev. 72, 399 (1947).
  * [87] T. Mori (MEG Collaboration), Nuovo Cimento C 39, 325 (2017).
  * [88] A. M. Baldini, E. Baracchini, C. Bemporad et al., Eur. Phys. J. C 78, 380 (2018).
  * [89] U. Bellgardt et al., Nucl. Phys. B 299, 1 (1988).
  * [90] F. Wauters (on behalf of the Mu3e Collaboration), SciPost Phys. Proc. 5, 20 (2021).
  * [91] M. Raidal and A. Santamaria, Phys. Lett. B 421, 250 (1998).
  * [92] V. Cirigliano, A. Kurylov, M. J. Ramsey-Musolf, and P. Vogel, Phys. Rev. D 70, 075007 (2004).
  * [93] V. Cirigliano, A. Kurylov, M. J. Ramsey-Musolf, and P. Vogel, Phys. Rev. Lett. 93, 231802 (2004).
  * [94] J. Hisano, T. Moroi, K. Tobe, and M. Yamaguchi, Phys. Rev. D 53, 2442 (1996).
  * [95] Evan Rule, Thesis, UC Berkeley, 2022 (unpublished).
  * [96] A. R. Edmonds, Angular Momentum in Quantum Mechanics, 3rd ed. (Princeton University Press, Princeton, 1974).
  * [97] B. A. Brown and H. Wildenthal, Annu. Rev. Nucl. Part. Sci. 38, 29 (1988).
  * [98] W. C. Haxton and C. Johnson, Phys. Rev. Lett. 65, 1325 (1990).
  * [99] M. Honma, T. Otsuka, B. A. Brown, and T. Mizusaki, Phys. Rev. C 65, 061301(R) (2002); Phys. Rev. C 69, 034335 (2004).
  * [100] J. B. McGrory, B. H. Wildenthal, and E. C. Halbert,. Phys. Rev. C 2, 181 (1970).
  * [101] B. A. Brown (unpublished); B. Cheal, E. Mane, J. Billowes, M. L. Bissell, K. Blaum, B. A. Brown et al., Phys. Rev. Lett. 104, 252502 (2010).
  * [102] M. Honma, T. Otsuka, T. Mizusaki, and M. Hjorth-Jensen, Phys. Rev. C 80, 064323 (2009).
  * [103] https://github.com/Berkeley-Electroweak-Physics/Mu2e.
  * [104] https://github.com/Berkeley-Electroweak-Physics/Elastic.