# Time series analysis of friction force at self-affine mode of ice surface
softening

A. V. Khomenko  $|$ ${}^{\textsf{\footnotesize{\mbox{?},\penalty
1000\mbox{?}}}}$, D. T. Logvinenko  $|$
${}^{\textsf{\footnotesize{\mbox{?}}}}$ Corresponding author:
<EMAIL_ADDRESS>
(Received April 23, 2021, in final form July 7, 2021)

###### Abstract

The self-affine mode of ice softening during friction is investigated within
the rheological model for viscoelastic medium approximation. The different
modes of ice rubbing, determined by formation of surface liquid-like layer,
are studied. The analysis of time series of friction force is carried out,
namely Fourier analysis, construction of autocorrelation and difference
autocorrelation functions. The spectral power law is detected for modes of
crystalline ice as well as of a mixture of stable ice and metastable
softening. The self-similarity and aperiodic character of corresponding time
series of friction force are proved.

Key words: ice surface softening, viscoelastic medium, self-similarity,
fluctuations intensity, autocorrelation function, stick-slip friction

###### Abstract

Â ðàìêàõ ðåîëîãíî¿ ìîäåë äëÿ àïðîêñèìàö¿ â’ÿçêîïðóæíîãî ñåðåäîâèùà äîñëäæóòüñÿ
ñàìîïîäáíèé ðåæèì ðîçì’ÿêøåííÿ ëüîäó ïä àñ òåðòÿ. Âèâåíî ðçí ðåæèìè òåðòÿ
ëüîäó, ùî âèçíààþòüñÿ óòâîðåííÿì ïîâåðõíåâîãî ðäèíîïîäáíîãî øàðó. Ïðîâîäèòüñÿ
àíàëç àñîâèõ ðÿäâ ñèëè òåðòÿ, à ñàìå Ôóð’-àíàëç, ïîáóäîâà àâòîêîðåëÿöéíî¿ òà
ðçíèöåâî¿ àâòîêîðåëÿöéíî¿ ôóíêöé. Ñïåêòðàëüíèé ñòåïåíåâèé çàêîí âèÿâëÿòüñÿ äëÿ
ðåæèìâ êðèñòàëíîãî ëüîäó, à òàêîæ äëÿ ñóìø ñòéêîãî ëüîäó òà ìåòàñòàáëüíîãî
ðîçì’ÿêøåííÿ. Äîâåäåíî ñàìîïîäáíñòü òà àïåðèîäèíñòü âäïîâäíèõ àñîâèõ ðÿäâ ñèëè
òåðòÿ.

Ключов слова: ðîçì’ÿêøåííÿ ïîâåðõí ëüîäó, â’ÿçêîïðóæíå ñåðåäîâèùå,
ñàìîïîäáíñòü, íòåíñèâíñòü ôëóêòóàöé, àâòîêîðåëÿöéíà ôóíêöÿ, ïåðåðèâàñòå òåðòÿ

##  1 Introduction

Ice friction is an important fundamental and applied problem, involving
processes such as creep, fracture and melting [1, 2, 3, 4, 5, 6]. The solution
of this problem is directly related to the definition of models to describe
these processes, the description of the modes which arise while ice sliding,
and their properties. In particular, the motion on ice may lead to the melting
or softening of the ice surface by friction [5, 6]. It was assumed that the
surface layer of the ice that softened during friction is characterized by a
high density of defects. Such property resembles the behavior of ultrathin
films of a lubricant confined between solid surfaces [7, 8, 9]. Being two-
dimensional systems, they differ in properties from ordinary bulk liquids and
may behave abnormally. Accordingly, a nonlinear model of a viscoelastic medium
with thermal conductivity has been developed to describe the ice surface
premelting within the framework of the theory of phase transitions [10]. The
system of kinetic equations was suggested representing the mutually
coordinated behavior of shear strain $\varepsilon$ and stress $\sigma$, and
temperature $T$ in the surface ice layer during friction. Experimental data
and calculations [11, 12, 13, 14, 15, 16, 17] have shown that intermittent
(stick-slip) friction can be a stochastic regime in which the static and
kinetic friction forces change randomly in time. Therefore, the previously
obtained transition model was supplemented by taking into account Gaussian
fluctuations that play the role of external forces and lead to stick-slip
friction [18]. In the basic equations, the intensities of shear strain
$\varepsilon$ and stress $\sigma$, and temperature $T$ noises
$I_{\varepsilon}$, $I_{\sigma}$, and $I_{T}$ were introduced.

Moreover, in the study [19] it was shown that the intermittent regime is also
self-affine for most natural and practically important surfaces. A new
stochastic model has helped to find out the reasons for the appearance of this
mode and its features [20]. We used the measure units

$\displaystyle\varepsilon_{s}=\sigma_{s}/G_{\varepsilon},\quad\sigma_{s}=\left(c_{p}\eta_{\varepsilon}T_{c}/\tau_{T}\right)^{1/2},\quad
T_{c},$
$\displaystyle\left(\varepsilon_{s}/\tau_{\varepsilon}\right)^{2},\quad\quad\quad\quad\sigma_{s}^{2},\quad\quad\quad\quad\quad\quad
T_{c}^{2},$ (1.1)

for magnitudes $\varepsilon$, $\sigma$, $T$, $I_{\varepsilon}$, $I_{\sigma}$,
and $I_{T}$, respectively, where
$G_{\varepsilon}\equiv\eta_{\varepsilon}/\tau_{\varepsilon}\equiv
G(\omega)|_{\omega\to 0}$ is the relaxed value of the ice shear modulus
($\omega$ is the circular frequency of a periodic external influence),
$\eta_{\varepsilon}$ is the effective value of shear viscosity $\eta$,
$\tau_{\varepsilon}$ is the relaxation time of ice strain, $c_{p}$ is the heat
capacity, $T_{c}$ is the characteristic temperature, $\tau_{T}\equiv
l^{2}c_{p}/\kappa$ is the time of heat conductivity, $l$ is the distance into
which heat penetrates the ice, $\kappa$ is the thermal conductivity. The
nonlinear relaxation term of strain and feedbacks were considered in the
governing equations with the fractional exponent $0<a<1$ [20]:

$\displaystyle\tau_{\varepsilon}\dot{\varepsilon}=-\varepsilon^{a}+\sigma+\sqrt{I_{\varepsilon}}\xi_{1}(t),$
(1.2)
$\displaystyle\tau_{\sigma}\dot{\sigma}=-\sigma+g(T-1)\varepsilon^{a}+\sqrt{I_{\sigma}}\xi_{2}(t),$
(1.3)
$\displaystyle\tau_{T}\dot{T}=(\tau_{T}Q-T)-\sigma\varepsilon^{a}+\sqrt{I_{T}}\xi_{3}(t).$
(1.4)

Here, the stress relaxation time $\tau_{\sigma}$, the constant
$g={G_{0}}/{G_{\varepsilon}}<1$ ($G_{0}\equiv G(T=2T_{c})$ is a typical value
of modulus) and the thermostat temperature (the background sliding block
temperature) $T_{e}=\tau_{T}Q$ ($Q$ is a heat flow from the sliding block to
the surface film) are presented. The function $\xi_{i}(t)$ describes
$\delta$-correlated Gaussian source (white noise) [21] with moments:

$\langle\xi_{i}(t)\rangle=0,\quad\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2\delta_{ij}\delta(t-t^{\prime}),$
(1.5)

where $\langle\ldots\rangle$ denotes a statistical ensemble average value,
$\delta(x)$ is the Dirac delta function and $\delta_{ij}$ stands for the
Kronecker delta symbol. The relaxation pattern of a viscoelastic ice during
friction is construed by the Kelvin-Voigt equation (1.2) ($a=1$), which is
commonly recognized in the rheology of ice friction [7, 8, 22]. The Landau-
Khalatnikov-type relaxation equation for stress (1.3) [7, 8, 23] takes into
account the dependence of the shear modulus on the dimensionless temperature
$G(T)=G_{0}(T-1)$. According to our approach [10], equation (1.3) is obtained
from the theory of phase transitions rather than from rheology. On the other
hand, replacing $\varepsilon/\tau_{\sigma}$ by $\partial\varepsilon/\partial
t$ in equation (1.3) at $G(T)={\rm const}$ and $a=1$ we reduce it to the
Maxwell-type equation for a viscoelastic medium that is widely approved for
description of ice dynamics [22, 24]. Equation (1.4) is the heat conductivity
expression specifying the thermal transmission from thermostat to the ice
surface, the dissipative heating of the stress-induced viscous flow and a heat
source due to dilatation and deformation energy [10]. The classification of
different regions of the phase diagram is based on the following property.
When shear strain $\varepsilon=0$, the ice is not premelted. The case
$\varepsilon\neq 0$ meets its softening. The $\varepsilon\neq 0$ appears at
temperature $T_{e}$ above critical $T_{c0}=1+g^{-1}$, but at small
$T_{e}\leqslant T_{c0}$, the ice solidifies [10, 18, 20].

It is believed that this paper will provide an opportunity to go forward in
establishing the empirical laws of ice rubbing. We are aiming at investigating
the time series of friction force to find their frequency and correlational
characteristics.

##  2 Langevin and Fokker-Planck equations

From the system (1.2)–(1.4), the Langevin equation is acquired within the
framework of the adiabatic approximation
$\tau_{\sigma},\tau_{T}\ll\tau_{\varepsilon}$ in [20]

$\tau_{\varepsilon}\dot{\varepsilon}=f_{a}(\varepsilon)+\sqrt{I_{a}(\varepsilon)}~{}\xi(t),$
(2.1)

where the deterministic force and the effective noise intensity read

$f_{a}(\varepsilon)=-\varepsilon^{a}+g\varepsilon^{a}\left(T_{e}-1\right)d_{a}(\varepsilon),~{}d_{a}(\varepsilon){\equiv}(1+g\varepsilon^{2a})^{-1},$
(2.2) $I_{a}(\varepsilon)\equiv
I_{\varepsilon}+\left[I_{\sigma}+I_{T}(g\varepsilon^{a})^{2}\right]d_{a}^{2}(\varepsilon).$
(2.3)

Since Langevin equation is a stochastic differential equation, let us write
the corresponding Fokker-Planck equation in the Stratonovich form [20]:

$\displaystyle\frac{\partial P_{a}(\varepsilon,t)}{\partial
t}=-\frac{\partial}{\partial\varepsilon}\left[f_{a}(\varepsilon)P_{a}(\varepsilon,t)\right]+\frac{\partial}{\partial\varepsilon}\left[\sqrt{I_{a}(\varepsilon)}\frac{\partial}{\partial\varepsilon}\sqrt{I_{a}(\varepsilon)}P_{a}(\varepsilon,t)\right].$
(2.4)

The stationary ($\partial P_{a}(\varepsilon,t)/\partial t=0$) probability
density of equation (2.1) solutions is derived from equation (2.4) [25]:

$P_{a}(\varepsilon)=Z^{-1}\exp\\{-U_{a}(\varepsilon)\\},$ (2.5)

where $Z$ is a normalization constant and the effective potential reads

$U_{a}(\varepsilon)=\frac{1}{2}\ln
I_{a}(\varepsilon)-\int\limits^{\varepsilon}_{0}{f_{a}(\varepsilon^{\prime})\over
I_{a}(\varepsilon^{\prime})}\mathrm{d}\varepsilon^{\prime}.$ (2.6)

The subsequent analysis consists in a numerical solution of the Langevin
stochastic differential equation (2.1) by the Euler method [20, 7, 26, 27].
With this aim, the Stratonovich’s stochastic differential equation (2.1) is
transformed into the Ito’s one:

$\mathrm{d}\varepsilon=\left[f_{a}(\varepsilon){+}\sqrt{I_{a}(\varepsilon)}\frac{\partial}{\partial\varepsilon}\sqrt{I_{a}(\varepsilon)}\right]\mathrm{d}t{+}\sqrt{I_{a}(\varepsilon)}\mathrm{d}W(t),$
(2.7)

where $\mathrm{d}W(t)$ is the Wiener process [21]. The determination of the
discrete analog of stochastic force differential
$\mathrm{d}W(t)\equiv\sqrt{\triangle t}W_{i}$, equations (2.1) and (2.3) lead
to the iteration procedure:

$\displaystyle\varepsilon_{i+1}{=}\varepsilon_{i}{+}\left(f_{a}(\varepsilon_{i}){+}\frac{ga\varepsilon_{i}^{2a-1}[gI_{T}(1{-}g\varepsilon_{i}^{2a}){-}2I_{\sigma}]}{(1{+}g\varepsilon_{i}^{2a})^{3}}\right)\triangle
t+\sqrt{I_{a}(\varepsilon_{i})\triangle t}W_{i}.$ (2.8)

Solution of (2.8) is carried out on the time range $t\in[0,T_{m}]$ with
increment $\triangle t=T_{m}/N$ ($N$ is the number of iterations). The
stochastic force $W_{i}$ with properties [21]

$\langle W_{i}\rangle=0,\quad\langle W_{i}W_{i^{\prime}}\rangle=0,\quad\langle
W^{2}_{i}\rangle\to 2$ (2.9)

can be written using the Box-Muller model [28]:

$W_{i}=\sqrt{\kappa^{2}}\sqrt{-2\ln r_{1}}\cos(2\piup r_{2}),\quad
r_{n}\in(0,1].$ (2.10)

Here, $\kappa^{2}=2$ is the variance, $W_{i}$ is the numbers with properties
(2.9) and $r_{1}$, $r_{2}$ are pseudo-random numbers characterized by
homogenous and periodical distribution.

##  3 Analysis of friction force time series

In figure 1, the time dependencies of frictional force are depicted, which are
built on the basis of approximation $F(t)=AG_{\varepsilon}|\varepsilon(t)|$ in
accordance with [20, 18]. Therein, we have shown that $F(t)$ corresponds to
the experiments described in [11, 12, 13, 14, 15, 16, 17] for the following
quantities: 1) the strain relaxation time $\tau_{\varepsilon}\approx 0.1-5$ s;
2) the contact area $A\approx 10^{-6}-10^{-1}$ m2 that is typical of friction
of a polymer (poly(methyl methacrylate) PMMA), rubber and steel on ice [11,
16, 17] as well as for polycrystalline freshwater and saline ice over itself
[12, 13, 14, 15]; 3) the shear modulus of ice $G_{\varepsilon}\approx 0.1-10$
MPa, i.e., yield stress [11].

Figure 1: Friction force time series $F(t)$ (in arbitrary units) and its
distribution $P_{a}$ obtained by solution of (2.8) and (2.5)
($g=0.8,T_{e}=1.2,a=0.9,N=10^{6},T_{m}=50,dt=5\cdot 10^{-5}$). The panels meet
the parameters: $I_{\sigma}=4.25,I_{T}=2$ [softened ice $(S)$];
$I_{\sigma}=0,I_{T}=5$ [crystalline ice $(I)$]; $I_{\sigma}=0,I_{T}=55$
[stable ice and metastable softening $(SI+MS)$]; $I_{\sigma}=5,I_{T}=12$
[metastable ice and stable softening $(MI+SS)$].

Probability density (2.5) shown in right-hand parts of figure 1 specifies the
friction modes represented by left-hand sides of this figure. Panel, where
only one maximum of (2.5) exists at most probable $F_{0}\neq 0$ is related to
the softened ice friction ($S$). It is assumed that the total friction force
$F$ involves the static and kinetic components. Ice friction region ($\it I$)
of the phase diagram is characterized by a single maximum of probability
density at $F_{0}=0$. Two-phase region $\it SI+MS$ is defined by coexistence
of distribution maxima $P_{a}(\varepsilon)$ at zero and nonzero $F_{0}$
corresponding to the stable ice and its metastable premelting. Probability
density $P_{a}(\varepsilon)$ in the area $\it MI+SS$ has two maxima at
$F_{0}=0$ and $F_{0}\neq 0$ and construes a stick-slip mode when transitions
between metastable ice and stable softened ice are realized. It has been
revealed that the crystalline ice and $SI+MS$ friction may proceed in the
self-affine regime [20]. The existence of heterogeneous quasi-liquid layers
containing two phases was confirmed experimentally [29].

The results of spectral analysis of the time dependencies of the friction
force $F(t)$ are presented in figure 2. In particular, let us test the
evolutionary dependencies for the presence of harmonic components that can be
formed due to the influence of additive noncorrelated noise. The spectrum
determination procedure is reduced to the application of the Fourier transform
[28, 30] that is determined by

$S_{p}=\frac{1}{P}\int_{0}^{P}f(t)\exp(-\mathrm{i}\omega t)\mathrm{d}t,$ (3.1)

where $P$ is function $f(t)$ period, $\omega=2\piup/P$ is fundamental
frequency. The calculation of the coefficients $S_{p}$ is usually carried out
by using the fast Fourier transform method, since it is a recursive algorithm,
which allows us to significantly save time during each iteration calculus
[28].

This paragraph shortly describes the stages on what the algorithm of the fast
Fourier transform is based. At first, we need to convert the problem into a
discrete form. To do this, the discrete instants of time
$t_{n}=n\cdot\triangle t$ are introduced, where $\triangle t$ is the sampling
period. Further, we calculate the discrete values of the function at these
instants of time $x_{n}=f(n\cdot\triangle t)$. In this case, the total period
of the function is given by the multiplication of the total number of points
and the sampling period $P=N\cdot\triangle t$. Correspondingly, the sampling
rate of the signal equals $\omega=2\piup/P=2\piup/(N\cdot\triangle t)$. As a
result, equation (3.1) leads to the definition of a discrete Fourier transform
for a one-dimensional array $x_{n}$ of length $N$

$S_{p}=\frac{1}{N}\sum_{n=0}^{N-1}x_{n}\exp\left(\frac{-2\piup\mathrm{i}nk}{N}\right),$
(3.2)

where $k=0,\>...,\>N-1$ is the index of discrete Fourier transform over the
frequency [28].

Figure 2: Power spectrum of the signal $S_{p}(\omega)$ (in conventional units)
for the $I$ (a) and $SI+MS$ (b) friction modes, which are depicted in figure 1
at the sampling rate $\omega_{n}=1/\triangle t\>\textrm{Hz}$ and the $1024$
frequency bands (Fourier lines). Dashed lines are determined by
$S_{p}(\omega)\propto 1/\omega^{0.9}$ (a) and $S_{p}(\omega)\propto
1/\omega^{0.91}$ (b).

The time dependencies of the friction force (figure 1) are analyzed on the
basis of the algorithm of the fast Fourier transform [28], which is described
above. It calculates the spectrum of $F(t)$ evolution using (3.2). Thus, the
function (3.2) computes a vector of complex numbers. In fact, the modules of
these values represent the amplitudes of the corresponding frequencies of the
harmonic components of the signal, while the arguments give their initial
phases. The entire spectrum equals a positive part with a doubled amplitude.
Therefore, negative frequencies are neglected since a mirror image has no
meaning. After normalization, we get the power spectrum $S_{p}(\omega)$
(figure 2). In the friction force time series, there are no expressed
(selected) frequencies of their regular (periodic) components, since
$S_{p}(\omega)$ has no maxima. Let us note that $F(t)$ series are built using
iterations (2.8) at $P=100$, $\triangle t=10^{-4}$ ($I$), $P=50$, $\triangle
t=5\cdot 10^{-5}$ ($SI+MS$), $N=10^{6}$.

It is apparent from figure 2 that the power spectrum $S_{p}(\omega)$ is mainly
the same for both friction modes, i.e., the noise intensities indicated in
figure 1. Spectral densities descent monotonously in logarithmic axes with an
increase in the frequency. This implies that fluctuations $I_{a}(\varepsilon)$
are more powerful at low $\omega$. Correspondingly, the approximating dashed
lines $S_{p}(\omega)=0.21/\omega^{0.9}$ (see figure 2a) and
$S_{p}(\omega)=0.37/\omega^{0.91}$ (see figure 2b), also decrease with the
ascent of frequency. The spectrum (energy) of fluctuations of the evolutionary
variables is inversely proportional to the frequency for both ice rubbing
modes. It means that various time correlations are present. This behavior is
contrary to the properties of white noise, because in the case of white noise,
the spectral density of the random signal assumes a constant value
$S_{p}(\omega)={\rm const}$. It corresponds to the range of all possible
frequencies, i.e., the signal power at all frequencies is the same.
Consequently, there are no correlations in the system [21].

It is commonly recognized that natural systems, for which the phase
transitions are inherent, have stochastic (fluctuation) processes with
spectral density $\propto\omega^{-1}$ known as $1/\omega$, flicker or
fractional (fractal) noise [21, 31, 32]. Particularly, such a pattern takes
place in many nonequilibrium systems that are characterized by Langevin
stochastic source [33, 34, 35]. Thus, the spectrum cut-off is realized owing
to the nonlinear form of Langevin equation (2.1), describing mutually
coordinated behavior of $\varepsilon,\sigma,T$, that counteracts the increase
in frequency. Thus, white noise $\xi(t)$, which is typical of most natural
systems, is transformed into a colored (correlated) one. The power exponents
$\alpha=0.9$ and $\alpha=0.91$, which fit the graphs in figure 2, demonstrate
the presence of a “pink” noise. The noise strongly influences the pattern of a
premelted film on the surface of ice, since it is very thin ($\sim$ 10 nm).
However, similar spectral densities of the signal have been obtained for the
refinement of grain structures during severe plastic deformation of volume
specimens without spatial restrictions [30].

As it is known, “pink” noise has a power spectrum fixed by expression
$S_{p}(\omega)\propto 1/\omega^{\alpha}$, $0<\alpha<2$. These fluctuations are
widespread in nature and intermediate between “white” with a spectrum
$S_{p}(\omega)\propto 1/\omega^{0}$ and “brown” or “red” with
$S_{p}(\omega)\propto 1/\omega^{2}$ [21, 31, 33, 34, 35, 36, 37, 38, 39, 40].
Besides, “pink” noise is a process with memory similar to the “red” one, i.e.,
the prehistory of the evolution of a casual process is considered. Contrary to
that, “white” noise has no memory (correlations). We can conclude that the
fluctuational spectrum at ice friction represents $1/\omega^{\alpha}$ or
“pink” noise, which displays the correlated fluctuations. The power-law
relationships for approximations of spectral power density $S_{p}(\omega)$
confirm the self-affinity of friction force time series [19].

Further, we construct a friction force-difference autocorrelation function
(ACF) for time series analysis, using the formula [41]:

$DA_{cf}(\triangle t)=\frac{1}{2}\left\langle\left\\{F(t)-F(t+\triangle
t)\right\\}^{2}\right\rangle_{t}.$ (3.3)

Here, $F(t)$ is the friction force as a function of time, while
$\langle\ldots\rangle_{t}$ denotes a time average value taken either over a
fairly large range or over a periodically repeated regions. The coefficient
1/2 is introduced into the formula (3.3) in order to convert
$DA_{cf}(\triangle t)$ into a scaling form for the $F(t)$ variance that gives
the mean-square $F(t)$ value at large $\triangle t$. Using a difference ACF
permits to detect the self-affine properties of friction force time series
additionally to an ordinary fast Fourier transform. Figure 3 shows the power
type of the difference ACF. It is important that the graphs are constructed in
logarithmic coordinates, and this in turn allows us to trace the linear
approximation. Hence, the power-law relationship for the difference ACF, which
directly indicates $F(t)$ self-similarity. Besides, the values of power
exponent smaller than $1$ are typical of the ice surface since it may have
large random mean square slopes, i.e., fragility, and may be flattened under
the external action [19].

Figure 3: Friction force-difference ACF $DA_{cf}(\triangle t)$ (in
conventional units) for the $I$ (a) and $SI+MS$ (b) friction modes, which are
represented in figures 1 and 2. Dashed lines are defined by $DA_{cf}(\triangle
t)\propto\triangle t^{0.28}$ (a) and $DA_{cf}(\triangle t)\propto\triangle
t^{0.26}$ (b).

Let us conduct a supplemental analysis of the random fluctuations of the
rubbing force $F(t)$ by computing an ACF [41, 42]

$A_{cf}(\tau)=\kappa^{2}-DA_{cf}(\tau),$ (3.4)

where $\kappa^{2}$ is the variance. The above formula shows the relationship
between two random processes shifted relative to each other by several time
$\tau$ (a signal delay). Autocorrelation function characterizes the complex
oscillations, permitting to find their periodic components. The study of the
extreme values of the obtained correlograms allows us to define the
correlation times. It is supposed that the ACF is periodic, if the generic
function is periodic [42]. For a discrete time series with mathematical
expectation $\lambda$, autocorrelation dependence is computed with the help of
equation [42]

$A_{cf}(\tau)=\frac{1}{N\kappa^{4}}\sum_{t=1}^{N-\tau}(F_{t}-\lambda)(F_{t-\tau}-\lambda),$
(3.5)

where $N$ is the number of points in time series, $\tau$ is the lag (delay)
time.

Figure 4 demonstrates the ACFs of the basic evolutionary dependence and the
shifted one by an amount $\tau$ (shift value) in time, computed by the formula
(3.5). Figures 4a and 4b are calculated for the corresponding time series of
the friction force in figure 1. It is revealed that at the beginning, the ACFs
decrease exponentially to zero (see inset) and, then, present damped
oscillations around a certain value till approximately $5\cdot 10^{5}$ and
$8\cdot 10^{5}$ for $I$ and $SI+MS$ modes, respectively [20]. Subsequently,
ACF becomes zero at larger $\tau$, as it should be for the stationary process.
Thus, the periodical oscillations are absent. It is worth noting that in
accordance with the Cheddok scale in the range $0<A_{cf}<0.3$, the connection
strength between the time series is weak. Hence, the evolution of friction
force and autocorrelation occur by casual manner. The distribution functions
$P_{a}(\varepsilon)$ (2.5) for friction modes $\it I$ ($I_{T}=5$) and $\it
SI+MS$ ($I_{T}=55$) are described by a power law corresponding to self-affine
regime [20]. This mode is restricted by conditions $\varepsilon\ll 1$ and
$I_{\varepsilon},I_{\sigma}\ll I_{T}$. Thus, the ACFs have been obtained for
different intensities of thermal fluctuations.

Figure 4: Autocorrelation functions $A_{cf}(\tau)$ (3.5) of time series. The
correlations, depicted in figures (a) and (b), are built for the time
dependencies of the friction force at the conditions of figures 1–3 and values
$P=100$, $\triangle t=10^{-4}$ ($I$), $P=50$, $\triangle t=5\cdot 10^{-5}$
($SI+MS$), $N=10^{6}$. Curves with circles are determined by equations
$A_{cf}(\tau)=0.707\exp(-8.514\cdot 10^{-5}\tau)$ (a) and
$A_{cf}(\tau)=1.239\exp(-6.323\cdot 10^{-5}\tau)$ (b).

As it is visible from figures, in both cases ACFs have the maximum at
$\tau=0$. This implies that the correlation time equals zero as for absolutely
random process (“white noise”), which has ACF in the form of the Dirac
function. However, comparison of the insets in figures 4a and 4b with the
approximation of $\delta$-function [20, 30] shows that an autocorrelation
takes place within a certain correlation time. Thus, in the studied system
correlated (colored) noise is realized. The colored Gaussian noise $\zeta(t)$
can be obtained as the solution of equation [43]:

$\frac{\mathrm{d}\zeta(t)}{\mathrm{d}t}=-\frac{1}{\tau_{\zeta}}\zeta(t)+\frac{1}{\tau_{\zeta}}\xi(t),$
(3.6)

where $\xi(t)$ represents the Gaussian white noise with intensity $I$ and
moments

$\langle\xi(t)\rangle=0,\quad\langle\xi(t)\xi(t^{\prime})\rangle=2I\delta(t-t^{\prime}).$
(3.7)

The correlation function of the process $\zeta(t)$ can be written as

$\langle\zeta(t)\zeta(t^{\prime})\rangle=\frac{I}{\tau_{\zeta}}\exp\left(-\frac{|t-t^{\prime}|}{\tau_{\zeta}}\right),$
(3.8)

where $\tau_{\zeta}$ is autocorrelation time.

The curves with circles (approximations) in the insets of figures 4a and 4b
are determined respectively by: $A_{cf}(\tau)=0.707\exp(-8.514\cdot
10^{-5}\tau)$ and $A_{cf}(\tau)=1.239\exp(-6.323\cdot 10^{-5}\tau)$. Thus, the
evolutionary dependencies of the friction force $F(t)$ demonstrated in figure
1 have autocorrelation times $\tau\equiv 1/8.514\cdot 10^{-5}=1.175\cdot
10^{4}$ and $\tau\equiv 1/6.323\cdot 10^{-5}=1.582\cdot 10^{4}$. The given
values imply that in the course of evolution the rubbing force preferably
conserves the trend during a certain time $\tau$. Dependently on the
direction, the correlation can be either positive (direct) or negative
(reverse), for example, the increasing and decreasing parts of the dependence.
At a positive correlation, the higher values of one feature correspond to the
higher values of another one, but the lower values of one feature meet the
lower values of another one. The connection of the present situation with the
history of the process is more expressed at larger autocorrelation time
$\tau$. Physically, we can predict which value of $F(t)$ will be there in the
future. That is, there appears a possibility of predicting the ice friction
type.

In the case of evolution of the friction force according to figure 4b, it is
assumed that autocorrelation time $\tau\equiv 1.582\cdot 10^{4}$ is
“effective”, since the ice surface undergoes repeated structural and kinetic
transformations. These transitions consist in the variations between two ice
rubbing modes (phases) with fixed $F(t)$ values [4, 5, 6, 29]. The presence of
two stable stationary modes in the region $SI+MS$ means that under the same
conditions during friction, both these stationary modes can be realized
simultaneously, but in different areas of the contact with the probability
displayed in figure 1. Moreover, due to the fluctuation nature of the process,
these areas can randomly transit from one stationary state to another and vice
versa. It should be noted that the formation of two-mode friction force
distributions is actually observed in experiments [5, 6, 29]. In addition, one
of the authors proposed a theory of the formation of a two-mode rubbing force
distribution, based on a stochastic approach [18, 20]. In these models, small
and large friction were formally considered as different types of surface
structures.

It is supposed that, under certain conditions and for a determined time
$\tau$, the ice slides most preferably in a steady state. According to the
probability densities of friction force realizations (right-hand side of
figure 1), it implies that the ice surface layer has different composition.
Thus, it is possible to provide the generation of a preferred stable
structural state with a necessary rubbing at subsequent times. So, the above
analysis shows that at ice friction there is a memory of the previous ice
stable states during time $\tau$. Therefore, we can conclude that the
magnitude of autocorrelation time influences the deviation of the scope of
rubbing force from the average value. Consequently, the $F(t)$, in the formed
ice structure, can take different values under the same conditions. On the
other hand, at larger $\tau$, more time is needed to obtain the desired
outcome. Thus, our study can be useful for applications for finding the
frictional conditions to reach the sought result, that is stable ice rubbing
mode.

##  4 Summary

Using the Langevin equation, the time series of friction force are calculated
for the following modes numerically: softened ice; crystalline ice; stable ice
and metastable softening; metastable ice and stable softening. The analysis of
time series for different modes of friction was performed using the fast
Fourier transform. Oscillations were obtained with a spectral power density of
the signal that is inversely proportional to the frequency and demonstrates
the realization of $1/\omega^{\alpha}$ or “pink” noise, which shows the
presence of correlated fluctuations in the system. It is concluded that at $I$
and $SI+MS$ regimes, the time series are self-affine. The analysis of the
friction force-difference autocorrelation function indicates the power type
and, correspondingly, confirms the self-similarity of the series. The study of
the autocorrelation function of the rubbing force random oscillations made it
possible to determine the type of correlation and to identify the frequency
characteristics of the process. It is revealed that the autocorrelation
function is determined by an exponential dependence and demonstrates aperiodic
behavior. Thus, these results can be useful in technical applications, since
they allow us to predict the mode of friction during a certain time.

## Acknowledgements

This work is supported by the Ministry of Education and Science of Ukraine
(Project “Atomistic and nonlinear models of formation and friction of
nanosystems”) and visitor grant of Forschungszentrum-Jülich, Germany. A. K.
thanks Dr. Bo N. J. Persson for hospitality during his stay in
Forschungszentrum-Jülich.

## References

  * [1] Kietzig A.-M., Hatzikiriakos S. G., Englezos P., J. Appl. Phys., 2010, 107, No. 8, 081101, doi:10.1063/1.3340792.
  * [2] Heorton H. D. B. S., Feltham D. L., Tsamados M., Philos. Trans. R. Soc. London, Ser. A, 2018, 376, 20170349,
doi:10.1098/rsta.2017.0349.

  * [3] Buzzard S. C., Feltham D. L., Flocco D., J. Adv. Model. Earth Syst., 2018, 10, 262283,
doi:10.1002/2017MS001155.

  * [4] Persson B. N. J., J. Chem. Phys., 2018, 149, 234701, doi:10.1063/1.5055934.
  * [5] Lahayne O., Pichler B., Reihsner R., Eberhardsteiner J., Suh J., Kim D., Nam S., Paek H., Lorenz B., Persson B. N. J., Tribol. Lett., 2016, 62, No. 2, 1–19, doi:10.1007/s11249-016-0665-z.
  * [6] Persson B. N. J., J. Chem. Phys., 2015, 143, No. 22, 224701, doi:10.1063/1.4936299.
  * [7] Khomenko A. V., Lyashenko I. A., J. Frict. Wear, 2010, 31, No. 4, 308–316, doi:10.3103/S1068366610040100.
  * [8] Khomenko A. V., Lyashenko I. A., Phys. Lett. A, 2007, 366, No. 1–2, 165–173.
  * [9] Metlov L., Myshlyaev M., Khomenko A., Lyashenko I., Tech. Phys. Lett., 2012, 38, 972–974,
doi:10.1134/S1063785012110107.

  * [10] Khomenko A. V., Khomenko K. P., Falko V. V., Condens. Matter Phys., 2016, 19, No. 3, 33002,
doi:10.5488/CMP.19.33002.

  * [11] Blackford J. R., Skouvaklis G., Purser M., Koutsos V., Faraday Discuss., 2012, 156, 243–254.
  * [12] Kennedy F. E., Schulson E. M., Jones D. E., Philos. Mag. A, 2000, 80, No. 5, 1093–1110,
doi:10.1080/01418610008212103.

  * [13] Lishman B., Sammonds P., Feltham D., Wilchinsky A., In: Proceedings of the 20th International Conference on Port and Ocean Engineering under Arctic Conditions, POAC09–66.
  * [14] Schulson E. M., Fortt A. L., J. Geophys. Res. Solid Earth, 2012, 117, No. B12, B12204,
doi:10.1029/2012JB009219.

  * [15] Sukhorukov S., Loset S., Cold Reg. Sci. Technol., 2013, 94, 1–12, doi:10.1016/j.coldregions.2013.06.005.
  * [16] Klapproth C., Kessel T., Wiese K., Wies B., Tribol. Int., 2016, 99, 169–181, doi:10.1016/j.triboint.2015.09.012.
  * [17] Marmo B. A., Blackford J. R., Jeffree C. E., J. Glaciol., 2005, 51, No. 174, 391–398,
doi:10.3189/172756505781829304.

  * [18] Khomenko A., Khomenko M., Persson B. N. J., Khomenko K., Tribol. Lett., 2017, 65, No. 2, 71,
doi:10.1007/s11249-017-0853-5.

  * [19] Persson B. N. J., Tribol. Lett., 2014, 54, No. 1, 99–106, doi:10.1007/s11249-014-0313-4.
  * [20] Khomenko A., Tribol. Lett., 2018, 66, No. 3, 82, doi:10.1007/s11249-018-1034-x.
  * [21] Gardiner C. W., Stochastic Methods, Springer-Verlag, Berlin, 4 edn., 2009.
  * [22] Kozin V., Zhestkaja V., Pogorelova A., Chizhiumov S., Dzhabrailov M., Morozov V., Kustov A., Applied Problems of Ice Cover Dynamics, Natural Sciences Academy Publishing, Moscow, 2008, (in Russian).
  * [23] Lifshits E. M., Pitaevskii L. P., Physical Kinetics, Course of Theoretical Physics, Vol. 10, Pergamon Press, Oxford, 1 edn., 1981.
  * [24] Persson B. N. J., Sliding Friction. Physical Principles and Applications, Springer-Verlag, Berlin, Heidelberg, New York, 2 edn., 2000.
  * [25] Landau L. D., Lifshitz E. M., Statistical Physics, Course of Theoretical Physics, Vol. 5, Butterworth, London, 1999.
  * [26] Khomenko A. V., Lyashenko I. A., Borisyuk V. N., Ukr. J. Phys., 2009, 54, No. 11, 1139–1148.
  * [27] Khomenko A., Lyashenko I., Tech. Phys., 2007, 52, No. 9, 1239–1243, doi:10.1134/S1063784207090241.
  * [28] Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., Numerical Recipes: the Art of Scientific Computing, Cambridge University Press, New York, 3 edn., 2007.
  * [29] Sazaki G., Zepeda S., Nakatsubo S., Yokomine M., Furukawa Y., Proc. Natl. Acad. Sci., 2012, 109, No. 4, 1052–1055, doi:10.1073/pnas.1116685109.
  * [30] Khomenko A., Troshchenko D., Metlov L., Phys. Rev. E, 2019, 100, 022110, doi:10.1103/PhysRevE.100.022110.
  * [31] Ovchinnikov I. V., Schwartz R. N., Wang K. L., Mod. Phys. Lett. B, 2016, 30, No. 08, 1650086,
doi:10.1142/S021798491650086X.

  * [32] Khomenko A. V., Prodanov N. V., Khomenko K. P., Troshchenko D. S., J. Nano- Electron. Phys., 2014, 6, No. 1, 01012.
  * [33] Lebyodkin M. A., Shashkov I. V., Lebedkina T. A., Gornakov V. S., Phys. Rev. E, 2017, 95, 032910, doi:10.1103/PhysRevE.95.032910.
  * [34] Mudrock R. N., Lebyodkin M. A., Kurath P., Beaudoin A. J., Lebedkina T. A., Scr. Mater., 2011, 65, No. 12, 1093–1096, doi:10.1016/j.scriptamat.2011.09.025.
  * [35] Khomenko A. V., Zakharov M. V., Khomenko K. P., Khyzhnya Y. V., Trofimenko P. E., In: Proceedings of the IEEE 8th International Conference on Nanomaterials: Application Properties (NAP’18) (Zatoka, Ukraine, 2018), Vol. 4, IEEE, USA, 2019, 04NNLS15–1–4, doi:10.1109/NAP.2018.8915102.
  * [36] Pogrebnjak A. D., Bagdasaryan A. A., Yakushchenko I. V., Beresnev V. M., Russ. Chem. Rev., 2014, 83, No. 11, 1027–1061, doi:10.1070/RCR4407.
  * [37] Goncharov A. A., Yunda A. N., Komsta H., Rogalski P., Acta Phys. Pol. A, 2017, 132, No. 2, 270–273, doi:10.12693/APhysPolA.132.270.
  * [38] Yushchenko O. V., Badalyan A. Y., J. Nano-Electron. Phys., 2012, 4, No. 3, 03009.
  * [39] Mazur O., Stefanovich L., Yurchenko V. M., Phys. Solid State, 2016, 58, No. 8, 1596–1604,
doi:10.1134/S1063783416080205.

  * [40] Khomenko A. V., Boyko D. V., Zakharov M. V., Khomenko K. P., Khyzhnya Y. V., In: Proceedings of the IEEE 7th International Conference on Nanomaterials: Application Properties (NAP’17) (Odessa, 2017), Vol. 6, IEEE, USA, 2017, 01NNPT01–1–4, doi:10.1109/NAP.2017.8190181.
  * [41] Wang A., Muser M. H., Tribol. Int., 2018, 123, 224–233, doi:10.1016/j.triboint.2018.02.002.
  * [42] Box G. E. P., Jenkins G. M., Reinsel G. C., Ljung G. M., Time Series Analysis: Forecasting and Control, 5 edn., Wiley Series in Probability and Statistics, Wiley, Hoboken, 2015.
  * [43] Risken H., The Fokker-Planck-Equation. Methods of Solution and Applications, Springer-Verlag, Berlin etc., 2 edn., 1989.

Ukrainian 0

Àíàëç àñîâèõ ðÿäâ ñèëè òåðòÿ ïðè ñàìîïîäáíîìó ðåæèì ðîçì’ÿêøåííÿ ïîâåðõí ëüîäó
Î. Â. Õîìåíêî${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Ä.
Ò. Ëîãâèíåíêî${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ Ñóìñüêèé äåðæàâíèé óíâåðñèòåò, âóë. Ðèìñüêîãî-
Êîðñàêîâà, 2, 40007 Ñóìè, Óêðà¿íà ${}^{\textsf{\footnotesize 2}}$ íñòèòóò
Ïåòåðà Ãðþíáåðãà-1, Äîñëäíèöüêèé öåíòð Þëõà, 52425, Þëõ, Íìåèíà


# Muon spin force

Yohei Ema<EMAIL_ADDRESS>William I. Fine Theoretical Physics Institute,
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA  School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Ting Gao<EMAIL_ADDRESS>School of Physics and
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  Maxim Pospelov
<EMAIL_ADDRESS>William I. Fine Theoretical Physics Institute, School of
Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA

###### Abstract

Current discrepancy between the measurement and the prediction of the muon
anomalous magnetic moment can be resolved in the presence of a long-range
force created by ordinary atoms acting on the muon spin via axial-vector
and/or pseudoscalar coupling, and requiring a tiny,
$\mathcal{O}(10^{-13}\,{\rm eV})$ spin energy splitting between muon state
polarized in the vertical direction. We suggest that an extension of the muon
spin resonance ($\mu$SR) experiments can provide a definitive test of this
class of models. We also derive indirect constraints on the strength of the
muon spin force, by considering the muon-loop-induced interactions between
nuclear spin and external directions. The limits on the muon spin force
extracted from the comparison of 199Hg/201Hg and 129Xe/131Xe spin precession
are strong for the pseudoscalar coupling, but are significantly relaxed for
the axial-vector one. These limits suffer from significant model
uncertainties, poorly known proton/neutron spin content of these nuclei, and
therefore do not exclude the possibility of a muon spin force relevant for the
muon $g-2$.

††preprint: UMN-TH-4219/23, FTPI-MINN-23-12

Introduction — The next few years will bring an unprecedented increase in
intensity of muon sources Aiba:2021bxe ; CGroup:2022tli ; Kawamura:2018apy ;
Thomason:2019fwe . There are many practical applications of muons such as muon
spin resonance ($\mu$SR), muon beams as calibrated sources of neutrinos, and
in the far future, the muon collider. At the same time, the studies of muon
properties, such as lepton-flavour violation and muon spin dynamics ($g-2$,
muon EDM $d_{\mu}$, muonium spectroscopy), will be brought to a new level of
accuracy. Currently, the measurement of the muon $g-2$ has put in question the
consistency with the Standard Model (SM) theoretical predictions
Aoyama:2020ynm ; Borsanyi:2020mff ; Muong-2:2021ojo . The experimental result
Muong-2:2021ojo  for the anomalous magnetic moment is larger than data-driven
value predicted in the SM Aoyama:2020ynm  by about $(3\,\mathchar
45\relax\mathchar 45\relax\,4)\sigma$, or

$\Delta a_{(\mu)}\simeq+2.5\times 10^{-9}.$ (1)

It is tempting to assign this discrepancy to effects of new physics. However,
one has to be cognisant of the fact that recent lattice QCD results
Borsanyi:2020mff ; Ce:2022kxy ; ExtendedTwistedMass:2022jpw ;
FermilabLatticeHPQCD:2023jof , as well as some re-measurement of the di-pion
production on $\gamma^{*}\to 2\pi$ processes CMD-3:2023alj  point to a larger
predicted $a_{\mu}$, that is perhaps consistent with the experimental value.
More analyses as well as independent determinations of the hadronic vacuum
polarization are needed to clarify this issue. Nevertheless, it is important
to investigate various scenarios under which the muon anomalous magnetic
moment can receive a sizeable beyond-SM correction (see e.g. Czarnecki:2001pv
).

In this paper, we will examine an alternative suggestion of a muon spin force
mediated by a very light force carrier, acting on muon spin at the tree level
and affecting its spin dynamics Janish:2020knz ; Agrawal:2022wjm ;
Davoudiasl:2022gdg . Such models almost inevitably operate with sub-eV fields,
very weakly coupled to the ordinary atoms and muons. Additional force acting
on the muon spin may alter its precession frequency, and therefore can be
interpreted as a contribution to $\Delta a_{\mu}$. Despite smallness of
couplings, the long range of a force can lead to a coherent enhancement due to
e.g. all atoms in the Earth enhancing the gradient of the scalar field at
Earth’s surface where muon spin experiments are performed. The most useful
quantity, capable of explaining discrepancy (1) is the muon-force-associated
energy splitting between up and down projections of the muon spin,

$\Delta E_{(\mu)}=6\times 10^{-14}\,{\rm eV}$ (2)

While significant parts of the available parameter space can be
constrained/excluded using the combination of the 5th force constraints and
astrophysical limits, it is important to investigate indirect constraints
imposed on muon spin force by spin dynamics of ordinary atoms. This has not
been performed before, and our paper fills this gap. The studies of spins
coupled to external directions are well known (see e.g. Pospelov:2004fj  and
references therein). One class of searches tests the coupling of atomic spins
to directions created by the external masses, that can be induced by e.g.
$CP$-odd couplings of axions Moody:1984ba . An important subset of such
experiments searches for the exotic couplings of spin to the vertical
direction, approximately parallel to the vector of the Earth’s gravitational
force, a setup that is most relevant for the muon $g-2$. These experiments can
be used, indirectly, to set bounds on exotic interactions of the muon spin. In
some sense, the problem is conceptually similar to the problem of muon EDM
inducing atomic EDMs, that was shown recently Ema:2021jds ; Ema:2022wxd  to
place more stringent bounds on $d_{\mu}$ compared to direct measurements
Semertzidis:1999kv ; Muong-2:2008ebm . In this paper, we investigate the
constraints imposed by the atomic experiments on exotic muon interactions. We
find that the constraints are generally stronger by about an order of
magnitude than the level needed to explain the discrepancy (1). However, we
also note that the calculations are beset by significant nuclear uncertainties
and near-accidental cancellations, that prevent us in deriving stronger
limits.

The second goal of this work is to point out that new dedicated experiments
with muons, similar to $\mu$SR in low/medium magnetic field, can be used to
detect/constrain the anomalous muon spin rotation. We show that already
existing muon facilities can be used to test values of the spin splitting (2)
relevant for the $g-2$ discrepancy explanation.

Muon spin force and $g-2$ — We assume the existence of the muon spin force
mediated by a very light scalar field:

$\displaystyle{\cal L}_{\mathrm{eff}}=...$
$\displaystyle+\frac{1}{2}(\partial_{\nu}\phi)^{2}-\frac{1}{2}m_{\phi}^{2}\phi^{2}-g_{S}\phi(\bar{n}n+\bar{p}p)$
$\displaystyle+\frac{g^{(\mu)}_{A}}{2m_{\mu}}\partial^{\alpha}\phi\times\bar{\mu}\gamma_{\alpha}\gamma_{5}\mu-{g_{P}^{(\mu)}}\phi\times\bar{\mu}i\gamma_{5}\mu.$
(3)

Notice that as far as the muon spin dynamics is concerned, and in the lowest
order in $\phi$, both axial $g^{(\mu)}_{A}$ and pseudoscalar $g^{(\mu)}_{P}$
couplings give equivalent effects. However, they are not equivalent at the
loop level, as they would predict e.g. different coefficients for the
$F_{\alpha\beta}\tilde{F}^{\alpha\beta}$ coupling below the muon mass, leading
eventually to different strength of the coupling to the atomic spin
Pospelov:2008jk ; Flambaum:2009mz . In particular, $g^{(\mu)}_{A}$ does not
generate couplings to on-shell photons at low energy, and therefore is immune
to searches of photon birefringence along the lines of recent proposal
Fedderke:2023dwj .

If we take one nucleon and one muon separated by distance $r$, this will
result in a potential acting on spin,

$H^{(\mu)}_{s\,\mathchar 45\relax\mathchar
45\relax\,m}=\frac{g_{S}(g^{(\mu)}_{A}+g^{(\mu)}_{P})}{4\pi\times
2m_{\mu}}(\bm{\sigma}^{(\mu)}\cdot\bm{\nabla})\frac{\exp(-m_{\phi}r)}{r}$ (4)

Here “$s\,\mathchar 45\relax\mathchar 45\relax\,m$” stands for “spin-mass”,
and superscript $(\mu)$ refers to the type of particles. $g_{S}$ is limited by
gravitational fifth-force searches; $g_{A}$ is limited by stellar energy loss
(Agrawal:2022wjm  and references therein).

First, we would like to establish a tentative size of the $g_{S}g^{(\mu)}_{P}$
suggested by the results of the muon $g-2$ experiment, if one interprets the
current discrepancy (1) as a real effect. Suppose that the deviation in the
muon $g-2$ experiment is caused by the muon spin interacting with the mostly
vertical gradient created by $\phi$ which should be the case for all ranges of
$\phi$ force, $m_{\phi}^{-1}$, exceeding several meters. We are going to
combine all constants within one quantity that has dimension of energy, and
introduce a vertical unit vector $\mathbf{n}=(0,0,1)$ aligned with
$\bm{\nabla}\phi$. Then the additional term in the Hamiltonian for the muon
spin propagating in the $g-2$ ring can be written as

$H=\Delta E_{(\mu)}(\mathbf{s}\cdot\mathbf{n})=\Delta E_{(\mu)}\,s_{z},$ (5)

with $s_{z}=\sigma_{z}/2$. Using Heisenberg equations of motion, we get in the
rest frame of muons,
$\frac{ds_{i}}{d\tau}=i[H,s_{i}]=\epsilon_{ijk}n_{j}s_{k}\times\Delta
E_{(\mu)}.$ Here $\tau$ is the proper time. This implies that spin precesses
around the vertical direction with the cyclic frequency $\Delta\omega=\Delta
E_{(\mu)}/\hbar$. The equation needs to be transformed to the lab frame, using
$\gamma\simeq 29.3$,111 This magic momentum is used at the BNL/FNAL
experiments. The J-PARC experiment plans to use a smaller momentum Abe:2019thb
, which results in a larger effect of the muon spin force Davoudiasl:2022gdg .
$dt=\gamma d\tau$, reducing the precession frequency by a gamma factor,
$\gamma^{-1}\Delta E^{(\mu)}/\hbar$ (notice that the axial-vector
$\bm{\nabla}\phi$ and the tensor $\mathbf{B}$, $\mathbf{E}$ transform
differently under the boost). The size of the effect that would account for
the observed $\Delta a_{(\mu)}$ discrepancy (1) can be determined from the
following equation:

$\Delta E_{(\mu)}=(\hbar\omega_{a})\times\gamma\times(\Delta
a_{(\mu)}/a_{(\mu)}),$ (6)

where $\omega_{a}$ is the main precession frequency measured by the $g-2$
experiment and $a_{(\mu)}\simeq\alpha/(2\pi)$. Substituting relevant
quantities we arrive to the required value of the spin coupling given in Eq.
(2).

Indirect constraints — The two main experiments that set indirect constraints
on $\Delta E_{(\mu)}$ compare spin precession around the vertical directions
of mercury isotopes, 199Hg and 201Hg PhysRevLett.68.135 , and xenon isotopes,
129Xe and 131Xe Zhang:2023qmu  that all have closed electron shells and
nonzero nuclear magnetic moments. Defining the interaction of $i$’th isotope
via nuclear spin operator $\mathbf{I}$, nuclear magneton $\mu_{N}$ and the
$g$-factor $g_{i}$ as

$H_{i}=-\mu_{N}g_{i}({\bf I\cdot B})+\Delta E_{i}(\mathbf{I}\cdot\mathbf{n}),$
(7)

one can formulate the constraints on exotic spin interactions as follows:

$\displaystyle|\Delta E_{\rm Hg}|$ $\displaystyle=\left|\Delta E_{201}-\Delta
E_{199}\,\rho_{\mathrm{Hg}}\right|<3.0\times 10^{-21}\,{\rm eV},$ (8)
$\displaystyle|\Delta E_{\rm Xe}|$ $\displaystyle=\left|\Delta E_{129}-\Delta
E_{131}\,\rho_{\mathrm{Xe}}\right|<1.7\times 10^{-22}\,{\rm eV},$ (9)

with $\rho_{\mathrm{Hg}}=g_{201}/g_{199}=-0.369139$ and
$\rho_{\mathrm{Xe}}=g_{129}/g_{131}=-3.37337$ the gyromagnetic ratios. Note
that the sensitivity to exotic spin force disappears if the $\Delta E$ for two
atomic species scale exactly as $g$ factors. More generally, the status of
atomic magnetometry relevant for these types of experiments was reviewed
recently in Terrano:2021zyh . Other notable limits (sensitive to electron and
proton couplings to Earth gravity/vertical direction) are reported in
JacksonKimball:2017zky . Note that the electron-spin-oriented torsion-pendulum
experiment Heckel:2008hw  also sets the constraint on the electron spin
coupling to Earth’s gravity, which corresponds to $\Delta E_{(e)}<6\times
10^{-18}$ eV in our notations.

Figure 1: The Feynman diagrams corresponding to one- and two-loop effective
Lagrangians. The thick line corresponds to the muon, and the cross dot
indicates the axial-vector/pseudoscalar coupling insertion. The
electromagnetic fields of the one-loop effective Lagrangian are sourced by the
nucleus.

At the next step we construct one- and two-loop effective Lagrangians,
integrating out the muon loop, whose diagrams are shown in Fig. 1. We treat
parameters $g_{A}^{(\mu)}$ and $g_{P}^{(\mu)}$ in (3) as an input, and
calculate the induced couplings to the external electromagnetic and $\phi$
fields at low energy,

$\displaystyle{\cal L}_{eff}$
$\displaystyle=-\frac{e^{2}}{16\pi^{2}m_{\mu}}\left(g_{P}^{(\mu)}\phi
F_{\mu\nu}\tilde{F}_{\mu\nu}+\frac{g_{A}^{(\mu)}\partial_{\nu}\phi}{3m_{\mu}^{2}}\tilde{F}_{\alpha\nu}\partial_{\mu}F_{\alpha\mu}\right)$
$\displaystyle=\frac{e^{2}\bm{\nabla}\phi}{4\pi^{2}m_{\mu}}\cdot\left(g_{P}^{(\mu)}\mathbf{B}A_{0}+\frac{g_{A}^{(\mu)}}{12m_{\mu}^{2}}\mathbf{B}\,\nabla^{2}A_{0}\right).$
(10)

In this expression we retained the leading order terms in $1/m_{\mu}$
expansion for each coupling, and neglected higher derivative terms acting on
$\phi$, $\nabla_{i}\nabla_{j}\phi\to 0$. The $g_{A}^{(\mu)}$ coupling does not
generate corrections to the propagation of photons. In the second line of (10)
we used time independence of all background fields, and introduced
electrostatic potential $A_{0}$ via $\mathbf{E}=-\bm{\nabla}A_{0}$.

As discussed in Flambaum:2009mz , in the presence of interaction (10), the
internal nuclear electromagnetic fields will lead to the coupling of the
nuclear spin to $\bm{\nabla}\phi$. In order to calculate observables $\Delta
E_{\rm Hg,Xe}$ in (8) and (9) induced by (10), one needs detailed knowledge of
the EM field distribution inside these nuclei, which is relatively simple for
$A_{0}(\mathbf{r})$ as it is dominated by collective effects, $A_{0}\propto
Z$, and quite difficult for $\mathbf{B}(\mathbf{r})$.

In our evaluation, we employ the following phenomenological approach. We take
the magnetic field to be the sum of the one created by the valence nucleon
(which happens to be a neutron for all these nuclei) and core polarization
that is needed to describe observed magnetic moments for these nuclei. In
other words,

$\mathbf{B}(\mathbf{r})=\mathbf{B}_{n}(\mathbf{r})+\mathbf{B}_{\rm
core}(\mathbf{r})\xrightarrow{r\gg
R_{N}}\frac{1}{4\pi}\bm{\nabla}\times(\bm{\mu}\times\bm{\nabla})\frac{1}{r},$
(11)

where $\bm{\mu}$ is the observed magnetic moment of a given nucleus. Boundary
condition (11) does not allow full reconstruction of the magnetic field
profiles. In order to achieve that, we solve for the valence neutron wave
functions numerically, using one particle Schrödinger equation in a collective
nuclear potential, and treating neutrons in 129,131Xe and 199,201Hg nuclei as
$s_{1/2}$, $d_{3/2}$, $p_{1/2}$, $p_{3/2}$ orbitals, as described in
Dmitriev:2003sc . This determines magnetic spin density, the spin current, and
ultimately $\mathbf{B}_{n}(\mathbf{r})$ everywhere, including inside the
nucleus. For $\mathbf{B}_{\rm core}$, we make a phenomenological model that it
corresponds to a sphere with constant magnetic polarization, following an
empirical fact of near constant nuclear density inside the nuclear radius
$R_{N}$, and fix its normalization by requiring that the sum (11) reproduces
the observed nuclear magnetic moments ($\mathbf{B}_{n}$ alone does not quite
reproduce the observed values except for 199Hg).

Following the steps outlined above we calculate $\Delta E_{i}$’s entering (8)
and (9). The details of the calculation are provided in the Appendix.
Parametrically, we have the following size of the nuclear spin coupling to the
vertical direction,

$\displaystyle\Delta E_{i}$ $\displaystyle=\Delta E_{(\mu)}\times
b_{i}\frac{Z\alpha^{2}(\mu_{n}/e)}{4\pi R_{Ni}}~{}~{}{\rm
for}~{}g_{P}^{(\mu)},$ (12) $\displaystyle\Delta E_{i}$ $\displaystyle=\Delta
E_{(\mu)}\times c_{i}\frac{Z\alpha^{2}(\mu_{n}/e)}{4\pi
R_{Ni}^{3}m_{\mu}^{2}}~{}~{}{\rm for}~{}g_{A}^{(\mu)},$ (13)

where $\mu_{n}$ is the neutron magnetic moment, and $R_{Ni}$ are corresponding
nuclear radii (related to the charge radii $r_{c}$ as
$R_{N}=\sqrt{5/3}\,r_{c}$). Constants $b_{i},c_{i}$ encode numerical
coefficients that arise from the numerical evaluation of the matrix elements
needed to convert (10) to $\Delta E_{i}$. Substituting numerical results, we
get the following prediction for the observables:

$\displaystyle\left\lvert\frac{\Delta E_{\rm Hg}}{\Delta
E_{(\mu)}}\right\rvert$ $\displaystyle=1\times
10^{-6};~{}\left\lvert\frac{\Delta E_{\rm Xe}}{\Delta
E_{(\mu)}}\right\rvert=3\times 10^{-6}~{}{\rm for}~{}g_{P}^{(\mu)},$ (14)
$\displaystyle\left\lvert\frac{\Delta E_{\rm Hg}}{\Delta
E_{(\mu)}}\right\rvert$ $\displaystyle=2\times
10^{-8};~{}\left\lvert\frac{\Delta E_{\rm Xe}}{\Delta
E_{(\mu)}}\right\rvert=4\times 10^{-9}~{}{\rm for}~{}g_{A}^{(\mu)}.$ (15)

Coupling $g_{P}^{(\mu)}$ is constrained stronger than $g_{A}^{(\mu)}$, and if
pseudoscalar coupling dominates the muon spin force, the muon $g-2$ anomaly
speculation is ruled out by the indirect effect, as the constraint from Xe
implies $\left\lvert\Delta E_{(\mu)}\right\rvert<6\times
10^{-17}\,\mathrm{eV}$, which is three orders of magnitude stronger than (2).
This is not the case for $g_{A}^{(\mu)}$ coupling. When interpreted as an
upper bound, the results for $g_{A}^{(\mu)}$ formally imply $|\Delta
E_{(\mu)}|<4\times 10^{-14}$ eV from Xe which is somewhat stronger than the
re-interpretation of the muon $g-2$ discrepancy (2). Notice, however, that the
results are extremely sensitive to the assumptions about the nuclear magnetic
fields, so that the coefficients in (15) are much smaller than the parametric
estimate in (13) and result from strong $\mathcal{O}(1\%)$ numerical
cancellations. Therefore, the re-interpretation of (8) and (9) constraints
into a firm limit on $\Delta E_{(\mu)}$ is problematic.

Now we turn to the calculation of the two-loop induced spin couplings. The
two-loop effects induced by the $g_{P}^{(\mu)}$ coupling is UV-finite (if the
Higgs boson contributions are neglected). The calculation of the induced
electron coupling is well defined (see e.g. Flambaum:2019ejc ; Ema:2021jds ),

$g_{P}^{(e)}=g_{P}^{(\mu)}\times\frac{3\alpha^{2}m_{e}}{2\pi^{2}m_{\mu}}\times\log(m_{\mu}/m_{e}),$
(16)

but is small and does not lead to competitive bounds. The
$g_{P}^{(\mu)}$-induced coupling to nucleons is also suppressed. The
calculation cannot be done with free quarks, and has to be performed at a
nucleon level by combining virtual E1 and M1 transitions. The most important
nucleon, on account of PhysRevLett.68.135 ; Zhang:2023qmu , is the neutron
that has a very suppressed coupling to the electric field. Parametrically, the
two-loop effect will then be suppressed by
$\alpha^{2}r_{n}^{2}m_{\rho}^{2}m_{\mu}(\mu_{n}/e)$, where
$r_{n}^{2}m_{\rho}^{2}$ is the product of the neutron charge radius and
$\rho$-mass i.e. the approximate hadronic scale of convergence for
electromagnetic loops. This combination is small, and subdominant to (12).
Therefore, we conclude that the main effect of $g_{P}^{(\mu)}$ is given by one
loop, in Eq. (14).

In contrast, the two-loop effects induced by $g_{A}^{(\mu)}$ can be important
for two reasons. It is well known that $g_{A}$ couplings will give a
logarithmic divergence at two loops, requiring the introduction of an
ultraviolet cutoff (i.e. the scale of UV completion) $\Lambda_{\rm UV}$ which
can partially compensate for loop smallness. The second reason is that the
one-loop effect in (13) are suppressed by
$(R_{N}^{3}m_{\mu}^{2}m_{p})^{-1}\sim 3\times 10^{-3}$ while two-loop effects
are not. The two-loop result at the elementary particle level is given by

$\displaystyle{\cal L}_{\mathrm{eff}}$
$\displaystyle=\sum_{i=u,d,s,e}\frac{g_{A}^{(i)}}{2m_{\mu}}\partial_{\alpha}\phi\times\bar{\psi}_{i}\gamma^{\alpha}\gamma_{5}\psi,$
$\displaystyle g_{A}^{(i)}$
$\displaystyle=-g_{A}^{(\mu)}\times\frac{3}{4}\left(\frac{\alpha}{\pi}\right)^{2}Q_{i}^{2}\log(\Lambda_{\rm
UV}^{2}/\Lambda_{\rm IR}^{2}),$ (17)

where $Q_{e}=-1,\,Q_{u}=2/3$ etc. $\Lambda_{\rm IR}$ should be identified with
the hadronic scale for quarks, $\Lambda_{\rm IR}\simeq m_{p}$, while for the
electrons the integration can be extended further down to $\Lambda_{\rm
IR}\simeq m_{\mu}$. UV cutoff depends on details of UV completion, but should
not parametrically exceed $m_{\mu}/g_{A}^{(\mu)}$. It is another fundamental
uncertainty of the low-energy model (3).

Further model dependence is exacerbated by the fact that other couplings and
other particles can contribute to $g_{A}^{(i)}$ at the same loop level. For
example, above the EW scale, $W^{+}W^{-}$ and $\gamma Z$ loops will contribute
to the running of the couplings, and one needs a proper $SU(2)\times U(1)$
completion of (3) to extend the treatment above $m_{W}$. Thus, for example, a
low energy vector coupling of muons,
$g^{(\mu)}_{V}/(2m_{\mu})\partial^{\alpha}\phi\times\bar{\mu}\gamma_{\alpha}\mu$
will contribute into $g_{A}^{(i)}$ due to running above the electroweak scale.

At the same time, it is easy to see that $\Lambda_{\rm UV}$ can be much lower
than the electroweak scale. For example $\mu\,\mathchar 45\relax\mathchar
45\relax\,\tau$ symmetry Karshenboim:2014tka  can be employed to have a full
cancellation of muon and tau loops above $\Lambda_{\rm UV}=m_{\tau}$. Given
that the couplings of $\phi$ to tau leptons are unconstrained, this choice is
perhaps the most conservative. In that case, the renormalization group
logarithm, as far as the quark couplings are considered, is given by
$\log(m_{\tau}^{2}/m_{p}^{2})\simeq 1$.

Now we would like to estimate the sensitivity of Xe and Hg experiments to
$g_{A}^{(i)}$ and ultimately to $g_{A}^{(\mu)}$. Using the neutron beta decay
measurements and lattice calculations FlavourLatticeAveragingGroupFLAG:2021npn
, one can translate quark axial couplings to those of the nucleons:

$\displaystyle
g_{A}^{(p)}=0.777g_{A}^{(u)}-0.438g_{A}^{(d)}-0.053g_{A}^{(s)};$
$\displaystyle
g_{A}^{(n)}=-0.438g_{A}^{(u)}+0.777g_{A}^{(d)}-0.053g_{A}^{(s)},$ (18)

where we have used (2+1+1)-light fermion results
FlavourLatticeAveragingGroupFLAG:2021npn .

At the next step, we need to combine neutron and proton spins into the total
spins of Xe and Hg isotopes. The isotopes of Xe and Hg all have unpaired
neutrons over the closed shells, but the accuracy of such a model in
predicting observable $\mu_{i}$ is rather poor for some isotopes (for example,
for 201Hg). Moreover, the method where a difference between observed magnetic
moment and the shell model prediction (Schmidt model) of $\mu_{i}$ is
attributed to the proton spin magnetism (see e.g. Flambaum:2009mz ) seems to
over-predict the nuclear proton spin content JacksonKimball:2014vsz .

In what follows, we consider the Xenon results (9) using the nuclear
calculations of the spin content from Klos:2013rwa , as well as two other
nuclear models quoted in the same paper. Combining the dominant neutron, and
subdominant proton contributions using additional inputs from (17) and (18),
we arrive at

$\Delta E_{\rm Xe}=(3\,\mathchar 45\relax\mathchar 45\relax\,7)\times
10^{-8}\times\Delta E_{(\mu)}\times\frac{\log(\Lambda_{\rm
UV}^{2}/\Lambda_{\rm IR}^{2})}{\log(m_{\tau}^{2}/m_{p}^{2})}.$ (19)

The uncertainty in this coefficient comes mostly from the variation of the
nucleon spin content of Xe isotopes, from model to model. Still, it appears to
be dominant over one loop result of (15) by a factor of a few, implying a
limit on muon spin force, $|\Delta E_{(\mu)}|<6\times 10^{-15}$ eV.

We conclude that the results of the Xe and Hg experiments set indirect
constraints on $\Delta E_{(\mu)}$ about one order of magnitude stronger than
the value needed to reconcile the predicted and measured values of the muon
$g-2$, Eq. (2). However, we note that there are sources of multiple
cancellations at the hadronic and nuclear level, as well as strong sensitivity
to UV completions of such models that negates the utility of such bounds. For
example, a $\mathcal{O}(10\%)$ accidental cancellation of different sources
inside $g_{A}^{(n)}$ would completely negate the utility of these bounds to
challenge the $g-2$ suggested value. Therefore, one should check for the
presence/absence of $\Delta E_{(\mu)}$ directly, using the high-intensity muon
sources.

Anomalous spin precession of stopped muons — The effect of an exotic spin
force (5) induces spin-precession of the stopped $\mu^{+}$ particles. If
magnetic field is completely shielded, the muonic spin force will introduce
the rotation angle during one muon lifetime

$\displaystyle\Delta\psi=\frac{\Delta
E_{(\mu)}t}{\hbar}~{}\Longrightarrow~{}\Delta\psi(t=\tau_{\mu})=2\times
10^{-4},$ (20)

where $\tau_{\mu}=2.2\,{\rm\mu s}$ is muon lifetime. The minimum number of
stopped muons required to detect such an angle, assuming large efficiency of
the positron detection, is given by $N_{(\mu)}>1/(\Delta\psi)^{2}\sim 10^{8}$.
With modern $\mu$SR beam lines capable of delivering over $10^{5}$ muons per
second, this appears to be entirely feasible.

On the other hand, the equivalent magnetic field that would introduce same
strength rotation is $B^{eq}=\Delta E_{(\mu)}/\mu_{\mu}\simeq 1.1\,{\rm mGs}$,
where $\mu_{\mu}$ is the muon intrinsic magnetic moment. The smallest DC
magnetic fields detected with $\mu$SR techniques are about 0.1 Gs
RevModPhys.69.1119 . It is then clear that the muon spin force effect needs to
be measured as the shift of muon spin precession frequency at low/moderate
values of the magnetic field. In the past, studies of the muon magnetic moment
stopped in liquid helium in the field of 66 Gs have brought a relative
precision of $\mathcal{O}(10^{-4})$ PhysRevLett.33.572 . One order of
magnitude improvement in accuracy would reveal or exclude $\Delta E_{(\mu)}$.

A suggested $\mu$SR setup for an experiment should involve the following
ingredients:

  * •

Uniform vertical magnetic field in the place where muons are stopped is
required. A high degree of uniformity across the stopping target should be
possible to achieve. The optimal magnetic field values that will afford
precision measurement of the muon spin precession frequency $\omega_{(\mu)}$
will correspond to $\Delta\psi\propto\mathcal{O}(1)$. This corresponds to the
magnetic fields in excess of 10 Gs, and we suggest the range of 5-to-50 Gs,
and a possibility of reversing polarity.

  * •

A co-magnetometer species that monitors the magnetic field and its stability
is needed. As in the muon $g-2$ experiment, one can use proton (water-based)
or ${}^{3}{\rm He}$ co-magnetometers Aoyama:2020ynm ; PhysRevLett.124.223001 .

  * •

The highest accuracy results can be achieved in the gaseous stopping target,
that would also afford atomic co-magnetometry. This will minimize the medium
feedback on the field and on the muon spin. Low densities of such targets
necessitate low-momentum of incoming muons. The upcoming high-intensity muon
beams at PSI Aiba:2021bxe  will provide $10^{5}$/sec muon beams with $\sim 4$
MeV momentum (reduced relative to a usual 29 MeV for surface muons) that will
stop within 1 cm inside a gaseous target.

  * •

For a continuous muon beam, $10^{5}$/sec intensity appears limiting to
minimize the muon pile-up, while the pulsed muon beam with a sub-$\mu$s pulse
duration can afford higher muon intensities.

If both muon and co-magnetometer spin precession frequencies are measured,
their ratio $r_{\mu/p}=\omega_{(\mu)}/\omega_{(p)}$ should not depend on the
value of the magnetic field, in the absence of the spin force. Therefore,
performing the measurement at several different values of the magnetic field
should isolate the signal of interest. The following quantity,

$-B\frac{dr_{\mu/p}}{dB}=\frac{\Delta E_{(\mu)}}{2\mu_{p}B}=3.4\times
10^{-4}\times\frac{\Delta E_{(\mu)}}{6\times 10^{-14}\,{\rm
eV}}\times\frac{{\rm 10\,Gs}}{B},$ (21)

provides a figure of merit that needs to be achieved in this suggested
experiment. Here we used proton spin precession as reference. The effect is
obviously enhanced at small values of the magnetic field, and reverses sign
with the reversal of $B$ direction. The effect also changes sign for
$\mu^{-}$, which can be tested experimentally by studying negative muon spin
precession in muonic atoms (if challenges due to lower fluxes of $\mu^{-}$ and
muon depolarization during the cascade can be overcome).

While $\mu$SR setting appears to be the most promising place for a muon spin
force search, it is not the only way of achieving this goal. Other
possibilities include precision studies of muonium hyperfine structure in
low/medium vertically directed magnetic field Kanda:2020mmc .

Conclusions — In this _Letter_ , we have investigated terrestrial probes of
the muon spin force mediated by a light $CP$-violating scalar field. In the
first part, we derived indirect constraints on the muon spin force from Hg and
Xe isotope experiments, with the muon appearing only inside quantum loops. We
found that the pseudoscalar coupling explaining the muon $g-2$ anomaly is
well-excluded by these experiments. Although the axial-vector coupling is also
excluded by an order of magnitude with our modeling of nuclear physics, this
constraint is subject to huge nuclear uncertainty, casting doubts on its
robustness. This motivates us to suggest, in the second part, a new $\mu$SR-
type experiment to look for the muon spin force directly. Based on our
estimation, an experiment testing the muon spin force required for the $g-2$
anomaly is entirely feasible within the current technology. Our study thus
offers an interesting experimental opportunity for future muon precision
physics facilities.

###### Acknowledgements.

Acknowledgments — Y.E. and M.P. are supported in part by U.S. Department of
Energy Grant No. desc0011842. We would like to acknowledge helpful
communications with Drs. D. Budker and M. Romalis. Y.E. would like to thank
the Aspen Center for Physics (supported by National Science Foundation grant
PHY-2210452) where part of this work was done. M.P. would like to thank
Perimeter Institute for hospitality during the completion of this project.
T.G. would like to thank the Institute for Nuclear Theory and the Perimeter
Institute for the visits and discussions about this project. The Feynman
diagrams in this paper are drawn with TikZ-Feynman Ellis:2016jkw .

## References

  * (1) M. Aiba et al., “Science Case for the new High-Intensity Muon Beams HIMB at PSI,” arXiv:2111.05788 [hep-ex].
  * (2) C. Group Collaboration, M. Aoki et al., “A New Charged Lepton Flavor Violation Program at Fermilab,” in Snowmass 2021.  3, 2022.  arXiv:2203.08278 [hep-ex].
  * (3) N. Kawamura et al., “New concept for a large-acceptance general-purpose muon beamline,” PTEP 2018 no. 11, (2018) 113G01.
  * (4) J. W. G. Thomason, “The ISIS Spallation Neutron and Muon Source—The first thirty-three years,” Nucl. Instrum. Meth. A 917 (2019) 61–67.
  * (5) T. Aoyama et al., “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887 (2020) 1–166, arXiv:2006.04822 [hep-ph].
  * (6) S. Borsanyi et al., “Leading hadronic contribution to the muon magnetic moment from lattice QCD,” Nature 593 no. 7857, (2021) 51–55, arXiv:2002.12347 [hep-lat].
  * (7) Muon g-2 Collaboration, B. Abi et al., “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm,” Phys. Rev. Lett. 126 no. 14, (2021) 141801, arXiv:2104.03281 [hep-ex].
  * (8) M. Cè et al., “Window observable for the hadronic vacuum polarization contribution to the muon g-2 from lattice QCD,” Phys. Rev. D 106 no. 11, (2022) 114502, arXiv:2206.06582 [hep-lat].
  * (9) Extended Twisted Mass Collaboration, C. Alexandrou et al., “Lattice calculation of the short and intermediate time-distance hadronic vacuum polarization contributions to the muon magnetic moment using twisted-mass fermions,” Phys. Rev. D 107 no. 7, (2023) 074506, arXiv:2206.15084 [hep-lat].
  * (10) Fermilab Lattice, HPQCD,, MILC Collaboration, A. Bazavov et al., “Light-quark connected intermediate-window contributions to the muon g-2 hadronic vacuum polarization from lattice QCD,” Phys. Rev. D 107 no. 11, (2023) 114514, arXiv:2301.08274 [hep-lat].
  * (11) CMD-3 Collaboration, F. V. Ignatov et al., “Measurement of the $e^{+}e^{-}\to\pi^{+}\pi^{-}$ cross section from threshold to 1.2 GeV with the CMD-3 detector,” arXiv:2302.08834 [hep-ex].
  * (12) A. Czarnecki and W. J. Marciano, “The Muon anomalous magnetic moment: A Harbinger for ’new physics’,” Phys. Rev. D 64 (2001) 013014, arXiv:hep-ph/0102122.
  * (13) R. Janish and H. Ramani, “Muon g-2 and EDM experiments as muonic dark matter detectors,” Phys. Rev. D 102 (2020) 115018, arXiv:2006.10069 [hep-ph].
  * (14) P. Agrawal, D. E. Kaplan, O. Kim, S. Rajendran, and M. Reig, “Searching for axion forces with precision precession in storage rings,” arXiv:2210.17547 [hep-ph].
  * (15) H. Davoudiasl and R. Szafron, “Muon g-2 and a Geocentric New Field,” arXiv:2210.14959 [hep-ph].
  * (16) M. Pospelov and M. Romalis, “Lorentz invariance on trial,” Phys. Today 57N7 (2004) 40–46.
  * (17) J. E. Moody and F. Wilczek, “NEW MACROSCOPIC FORCES?,” Phys. Rev. D 30 (1984) 130.
  * (18) Y. Ema, T. Gao, and M. Pospelov, “Improved Indirect Limits on Muon Electric Dipole Moment,” Phys. Rev. Lett. 128 no. 13, (2022) 131803, arXiv:2108.05398 [hep-ph].
  * (19) Y. Ema, T. Gao, and M. Pospelov, “Reevaluation of heavy-fermion-induced electron EDM at three loops,” Phys. Lett. B 835 (2022) 137496, arXiv:2207.01679 [hep-ph].
  * (20) Y. K. Semertzidis et al., “Sensitive search for a permanent muon electric dipole moment,” in KEK International Workshop on High Intensity Muon Sources (HIMUS 99).  12, 1999.  arXiv:hep-ph/0012087.
  * (21) Muon (g-2) Collaboration, G. W. Bennett et al., “An Improved Limit on the Muon Electric Dipole Moment,” Phys. Rev. D 80 (2009) 052008, arXiv:0811.1207 [hep-ex].
  * (22) M. Pospelov, A. Ritz, and M. B. Voloshin, “Bosonic super-WIMPs as keV-scale dark matter,” Phys. Rev. D 78 (2008) 115012, arXiv:0807.3279 [hep-ph].
  * (23) V. Flambaum, S. Lambert, and M. Pospelov, “Scalar-tensor theories with pseudoscalar couplings,” Phys. Rev. D 80 (2009) 105021, arXiv:0902.3217 [hep-ph].
  * (24) M. A. Fedderke, J. O. Thompson, R. Cervantes, B. Giaccone, R. Harnik, D. E. Kaplan, S. Posen, and S. Rajendran, “MAGPI: Measurement of Axion Gradients with Photon Interferometry,” arXiv:2304.11261 [hep-ph].
  * (25) M. Abe et al., “A New Approach for Measuring the Muon Anomalous Magnetic Moment and Electric Dipole Moment,” PTEP 2019 no. 5, (2019) 053C02, arXiv:1901.03047 [physics.ins-det].
  * (26) B. J. Venema, P. K. Majumder, S. K. Lamoreaux, B. R. Heckel, and E. N. Fortson, “Search for a coupling of the earth’s gravitational field to nuclear spins in atomic mercury,” Phys. Rev. Lett. 68 (Jan, 1992) 135–138. https://link.aps.org/doi/10.1103/PhysRevLett.68.135.
  * (27) S. Zhang, Z. Ba, D. Ning, N. Zhai, Z. Lu, and D. Sheng, “Search for Spin-Dependent Gravitational Interactions at Earth Range,” Phys. Rev. Lett. 130 no. 20, (2023) 201401, arXiv:2303.10352 [gr-qc].
  * (28) W. A. Terrano and M. V. Romalis, “Comagnetometer probes of dark matter and new physics,” Quantum Sci. Technol. 7 no. 1, (2022) 014001, arXiv:2106.09210 [hep-ph].
  * (29) D. F. Jackson Kimball, J. Dudley, Y. Li, D. Patel, and J. Valdez, “Constraints on long-range spin-gravity and monopole-dipole couplings of the proton,” Phys. Rev. D 96 no. 7, (2017) 075004, arXiv:1707.00745 [physics.atom-ph]. [Erratum: Phys.Rev.D 107, 019903 (2023)].
  * (30) B. R. Heckel, E. G. Adelberger, C. E. Cramer, T. S. Cook, S. Schlamminger, and U. Schmidt, “Preferred-Frame and CP-Violation Tests with Polarized Electrons,” Phys. Rev. D 78 (2008) 092006, arXiv:0808.2673 [hep-ex].
  * (31) V. F. Dmitriev and R. A. Sen’kov, “Schiff moment of the mercury nucleus and the proton dipole moment,” Phys. Rev. Lett. 91 (2003) 212303, arXiv:nucl-th/0306050.
  * (32) V. V. Flambaum, M. Pospelov, A. Ritz, and Y. V. Stadnik, “Sensitivity of EDM experiments in paramagnetic atoms and molecules to hadronic CP violation,” Phys. Rev. D 102 no. 3, (2020) 035001, arXiv:1912.13129 [hep-ph].
  * (33) S. G. Karshenboim, D. McKeen, and M. Pospelov, “Constraints on muon-specific dark forces,” Phys. Rev. D 90 no. 7, (2014) 073004, arXiv:1401.6154 [hep-ph]. [Addendum: Phys.Rev.D 90, 079905 (2014)].
  * (34) Flavour Lattice Averaging Group (FLAG) Collaboration, Y. Aoki et al., “FLAG Review 2021,” Eur. Phys. J. C 82 no. 10, (2022) 869, arXiv:2111.09849 [hep-lat].
  * (35) D. F. Jackson Kimball, “Nuclear spin content and constraints on exotic spin-dependent couplings,” New J. Phys. 17 no. 7, (2015) 073008, arXiv:1407.2671 [physics.atom-ph].
  * (36) P. Klos, J. Menéndez, D. Gazit, and A. Schwenk, “Large-scale nuclear structure calculations for spin-dependent WIMP scattering with chiral effective field theory currents,” Phys. Rev. D 88 no. 8, (2013) 083516, arXiv:1304.7684 [nucl-th]. [Erratum: Phys.Rev.D 89, 029901 (2014)].
  * (37) A. Amato, “Heavy-fermion systems studied by $\mu$sr technique,” Rev. Mod. Phys. 69 (Oct, 1997) 1119–1180. https://link.aps.org/doi/10.1103/RevModPhys.69.1119.
  * (38) T. W. Crane, D. E. Casperson, H. Chang, V. W. Hughes, H. F. Kaspar, B. Lovett, A. Schiz, P. Souder, R. D. Stambaugh, G. z. Putlitz, and J. R. Kane, “Behavior of positive muons in liquid helium,” Phys. Rev. Lett. 33 (Sep, 1974) 572–574. https://link.aps.org/doi/10.1103/PhysRevLett.33.572.
  * (39) M. Farooq, T. Chupp, J. Grange, A. Tewsley-Booth, D. Flay, D. Kawall, N. Sachdeva, and P. Winter, “Absolute magnetometry with ${}^{3}\mathrm{He}$,” Phys. Rev. Lett. 124 (Jun, 2020) 223001. https://link.aps.org/doi/10.1103/PhysRevLett.124.223001.
  * (40) S. Kanda et al., “New precise spectroscopy of the hyperfine structure in muonium with a high-intensity pulsed muon beam,” Phys. Lett. B 815 (2021) 136154, arXiv:2004.05862 [hep-ex].
  * (41) J. Ellis, “TikZ-Feynman: Feynman diagrams with TikZ,” Comput. Phys. Commun. 210 (2017) 103–123, arXiv:1601.05437 [hep-ph].
  * (42) V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Calculations in External Fields in Quantum Chromodynamics. Technical Review,” Fortsch. Phys. 32 (1984) 585.
  * (43) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii, QUANTUM ELECTRODYNAMICS, vol. 4 of Course of Theoretical Physics.  Pergamon Press, Oxford, 1982.
  * (44) L. D. Landau and E. M. Lifshits, Quantum Mechanics: Non-Relativistic Theory, vol. v.3 of Course of Theoretical Physics.  Butterworth-Heinemann, Oxford, 1991.
  * (45) N. J. Stone, “Table of nuclear magnetic dipole and electric quadrupole moments,” Atom. Data Nucl. Data Tabl. 90 (2005) 75–176.

Supplemental Material

In this supplemental material, we provide details on our one-loop, nuclear
magnetic field, and two-loop calculation.

##  Appendix S1 Coupling to the external electromagnetic field at one-loop

Here we start from the muon spin force (3) and derive the one-loop effective
Lagrangian (10).

We first calculate the effect induced by $g_{A}^{(\mu)}$. Following
Novikov:1983gd , the one-loop effective action in the external electromagnetic
field can be written as

$\displaystyle S_{A,\mathrm{eff}}$
$\displaystyle=-\frac{g_{A}^{(\mu)}}{2m_{\mu}}\partial_{\nu}\phi\times\text{Tr}\left[\frac{i}{i\not{D}-m_{\mu}}\gamma^{\nu}\gamma^{5}\right],$
(S1)

where the trace is over the spinor indices and the spacetime, and $D_{\alpha}$
is the covariant derivative, which is given by

$iD_{\alpha}=i\partial_{\alpha}+eA_{\alpha},$ (S2)

and $e$ is the charge of the positron. Using
$(i\not{D})(i\not{D})=(iD)^{2}+\frac{1}{2}e(\sigma F)$, where $(\sigma
F)=\sigma_{\alpha\beta}F^{\alpha\beta}$, we expand (S1) in terms of the
external field. Up to the order of our interest, we get

$S_{A,eff}=-i\frac{g_{A}^{(\mu)}}{2m_{\mu}}\partial_{\nu}\phi\times\left(T_{0}^{\nu}-\frac{e}{2}T_{1}^{\nu}+\frac{e^{2}}{4}T_{2}^{\nu}\right),$
(S3)

where

$\displaystyle
T_{0}^{\nu}=\text{Tr}\left[\frac{1}{(iD)^{2}-m_{\mu}^{2}}(i\not{D})\gamma^{\nu}\gamma^{5}\right],$
(S4) $\displaystyle
T_{1}^{\nu}=\text{Tr}\left[\frac{1}{(iD)^{2}-m_{\mu}^{2}}(\sigma
F)\frac{1}{(iD)^{2}-m_{\mu}^{2}}(i\not{D})\gamma^{\nu}\gamma^{5}\right],$ (S5)
$\displaystyle
T_{2}^{\nu}=\text{Tr}\left[\frac{1}{(iD)^{2}-m_{\mu}^{2}}(\sigma
F)\frac{1}{(iD)^{2}-m_{\mu}^{2}}(\sigma
F)\frac{1}{(iD)^{2}-m_{\mu}^{2}}(i\not{D})\gamma^{\nu}\gamma^{5}\right].$ (S6)

In the expression above, we omitted terms with $m_{\mu}$ in the numerator
because they vanish after taking the trace over spinor indices. $T_{0}^{\nu}$
also vanishes due to the trace property of gamma matrices.

To evaluate $T_{1}^{\nu}$, we use the fact that the quantity
$\text{Tr}\left[f(F)\frac{1}{(iD)^{2}-m_{\mu}^{2}}\right]$ does not depend on
the momentum $iD_{\alpha}$, where $f(F)$ is an arbitrary function of the
electromagnetic field. By shifting the momentum by $q_{\alpha}$ and focusing
on terms linear in $q_{\alpha}$, we get

$\text{Tr}\left[f(F)\frac{1}{(iD)^{2}-m^{2}}(iD)_{\alpha}\frac{1}{(iD)^{2}-m^{2}}\right]=0.$
(S7)

Therefore,

$T_{1}^{\nu}=\text{Tr}\left[F^{\alpha\beta}\frac{1}{(iD)^{2}-m_{\mu}^{2}}(iD)_{\gamma}\frac{1}{(iD)^{2}-m_{\mu}^{2}}\right]\text{tr}\left[\sigma_{\alpha\beta}\gamma^{\gamma}\gamma^{\nu}\gamma^{5}\right]=0.$
(S8)

To evaluate $T_{2}^{\nu}$ up to $(DF)F$ order, the following identities turn
out to be useful

$\displaystyle\left[(iD)^{2},\sigma F\right]$
$\displaystyle=\left[iD_{\alpha},\left[iD^{\alpha},\sigma
F\right]\right]+2\left[iD_{\alpha},\sigma
F\right]iD^{\alpha}=2iD_{\alpha}(\sigma F)iD^{\alpha}+\cdots,$ (S9)
$\displaystyle\left[(iD)^{2},iD_{\nu}\right]$
$\displaystyle=\left[iD_{\alpha},\left[iD^{\alpha},iD_{\nu}\right]\right]+2\left[iD_{\alpha},iD_{\nu}\right]iD^{\alpha}=-eD^{\alpha}F_{\alpha\nu}+2ieF_{\alpha\nu}iD^{\alpha}=0+\cdots,$

where terms that introduce extra electromagnetic fields or more than one
derivative on the electromagnetic field are not of our interest. Now
$T_{2}^{\nu}$ can be written as

$\displaystyle T_{2}^{\nu}$
$\displaystyle=\text{Tr}\left[\frac{1}{(iD)^{2}-m^{2}}(\sigma
F)\frac{1}{(iD)^{2}-m^{2}}(\sigma
F)\frac{1}{(iD)^{2}-m^{2}}(iD)^{\alpha}\gamma_{\alpha}\gamma^{\nu}\gamma^{5}\right]$
(S10)
$\displaystyle=\text{Tr}\left[\frac{1}{\left((iD)^{2}-m^{2}\right)^{3}}(\sigma
F)(\sigma F)(iD)^{\alpha}\gamma_{\alpha}\gamma^{\nu}\gamma^{5}\right]$
$\displaystyle\quad+\text{Tr}\left[\frac{1}{\left((iD)^{2}-m^{2}\right)^{4}}\left(2\left[(iD)^{2},\sigma
F\right](\sigma F)+(\sigma F)\left[(iD)^{2},\sigma
F\right]\right)(iD)^{\alpha}\gamma_{\alpha}\gamma^{\nu}\gamma^{5}\right]+\cdots.$

By looking at terms linear in $q$ in
$\text{Tr}\left[f(F)\frac{1}{\left((iD-q)^{2}-m^{2}\right)^{2}}\right]$, one
can show that the first term in (S10) vanishes up to a total derivative. The
second term already has enough powers of $F$ and $DF$, so everything commutes.
Using (S9), we arrive at the following result

$\displaystyle T_{2}^{\nu}$
$\displaystyle=\frac{ig^{\rho\zeta}}{2}\text{Tr}\left[\frac{(iD)^{2}}{((iD)^{2}-m^{2})^{4}}\right]\left[\left(D_{\rho}F^{\alpha\beta}F^{\gamma\delta}+D_{\rho}(F^{\alpha\beta}F^{\gamma\delta})\right)\right]\text{tr}\left[\sigma_{\alpha\beta}\sigma_{\gamma\delta}\gamma_{\zeta}\gamma^{\nu}\gamma^{5}\right]=-\frac{i}{6\pi^{2}m_{\mu}^{2}}\int
d^{4}x{\tilde{F}_{\alpha}}^{~{}\nu}D_{\beta}F^{\alpha\beta},$ (S11)

where we performed integration by parts and neglected terms proportional to
$\partial_{\alpha}\partial_{\beta}\phi$. Putting it back to (S3), we get

$S_{A,\mathrm{eff}}=-\frac{e^{2}g_{A}^{(\mu)}\partial^{\nu}\phi}{48\pi^{2}m_{\mu}^{3}}\int
d^{4}x\tilde{F}_{\alpha\nu}D_{\beta}F^{\alpha\beta}.$ (S12)

We have checked that we reproduce the same result by the direct Feynman
diagrammatic computation of the triangle diagrams. The effect induced by
$g_{P}^{(\mu)}$ is obtained directly from the anomaly equation:

$\partial_{\nu}(\bar{\mu}\gamma^{\nu}\gamma^{5}\mu)-2m_{\mu}\bar{\mu}i\gamma^{5}\mu=-\frac{e^{2}}{8\pi^{2}}F_{\alpha\beta}\tilde{F}^{\alpha\beta},$
(S13)

multiplying both sides by $\frac{g_{P}^{(\mu)\phi}}{2m_{\mu}}$ and performing
integration by parts leads to

$-g_{P}^{(\mu)}\phi\bar{\mu}i\gamma^{5}\mu=-\frac{e^{2}g_{P}^{(\mu)}}{16\pi^{2}m_{\mu}}\phi
F_{\alpha\beta}\tilde{F}^{\alpha\beta}+\frac{g_{P}^{(\mu)}\partial_{\nu}\phi}{2m_{\mu}}\bar{\mu}\gamma^{\nu}\gamma^{5}\mu.$
(S14)

As is shown in (S12), the contribution from
$\bar{\mu}\gamma^{\nu}\gamma^{5}\mu$ is suppressed by the muon mass, so at
leading order, only the first term on the right-hand side contributes, so we
get

$\mathcal{L}_{P,\mathrm{eff}}=-\frac{e^{2}g_{P}^{(\mu)}}{16\pi^{2}m_{\mu}}\phi
F_{\alpha\beta}\tilde{F}^{\alpha\beta}.$ (S15)

Again, we have checked that we obtain the same result by computing the
triangle diagrams directly.

##  Appendix S2 Nuclear magnetic field

Here we provide details on our model of the magnetic field in (11) which leads
to (14) and (15). The non-relativistic eigenfunction of a spin 1/2 particle in
spherical coordinates can be written as Berestetskii:1982qgu

$\psi(r,\theta,\varphi)=R_{nl}(r)\Omega_{jlm}(\theta,\varphi),$ (S16)

where

$\displaystyle\Omega_{l+1/2,l,m}=\left(\begin{array}[]{c}\sqrt{\frac{j+m}{2j}}Y_{l,m-1/2}\\\
\sqrt{\frac{j-m}{2j}}Y_{l,m+1/2}\end{array}\right),\quad\Omega_{l-1/2,l,m}=\left(\begin{array}[]{c}-\sqrt{\frac{j-m+1}{2j+2}}Y_{l,m-1/2}\\\
\sqrt{\frac{j+m+1}{2j+2}}Y_{l,m+1/2}\end{array}\right),$ (S17)

$Y_{l,m}$ are spherical harmonics, and the radial function $R_{nl}(r)$ is
determined by

$\frac{d^{2}(rR_{nl})}{dr^{2}}+2m(E_{n}-U)rR_{nl}-\frac{l(l+1)}{r^{2}}rR_{nl}=0.$
(S18)

They are normalized by

$\int
d\Omega\left(\Omega_{j^{\prime}l^{\prime}m^{\prime}}^{\dagger}\Omega_{jlm}\right)=\delta_{jj^{\prime}}\delta_{ll^{\prime}}\delta_{mm^{\prime}},\quad\int
drr^{2}R_{nl}R_{n^{\prime}l}=\delta_{nn^{\prime}}.$ (S19)

The magnetic field generated by the magnetic moment of the valance neutron is
given by Landau:1991wop

$\bm{B}_{n}(\bm{r})=-\frac{\mu_{n}}{4\pi}\int
d^{3}r^{\prime}\left[\bm{\nabla}^{\prime}\times\left(\psi^{\dagger}(\bm{r}^{\prime})\bm{\sigma}\psi(\bm{r}^{\prime})\right)\right]\times\bm{\nabla}\frac{1}{|\bm{r-r^{\prime}}|},$
(S20)

where $\mu_{n}$ is the magnetic moment of the neutron and $\bm{\sigma}$ are
Pauli matrices.

While it is possible to compute the magnetic field directly for any given
eigenstate, a much simpler expression can be obtained when we average out the
angular dependence, this is valid as long as the electric potential $A_{0}$ is
independent of the angle. The angular averaged magnetic field is calculated to
be

$\displaystyle\bm{\mathcal{B}}_{n}(r)$ $\displaystyle=\frac{1}{4\pi}\int
d\Omega\bm{B}_{n}(\bm{r})$ (S21) $\displaystyle=\frac{\mu_{n}}{(4\pi)^{2}}\int
d^{3}r^{\prime}\left[\bm{\nabla}^{\prime}\times\left(\psi^{\dagger}(\bm{r}^{\prime})\bm{\sigma}\psi(\bm{r}^{\prime})\right)\right]\times\int
d\Omega\frac{\bm{r-r^{\prime}}}{|\bm{r-r^{\prime}}|^{3}}$
$\displaystyle=-\frac{\mu_{n}}{4\pi}\int
d^{3}r^{\prime}\left[\bm{\nabla}^{\prime}\times\left(\psi^{\dagger}(\bm{r}^{\prime})\bm{\sigma}\psi(\bm{r}^{\prime})\right)\right]\times\hat{\bm{r}}^{\prime}\frac{\Theta(r^{\prime}-r)}{r^{\prime
2}},$

where $\Theta(x)$ is the Heaviside step function. The integration over
$d\Omega^{\prime}$ involves only spherical harmonics and the result can be
parameterized as

$\int
d\Omega^{\prime}\left[\bm{\nabla}^{\prime}\times\left(\psi^{\dagger}(\bm{r}^{\prime})\bm{\sigma}\psi(\bm{r}^{\prime})\right)\right]\times\hat{\bm{r}}^{\prime}=\left(c_{1}\frac{R_{nl}^{2}}{r^{\prime}}+c_{2}\frac{\text{d}R_{nl}^{2}}{\text{d}r^{\prime}}\right)\hat{\bm{I}},$
(S22)

where $\hat{\bm{I}}$ is the direction of the nuclear spin, and $c_{1}$ and
$c_{2}$ are constants which depends on the angular quantum numbers. For the
cases of our interest, they are

$(c_{1},c_{2})_{{}^{129}\text{Xe}}=\left(0,\frac{2}{3}\right),\quad(c_{1},c_{2})_{{}^{131}\text{Xe}}=\left(-\frac{6}{5},-\frac{4}{5}\right),\quad(c_{1},c_{2})_{{}^{199}\text{Hg}}=\left(-\frac{4}{3},-\frac{2}{3}\right),\quad(c_{1},c_{2})_{{}^{201}\text{Hg}}=\left(\frac{2}{5},\frac{4}{5}\right).$
(S23)

Now (S21) simplifies to

$\bm{\mathcal{B}}_{n}(r)=-\frac{\mu_{n}\hat{\bm{I}}}{4\pi}\left[c_{1}\left(\int_{r}^{\infty}dr^{\prime}\frac{R_{nl}^{2}(r^{\prime})}{r^{\prime}}\right)-c_{2}R_{nl}^{2}(r)\right]$
(S24)

The magnetic moment of the nucleus generated by the valance neutron is
obtained from the relation $\int
d^{3}r\bm{B}_{n}(\bm{r})=\frac{2}{3}\mu_{N}\hat{\bm{I}}$, it is related to
$c_{1}$ and $c_{2}$ by

$\mu_{N}=-\frac{\mu_{n}}{2}(c_{1}-3c_{2}).$ (S25)

The difference between the valance neutron contribution to nucleus magnetic
moment $\mu_{N}$ and the experimental value $\mu_{N,\text{exp}}$ (taken _e.g._
from Stone:2005rzh ) is modeled by a uniformly magnetic polarized core,
averaging over all directions, we get

$\bm{\mathcal{B}}_{\text{core}}(r)=\frac{\mu_{N,\text{exp}}-\mu_{N}}{2\pi
R_{N}^{3}}\hat{\bm{I}}\Theta(R_{N}-r).$ (S26)

To get numerical results, we numerically solve the radial wave function taking
$U$ to be a Woods-Saxon potential using the parameters in Dmitriev:2003sc .
The numbers in (14) and (15) are then obtained from (10) using

$H=-\int d^{3}r\mathcal{L}_{\mathrm{eff}}.$ (S27)

##  Appendix S3 Coupling to light fermions at two-loop

Figure S1: Two-loop diagrams leading to Eq. (17). The upper fermion line
represents the light fermion and the lower fermion loop is formed by the muon.

Here we discuss the computation that leads to (17). The contribution comes
from the two diagrams shown in Fig. S1. We are interested in the UV divergent
part of the diagram, so all the external momentum as well as the mass of the
light fermion are set to 0. Following the procedure described in Ema:2022wxd ,
we write the amplitude as

$\mathcal{M}=\frac{g_{A}^{(\mu)}\partial_{\alpha}\phi}{2m_{\mu}}\bar{\psi}_{i}\mathcal{M}^{\alpha}\psi_{i},$
(S28)

where $\mathcal{M}_{\alpha}$ can be factored into a scalar integral and a term
describing its Lorentz and spinor structure:

$\mathcal{M}^{\alpha}=\frac{1}{24d(d-1)(d-2)(d-3)}\text{tr}\left[\mathcal{M}_{\kappa}\epsilon^{\kappa\lambda\nu\xi}\gamma_{\lambda}\gamma_{\nu}\gamma_{\xi}\right]\times\epsilon^{\alpha\beta\gamma\delta}\gamma_{\beta}\gamma_{\gamma}\gamma_{\delta}.$
(S29)

After the integration-by-parts reduction, the two-loop integral reduces to two
one-loop integrals, which are computed using the standard Euclidean space
technique. The divergent part of our result is given by

$\mathcal{M}=-\frac{g_{A}^{(\mu)}}{2m_{\mu}}\frac{3}{4}\left(\frac{\alpha}{\pi}\right)^{2}Q_{i}^{2}\ln\left(\frac{\Lambda_{\text{UV}}^{2}}{\Lambda_{\mathrm{IR}}^{2}}\right)\partial_{\alpha}\phi\times\bar{\psi}_{i}\gamma^{\alpha}\gamma^{5}\psi_{i}.$
(S30)

In the case of the electron, the IR scale is given by
$\Lambda_{\mathrm{IR}}=m_{\mu}$. In the case of the up and down quarks, the
QCD scale kicks in before the energy scale reaches $m_{\mu}$, below which the
free quark picture is no longer valid. Therefore we may take the IR scale as
the hadronic scale, _e.g._ , $\Lambda_{\mathrm{IR}}=m_{p}$.