# Probing Scalar and Pseudoscalar Solutions of the $g-2$ Anomaly

Fayez Abu-Ajamieh

► To cite this version:

Fayez Abu-Ajamieh. Probing Scalar and Pseudoscalar Solutions of the $g-2$ Anomaly. Advances in High Energy Physics, 2020, 2020, pp.1751534. 10.1155/2020/1751534. hal-02885943

HAL Id: hal-02885943

https://hal.science/hal-02885943

Submitted on 30 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Probing Scalar and Pseudoscalar Solutions to the g-2 Anomaly

Fayez Abu-Ajamieh*

Department of Physics, UC Davis, Davis, California, 95616, USA

LUPM UMR5299, Université de Montpellier, 34095 Montpellier, France

CPPM, Aix-Marseille Université, Marseille, France

I investigate a class of models with scalar and pseudoscalar solutions to the g-2 anomaly for both the muon and the electron over the mass range of perturbativity ($m_\phi \lesssim 50$ GeV), with Yukawa couplings proportional to the lepton's mass. In particular, I investigate the constraints from BaBar, beam dump experiments, Z decay measured quantities, LEP mono-$\gamma$ searches, ee $\to \tau\tau(\gamma)$ searches, and solar and Horizontal Branch (HB) stars bounds. For a pseudoscalar, I find that no region in the parameters space can simultaneously provide a solution for both the electron and the muon anomalies while maintaining the required form of the couplings, and therefore the pseudoscalar solution is disfavored. On the other hand, I find for the scalar case that there is an open window above $\sim 30$ MeV in the allowed region, but with significant tension with experiment for the region $m_\phi \gtrsim 10$ GeV. In addition, there is a smaller window between $\sim 350$ KeV and 1 MeV that is not ruled out by cosmological observations. Part of the first open window is expected to be covered by the proposed NA64 experiment. Similar analysis can be readily applied to other proposed solutions to the anomaly, such as solutions with Z' or with the dark photon.

# I. INTRODUCTION

An exciting piece of evidence for the existence of physics Beyond the Standard Model (BSM) is the discrepancy between the predicted and the measured values of the muon anomalous magnetic moment $a_\mu \equiv (g_\mu - 2)/2$. The current measured value [1-4] shows a $3.5\sigma$ discrepancy compared with the SM prediction [5-7]:

$$ \Delta a_{\mu} = a_{\mu}^{\text{Exp}} - a_{\mu}^{\text{SM}} = 273 \pm 80 \times 10^{-11}. \quad (1) $$

A similar less significant discrepancy of about $1.1\sigma$ was also observed for the electron [8]:

$$ \Delta a_e = a_e^{\text{Exp}} - a_e^{\text{SM}} = -91 \pm 82 \times 10^{-14}. \quad (2) $$

Although both discrepancies fall short of the $5\sigma$ limit required to confirm their existence, they nonetheless pose tantalizing hints for physics BSM. In addition, current experiments at Fermilab [9, 10] and at the J-PARC E34 collaboration [11, 12] are expected to yield improved experimental results in the near future.

New physics explanations of this anomaly include (see [13] for a comprehensive review) supersymmetry (see [14] for a review), a light Z' boson [15-23] (also see [24] for a review), a scalar contribution within the framework of the 2 Higgs Doublet Model (2HDM) [25-32], additional fermions [33], leptoquarks [34, 35], and the dark photon [36].

Recently, there have been proposed solutions to this anomaly through a scalar [5] or a pseudoscalar Axion-Like Particle (ALP)[6] in a general framework. In this

short paper, I will investigate the viability of these solutions, explore the relevant experimental limits, and highlight the experimental probes for their discovery for the mass range of their validity.

For the case of a pseudoscalar, the effective interaction with photons and fermions can be parametrized by:<fnref n="1" />

$$ \mathcal{L} = \frac{1}{4} g_{\phi\gamma\gamma} \phi F_{\mu\nu} \tilde{F}^{\mu\nu} + i Y_{\phi ll} \phi \bar{\psi} \gamma_5 \psi, \quad (3) $$

where $g_{\phi\gamma\gamma}$ is a dimensionful coupling, $Y_{\phi ll}$ is a dimensionless Yukawa coupling and $F^{\mu\nu}$, $\tilde{F}^{\mu\nu}$ are the magnetic field strength tensor and its dual respectively. For a scalar, $\tilde{F}^{\mu\nu}$ is replaced with $F^{\mu\nu}$ and there is no $i\gamma_5$ in the second term. Since $|\frac{\Delta a_e}{\Delta a_\mu}| \sim \frac{Y_{\phi ee}^2}{Y_{\phi\mu\mu}^2}$ (see Eqs. 5 and 6 below), and we can see from Eqs. 1 and 2 that within the allowed range of uncertainties, we could obtain $|\frac{\Delta a_e}{\Delta a_\mu}| \approx (\frac{m_e}{m_\mu})^4$ $^2$, then we are motivated to define the Yukawa couplings to be proportional to the lepton mass:

$$ Y_{\phi ll} = \frac{m_l}{v} \equiv m_l g_{\phi ll}, \quad (4) $$

where $v \equiv g_{\phi ll}^{-1}$ is some model-dependent energy scale that is universal for all leptons, such as the axion decay

<footnote n="1">¹ This interaction can be viewed as an effective theory of a UV-complete model. One possible UV-completion that is consistent with the SM EW theory was introduced in [5]. Such a model could lead to lepton-flavor violation through terms like $Y_{\phi ij} \phi \bar{l}_i l_j$, $i \neq j$. We ignore this possibility in this paper as it will not affect the interaction in Eq. 3 and as it was studied in detail in [5]. The interested reader is instructed to refer to [5] for detailed analysis.</footnote>

<footnote n="2">² Notice that from Eqs. 5 and 6, $\Delta a_l \sim Y_{\phi ll}^2 r^{-2}$, where $r = m_\phi/m_l$. The assumption in Eq. 4 makes $\Delta a_l \sim m_l^4$, and thus $\frac{\Delta a_e}{\Delta a_\mu} \approx (\frac{m_e}{m_\mu})^4$.</footnote>

<footnote>* abuajamieh@ucdavis.edu, fayez.abu-ajamieh@umontpellier.fr</footnote>
constant or the radion constant. The coupling of the form given in Eq. 4 has the additional advantage in that it arises in many UV-comepltions as Psuedo Nambu Goldstone Bosons (PNGB), such as axion models [37], Left-Right Twin Higgs Models [38] and dark matter models with a scalar portal to the dark sector [39] (also see [40]). PNGBs have an approximate shift symmetry and their interactions are proportional to some universal symmetry breaking scale. In addition, such models have been widely discussed in literature as an effective theory for solving the g-2 anomaly (see for example [5], [6], [41]). I will focus on this form of Yukawa couplings throughout this paper.

It was shown in [5] that the discrepancy in $g_μ-2$ can be explained by a scalar with $Y_{φμμ} \sim O(10^{-3})$, while in [6], it was shown that an ALP pseudoscalar can explain both of the electron and the muon anomalies by considering the NLO contributions.

The LO and NLO contributions to the $\Delta a_{\mu,e}$ are shown in Fig. (1). The LO contribution for the scalar, as well as for the pseudoscalar, was calculated in [5]:

$$ \Delta a_l^p = -\frac{Y_{\phi ll}^2}{8\pi^2} r^{-2} \int_0^1 dz \frac{(1-z)^3}{r^{-2}(1-z)^2 + z}, \quad (5) $$

$$ \Delta a_l^s = \frac{Y_{\phi ll}^2}{8\pi^2} r^{-2} \int_0^1 dz \frac{(1+z)(1-z)^2}{r^{-2}(1-z)^2 + z}, \quad (6) $$

where $r \equiv \frac{m_\phi}{m_l}$. On the other hand, the NLO contribution includes the Barr-Zee (BZ) contribution (top right diagram in Fig. 1), the two-loop Light-By-Light (LBL) contribution (bottom left diagram in Fig. 1), and the Vacuum Polarization (VP) contribution (bottom right diagram in Fig. 1). These contributions are the same for both the scalar and the pseudoscalar cases and are given by [6]:

$$ a_{l,\phi}^{\text{BZ}} \approx \left(\frac{m_l}{4\pi^2}\right) g_{\phi\gamma\gamma} Y_{\phi ll} \ln \frac{\Lambda}{m_\phi}, \quad (7) $$

$$ a_{l,\phi}^{\text{LBL}} \approx \frac{3\alpha}{\pi} \left(\frac{m_l g_{\phi\gamma\gamma}}{4\pi}\right)^2 \ln^2 \frac{\Lambda}{m_\phi}, \quad (8) $$

$$ a_{l,\phi}^{\text{VP}} \approx \frac{\alpha}{\pi} \left(\frac{m_l g_{\phi\gamma\gamma}}{12\pi}\right)^2 \ln \frac{\Lambda}{m_\phi}, \quad (9) $$

where $\phi$ is either a scalar or a pseudoscalar, $g_{\phi\gamma\gamma}$ is the dimensionful coupling of $\phi$ to photons, and $\Lambda$ is some UV cutoff scale that is assumed to be much larger than $m_\phi$. I will set the cutoff scale $\Lambda = 1$ TeV throughout this paper. Notice that since the Lagrangian in Eq. 3 is CP-conserving, there will be no contribution to the lepton's Electric Dipole Moment (EDM). I will ignore the more general case where CP-violating terms are present.

<figure><img src="image_2.png" /><figcaption>FIG. 1. LO (top left diagram) and NLO scalar/pseudoscalar contributions to the lepton anomalous magnetic moment.</figcaption></figure>

## II. FAVORED REGION

In this section, I will investigate the parameter space and try to establish the favored region for both the scalar and the pseudoscalar cases. If we inspect equations 5 through 9, we notice the following:

* For a scalar, the LO contribution is always positive, while for a pseudoscalar it is negative. As for the NLO contributions, we can see that the LBL and the VP are always positive, while the BZ contribution depends on the sign of $g_{\phi\gamma\gamma}Y_{\phi ll}$.

* Since the central measured anomaly for the muon is positive, while for the electron it is negative, a scalar solution can easily accommodate the muon anomaly. However, yielding the central measured electron anomaly would require the assumption that $g_{\phi\gamma\gamma}Y_{\phi ll} < 0$ so that the BZ contribution can offset all other (positive) contributions. This would require large (nonperturbative) Yukawa couplings for the electron and therefore it is disfavored. Thus I will assume that all couplings are positive for the scalar case. This means that the central measured electron anomaly cannot be produced. However, it is possible to show that the contribution is within $\Delta a_e + 2\sigma_e^3$.

* As the LO contribution for the pseudoscalar case is negative, it can easily accommodate the measured electron anomaly, however, in order to yield the (positive) measured muon anomaly, one needs somewhat large couplings to photons while keeping the Yukawa coupling somewhat small in order for the NLO contributions to dominate over the LO. Nonetheless, it is possible to find such solutions while maintaining perturbative couplings as we shall see below.

* For the pseudoscalar case, if we assume that $g_{\phi\gamma\gamma} > 0$, then $Y_{\phi ee} < 0$ while $Y_{\phi\mu\mu}$ could be either positive or negative. On the other hand, assuming that

<footnote n="3">³ See added note at the end of this paper.</footnote>
$g_{\phi\gamma\gamma} < 0$ yields exactly the same solution but with opposite signs.

## A. Pseudoscalar

For concreteness, I will assume that $g_{\phi\gamma\gamma} < 0$. Fig. (2) shows the $2\sigma$ allowed regions for $e$ and $\mu$ with $g_{\phi\gamma\gamma} = -0.05 \text{ GeV}^{-1}$. Notice that there is no overlap between the two regions even for small masses, which means that there is no region in the parameter space where Eq. 4 is valid. As a matter of fact, there is no pseudoscalar solution in the whole parameter space where Eq. 4 is true. In addition, for both Yukawa couplings to have the same sign, one needs $|g_{\phi\gamma\gamma}| \gtrsim 0.03 \text{ GeV}^{-1}$, which is excluded by cosmological observations [42], therefore the pseudoscalar solution is disfavored. One can avoid these constraints by assuming that $Y_{\phi\mu\mu} > 0$ and $Y_{\phi ee} < 0$, however, one needs to justify this assumption. I will disregard the pseudoscalar solution in the remainder of this paper.

<figure><img src="image_3.png" /><figcaption>FIG. 2. The $2\sigma$ allowed region for $e$ (green) and $\mu$ (blue) of the parameter space for the pseudoscalar solution assuming $g_{\phi\gamma\gamma} = -0.05 \text{ GeV}^{-1}$. The lack of any overlap between the favored regions for the electron and the muon indicates that we cannot have a solution where Eq. 4 is valid. Therefore the pseudoscalar solution is disfavored.</figcaption></figure>

## B. Scalar

With the assumption that all couplings are positive, I attempt at finding the favored region in the parameter space for the scalar case. Here, I will only focus on solutions of the form in Eq. 4 and discard all other possibilities.

Notice that we have three parameters, namely $m_\phi$, $g_{\phi ll}$ and $g_{\phi\gamma\gamma}$. I will fix $g_{\phi\gamma\gamma}$ and keep $g_{\phi ll}$ and $m_\phi$ as free parameters. I will select appropriate benchmark points for $g_{\phi\gamma\gamma}$ by minimizing the $\chi^2$ of the electron and muon measurements:

$$ \chi^2 = \frac{(\Delta_{\mu}^{\text{Exp}} - \Delta_{\mu}^{\phi})^2}{\sigma_{\mu}^2} + \frac{(\Delta_{e}^{\text{Exp}} - \Delta_{e}^{\phi})^2}{\sigma_{e}^2}. \quad (10) $$

More concretely, I fix the value of $g_{\phi\gamma\gamma}$, then I use Eq. 10 to find $g_{\phi ll}$ and $m_\phi$ that minimize $\chi^2$, together with the value of $\chi^2$ at the minimum. Then I will scan through a wide range of $g_{\phi\gamma\gamma}$, and then set the benchmark points where $\chi^2_{min}$ is the smallest.

Fig. 3 shows $\chi^2_{min}$ for several values of $g_{\phi\gamma\gamma}$. As the plot shows, $g_{\phi\gamma\gamma} \lesssim 10^{-6} \text{ GeV}^{-1}$ yields the lowest values of $\chi^2$. Notice that $\chi^2$ becomes almost constant for smaller couplings. This is reasonable as when the coupling to photons becomes very small, the NLO contributions become negligible and the LO contribution is dominant. This high-level analysis seems to favor smaller couplings to photons, suggesting that the coupling to leptons is the dominant coupling. This is consistent with the cosmological constraints on ALPs (see for instance [42]). I will focus on this scenario and I will choose $g_{\phi\gamma\gamma} = 10^{-6}$ and $10^{-11} \text{ GeV}^{-1}$ as two benchmark points. Notice that to a good level of accuracy, the second benchmark point is representative of the entire region of $g_{\phi\gamma\gamma} \lesssim 10^{-8} \text{ GeV}^{-1}$, with very similar favored regions for the predicted mass and coupling.

<figure><img src="image_4.png" /><figcaption>FIG. 3. $\chi^2_{min}$ vs. $g_{\phi\gamma\gamma}$. The plot shows that smaller couplings to photons are favored.</figcaption></figure>

Now we can use Eqs. 6-9 in order to find the allowed region in the $m_\phi - Y_{\phi ll}$ parameter space corresponding to a $2\sigma$ deviation from the central values. In order to set an upper limit on $m_\phi$, we demand that all Yukawa couplings remain perturbative. Obviously, the most stringent bound comes from $Y_{\phi\tau\tau}$ as it has the largest value. Requiring that $Y_{\phi\tau\tau} \lesssim 1$, we obtain an upper bound on $m_\phi$ of $\sim 45(50) \text{ GeV}$ for $g_{\phi\gamma\gamma} = 10^{-11}(10^{-6} \text{ GeV})$. Fig. (4) shows the allowed region corresponding to the two benchmark points. The plots show the $2\sigma$ bands for the Yukawa couplings to electrons and muons assuming that Eq. 4 holds. In addition, the plots also show the region where the contribution to $\Delta a_e$ is within $2\sigma$ of the measured value assuming that the Yukawa couplings to leptons are independent of one another. Notice that the brown region corresponds to $\Delta a_e^\phi \in (0, \Delta a_e^{\text{Exp}} + 2\sigma_e]$ since we are assuming positive couplings. Thus as noted earlier, there is no point in the parameter space that can yield the central value of the measured $\Delta a_e = -91 \times 10^{-14}$.

We can carry the $\chi^2$ analysis further to find the favored region in the $m_\phi - g_{\phi ll}$ parameters space for the bench-
<figure><img src="image_5.png" /><figcaption>FIG. 4. The 2$\sigma$ bands corresponding to the allowed region for $e$ and $\mu$ for the scalar case at the two benchmark points with $g_{\phi\gamma\gamma} = 10^{-11}$ GeV$^{-1}$ (top) and $10^{-6}$ GeV$^{-1}$ (bottom). The brown region corresponds to $\Delta a_e^\phi \le a_e + 2\sigma_e$. Notice here that for the blue and green bands $\frac{\sqrt{Y_{\phi\mu\mu}}}{Y_{\phi ee}} = \frac{m_\mu}{m_e}$ in accordance with Eq. 4.</figcaption></figure>

mark points. Fig. (5) shows the scalar mass and coupling to leptons that minimize $\chi^2$, together with the 68% and 95% confidence level contours. Notice that this region is a subset of the allowed region in Fig. (4). For $g_{\phi\gamma\gamma} = 10^{-11}$ GeV$^{-1}$, we find a predicted scalar mass of $\sim 540$ MeV, with $g_{\phi ll} \simeq 1.45 \times 10^{-2}$ GeV$^{-1}$. For this point, one finds a predicted anomaly $\Delta a_{\mu(e)} = 273 \times 10^{-11}(8 \times 10^{-18})$. On the other hand, for the second benchmark point ($g_{\phi\gamma\gamma} = 10^{-6}$ GeV$^{-1}$), the predicted mass and coupling are 64 MeV and $5 \times 10^{-3}$ GeV$^{-1}$ respectively, with $\Delta a_{\mu(e)} = 272 \times 10^{-11}(5 \times 10^{-17})$. This is consistent with the results found in [5].

Notice that the predicted values of the electron anomaly for these benchmark points are much smaller than the (absolute) measured central value of $91 \times 10^{-14}$ (see Eq. 2), albeit they are still within the $2\sigma$ limit. This is logical as we can see from Eq. 6, $\Delta a_e \sim \Delta a_\mu \times (\frac{m_e}{m_\mu})^4$ if the Yukawa coupling has the form given in Eq. 4. This is an important prediction to test this model. That is, we claim that if indeed this model is correct, then more accurate measurements of electron anomaly should yield $\Delta a_e \sim \Delta a_\mu \times (\frac{m_e}{m_\mu})^4 \sim O(10^{-18})$. If future measurements of the electron anomaly are inconsistent with this, then the assumption in Eq. 4 would be ruled out and other explanations would be needed. Another important

prediction of this model is the tau anomaly. Given that Eq. 4 predicts the anomaly to be proportional to the mass, then $(g_\tau - 2)$ should be large enough to be measured. Specifically, we predict $(g_\tau - 2) = 8(1.5) \times 10^{-6}$ for $g_{\phi\gamma\gamma} = 10^{-11}(10^{-6})$ GeV$^{-1}$.

<figure><img src="image_6.png" /><figcaption>FIG. 5. The favored point in the $m_\phi - g_{\phi ll}$ parameter space of the scalar solution, with 68% and 95% confidence level contours for $g_{\phi\gamma\gamma} = 10^{-11}$ GeV$^{-1}$ (top) and $10^{-6}$ GeV$^{-1}$ (bottom).</figcaption></figure>

### III. EXPERIMENTAL PROBES AND LIMITS

For the mass range $\lesssim 50$ GeV, the most relevant constraints come from the BaBar experiment, beam dump experiments, the Z decay, the LEP mono-$\gamma$ searches, $ee \to \tau\tau$ and $ee \to \tau\tau + \gamma$ searches, and from the bounds on solar emission and the on emission from HB stars. Although there is overlap between the results of this chapter and [43], where they discuss a similar effective model and a UV completion through the 2HDM, there are several novel features in this work, including bounds from the Z decay, LEP searches and the solar and HB stars. In addition, here we attempt at explaining both the electron and the muon anomalies and we extend the range of the
mass over the entire range of validity. Thus our work should be viewed as complementary to theirs.

## A. BaBar

Recent results from the BaBar experiment [44] searching for the process $e^{+}e^{-} \to \mu^{+}\mu^{-}Z', Z' \to \mu^{+}\mu^{-}$ can be important for constraining the parameters space. The results can be used to extract the constraints on the process $e^{+}e^{-} \to \mu^{+}\mu^{-}Z', Z' \to \mu^{+}\mu^{-}$. The tree-level Feynman diagrams that contribute to this process are shown in Fig. (6), where I have neglected the diagrams where the scalar is radiated by the initial state particles or by the intermediate particle, since the coupling to muons dominates over the coupling to electrons and photons, and I am assuming a subleading coupling to the Z. Fig. (4) in [44] sets an upper limit $\sigma_{Max}$ on the process $e^{+}e^{-} \to \mu^{+}\mu^{-}Z', Z' \to \mu^{+}\mu^{-}$. Therefore, we can extract the limits on the $g_{\phi ll} - m_{\phi}$ parameter space by requiring:

$$ \sigma(e^{+}e^{-} \to \mu^{+}\mu^{-}\phi) \times \mathrm{Br}(\phi \to \mu^{+}\mu^{-}) < \sigma_{\mathrm{Max}}. \quad (11) $$

<figure><img src="image_7.png" /><figcaption>FIG. 6. The dominant tree-level contributions to $e^{+}e^{-} \to \mu^{+}\mu^{-}\phi$.</figcaption></figure>

The excluded part of the parameter space is shown in the light gray region in Fig. (9). Notice that it does not constrain the favored region corresponding to a $2\sigma$ deviation, although larger deviations would be constrained. This result also shows us that any new searches from BaBar (or similar experiments like Belle II) need to be at least an order of magnitude better in order to explore the favored region in the parameter space, or that searches should be made for a similar process with taus instead of muons, as the former has a larger coupling. Thus, there is a good motivation to search for processes like $ee \to 4\tau$.

Notice that in calculating the branching fraction, we are assuming that the scalar can only decay to leptons or photons. If the scalar is allowed decay to other SM particles, the limits will become even weaker. Also notice that the upper limit in BaBar's results found in Fig. (4) of [44] is spiky, therefore the upper limit used in Eq. 11 contains significant uncertainty of about a factor of a few. In extracting $\sigma_{Max}$, I was conservative and used the smallest cross-section. A less conservative estimate would relax the limits further.

## B. Z Decay<fnref n="5" />

The excellent measurements of the Z decay width and branching fractions present us with a potentially suitable tool for probing the parameter space. In particular, we can explore the limits associated with the scalar loop correction to the Z decay to a pair of leptons.

The scalar loop corrections can significantly affect the leptonic decay width of the Z boson. The NLO corrections to $Z \to l\bar{l}$ are shown in Fig. (7), where the coupling of $\phi$ to Z is assumed to be subdominant compared with the coupling to leptons. Notice here that UV divergences in the leg corrections cancel that in the vertex correction, and that for a massive $\phi$ the result is free of IR divergences.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Scalar NLO corrections to $Z \to l\bar{l}$.</figcaption></figure>

Dropping the lepton mass in the loops, and keeping $m_{\phi}$ only as an IR regulator, which is justified for $m_{\phi} \gg m_{l}$ where this bound is mostly relevant, the NLO correction is approximately given by:

$$ \delta\Gamma(Z \to l\bar{l}) \approx -\Gamma_0 \frac{m_l^2 g_{\phi ll}^2}{8\pi^2} \left[ \log\left(\frac{M_Z^2}{m_\phi^2}\right) - 2 \right], \quad (12) $$

where $\Gamma_0$ is the LO decay width given by:

$$ \Gamma_0(Z \to l\bar{l}) = \frac{g_V^2 + g_A^2}{48\pi \cos^2 \theta_W} M_Z \sqrt{1 - \frac{4m_l^2}{M_Z^2}}. \quad (13) $$

We can compare this correction with the branching fractions of the leptonic Z decays, which are given by [48]:

$$ \mathrm{Br}(Z \to e^{+}e^{-}) = (3363.2 \pm 4.2) \times 10^{-3} \%, \quad (14) $$

<footnote n="4">⁴ [45] extracts the constraints from BaBar's results for the pseudoscalar case.</footnote>

<footnote n="5">⁵ [46, 47] also discuss the 1-loop corrections to the Z-decay in the context of R-parity violating extensions to SUSY and in the context of type-II 2HDM.</footnote>
$$ \mathrm{Br}(Z \to \mu^+\mu^-) = (3366.2 \pm 6.6) \times 10^{-3} \%, \quad (15) $$

$$ \mathrm{Br}(Z \to \tau^+\tau^-) = (3369.6 \pm 8.3) \times 10^{-3} \%. \quad (16) $$

As it turns out, the decay to $\tau\tau$ provides the most stringent constraints in spite of the larger uncertainty in its branching fraction. This is due to its larger coupling given the assumption in Eq. 4.

The excluded region of the parameter space at a $2\sigma$ level is shown in red in Fig. (9). As can be seen from the plots, there is some tension between the Z-loop decay and the allowed region near $m_\phi \sim 10$ GeV, although the allowed region is not fully excluded. Notice here that since we are dropping the mass of the tau, the bound will be less reliable for the range $m_\phi < m_\tau$, however, this is unimportant as in that region the Z-decay constraint is far from the favored region and therefore our results are unaffected. On the other hand, in the region of interest where $m_\phi$ is larger than a few GeV, dropping the mass of the tau is justified and won't impact the results significantly. The same argument applies for the constraint obtain from $e^+e^- \to \tau^+\tau^-$ we discuss below.

### C. Constraints from $e^+e^- \to \tau^+\tau^-$

The analysis conducted in the previous section can be extended to the $e^+e^- \to \tau^+\tau^-$ searches conducted by the LEP experiment [49–64], the KEK collaboration [65–68], DESY-PETRA collaboration [69–74], and the SLAC-PEP experiment [75].

Here, the process of $ee \to \tau\tau$ proceeds through the s-channel with a photon or a Z propagator, and the NLO loop correction due to $\phi$ will be identical to the case of the Z decay shown in Fig. (7). So we can set the $2\sigma$ bound as:

$$ \delta\sigma_{\text{loop}}(ee \to \tau\tau) < 2\delta_{\sigma}^{\text{Measured}}, \quad (17) $$

where

$$ \sigma_{\text{loop}} = \sigma_0 \left[ 1 + \frac{m_{\tau}^2 g_{\phi ll}^2}{8\pi^2} \left( 2 - \log \left( \frac{s}{m_{\phi}^2} \right) \right) \right], \quad (18) $$

and $\sigma_0$ is the tree level cross-section. The bound is shown in Fig. (9) in light green. As the plot shows, and similar to the case of the Z decay, there is significant tension with experimental results for $m_\phi \gtrsim 10$ GeV, although this region is not fully excluded.

### D. Muon Beam Dump Experiments

Muon beam dump experiments provide a powerful tool for probing the mass range $\sim 1–200$ MeV. The relevant constraints come from Orsay [76] and the E137 experiment at SLAC [77] (also see [41] for a summary).

In the Orsay beam dump experiment, searches for the light Higgs boson in the 2HDM were conducted through looking for the process $eN \to eNH, H \to e^+e^-$. In [76], the coupling of the lighter Higgs to electrons is assumed to be:

$$ g_{\text{hee}} = \frac{m_e}{v} \tan \beta, \quad (19) $$

where $\tan \beta = v_1/v_2$, the ratio of the two doublets' VEVs. The results show the excluded region in the $\tan \beta - m_h$ parameter space. Therefore, they can readily be extrapolated to this model by setting $g_{\phi ll} \equiv \frac{\tan \beta}{v}$.

On the other hand, [77] presents the results of the E137 beam dump searches for axions produced via bremsstrahlung followed by the subsequent decay to $e^+e^-$. The results were extracted for the case of a scalar in [41], so I will just use their results.

I show these constraints in magenta (Orsay) and orange (E137) in Fig. (9). As can be seen from the plots, muon beam dump experiments exclude the region between $m_\phi \sim 1$ MeV and 30 MeV. On the other hand, the window between $m_\phi \sim 30–200$ MeV is still open. This window is projected to be explored by the proposed NA64 project at CERN [78, 79]. The NA64 experiment is a fixed-target experiment that can run in the muon mode with a beam energy of 160 GeV and is designed for searching for missing energy $\gtrsim 50$ GeV. This experiment can help probe this open window. The projected region in the parameter space is shown by the dashed line in Fig. (9).

Another proposed experiment is Fermilab's displaced decay search with a muon beam energy of 3 GeV [80]. However, the projected sensitivity of this experiment covers only a part of the projected sensitivity of the NA64 experiment, therefore I will not plot it here.

### E. Constraints from LEP Mono-$\gamma$ and $ee \to \tau\tau\gamma$ Searches

The LEP mono-$\gamma$ searches were conducted to set limits on the number of neutrinos via studying the process $e^-e^- \to \nu\bar{\nu}\gamma$. The results can be used to set limits on the $m_\phi - g_{\phi ll}$ parameter space by considering the process $e^+e^- \to \gamma\phi$ for the range $m_\phi \le 2m_e$. The tree-level process proceeds through the t and u channels, however, due to the smallness of the scalar's coupling to the electron, the constraints are weak. On the other hand, stronger constraints can be obtained through the triangle diagram shown in Fig. (8) where the $\tau$ runs in the loop.

Similar to the case of the Higgs, we can write the effective Lagrangian as:

$$ \mathcal{L}_{\text{eff}} = c_{\gamma} \frac{\alpha}{\pi v_{\phi}} A_{\mu\nu} A^{\mu\nu} \phi + c_{\gamma Z} \frac{\alpha}{\pi v_{\phi}} A_{\mu\nu} Z^{\mu\nu} \phi, \quad (20) $$

where $v_\phi = g_{\phi\tau\tau}^{-1}$. The effective couplings can be extracted from the results of the Higgs. For example, we
<figure><img src="image_9.png" /><figcaption>FIG. 8. $e^{+}e^{-} \rightarrow \gamma\phi$ triangle diagram.</figcaption></figure>

can use the result of the Higgs decay to $\gamma\gamma$ where instead of the top running in the loop, we have the tau. This gives:

$$c_{\gamma} = \frac{1}{8} |A_{\tau}^{\phi}(\tau_{\tau})|, \qquad (21)$$

where $\tau_{\tau} = \frac{4m_{\tau}^{2}}{m_{\phi}^{2}}$, and:

$$A_{\tau}^{\phi}(\tau) = 2\tau \left[ 1 + (1 - \tau)f(\tau) \right], \qquad (22)$$

$$f(\tau) = \begin{cases} \left[ \sin^{-1}\left(\frac{1}{\sqrt{\tau}}\right) \right]^2, & \tau \ge 1 \\ -\frac{1}{4} \left[ \ln\left(\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}\right) - i\pi \right]^2, & \tau < 1. \end{cases} \qquad (23)$$

Similarly, we find $c_{\gamma Z}$ from the Higgs decay to $Z\gamma$. We obtain:

$$c_{\gamma Z} = \frac{1}{4 \sin \theta_W} |B_{\tau}^{\phi}(\tau_{\tau}, \lambda_{\tau})|, \qquad (24)$$

where $\lambda_{\tau} = \frac{4m_{\tau}^{2}}{m_{Z}^{2}}$ and:

$$B_{\tau}^{\phi}(\tau, \lambda) = 2N_{c}^{\tau} \left[ \frac{Q_{\tau}(I_{3}^{\tau} - 2Q_{\tau}s_{w}^{2})}{c_{w}} \right] [I_{1}(\tau, \lambda) - I_{2}(\tau, \lambda)], \qquad (25)$$

$$I_{1}(\tau, \lambda) = \frac{\tau\lambda}{2(\tau - \lambda)} + \frac{\tau^{2}\lambda^{2}}{2(\tau - \lambda)^{2}}[f(\tau) - f(\lambda)] + \frac{\tau^{2}\lambda}{(\tau - \lambda)^{2}}[g(\tau) - g(\lambda)], \qquad (26)$$

$$I_2(\tau, \lambda) = -\frac{\tau\lambda}{2(\tau - \lambda)}[f(\tau) - f(\lambda)], \qquad (27)$$

$$g(\tau) = \begin{cases} \sqrt{\tau - 1} \sin^{-1}\left(\frac{1}{\sqrt{\tau}}\right), & \tau \ge 1 \\ \frac{\sqrt{\tau - 1}}{2} \left[ \ln\left(\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}\right) - i\pi \right], & \tau < 1. \end{cases} \qquad (28)$$

Here $Q_{\tau} = -1$, $N_{c}^{\tau} = 1$ and $I_{3}^{\tau} = -\frac{1}{2}$. Armed with this, we can find the cross-section of the triangle diagram in Fig. (8):

$$\sigma(e^{+}e^{-} \rightarrow \gamma\phi) = \left( \frac{e^2 g_{\phi l l}}{32\pi^2 \sqrt{2\pi} c_w} \right)^2 \left( \frac{1 - m_{\phi}^2/s}{1 - m_Z^2/s} \right)^2 \times \left[ (1 - 4s_w^2 + 8s_w^4)g^2 c_{\gamma Z}^2 + 32c_w^2 e^2 c_{\gamma}^2 (1 - m_Z^2/s)^2 - 8c_w g e c_{\gamma Z} (1 - 4s_w^2)(1 - m_Z^2/s) \right], \qquad (29)$$

where $s_w, c_w$ are the sine and cosine of the Weinberg angle respectively. The results from ALEPH [81, 82], L3 [83–85], OPAL [86] and DELPHI collaborations [87] can be used in order set constraints on the scalar by requiring that at a $2\sigma$ level:

$$\sigma(e^{+}e^{-} \rightarrow \gamma\phi) < 2\delta_{\sigma}^{\text{Measured}}, \qquad (30)$$

for $m_{\phi} \le 2m_e$. Notice that for $m_{\phi} \le 2m_e$, $\phi$ could decay to $\gamma\gamma$, however, for both benchmark points, the decay length is orders of magnitude larger than the dimensions of the detector, so $\phi$ appears as missing energy. On the other hand, for $m_{\phi} > 2m_l$, where $l = e, \mu, \tau, \phi$ can decay to a lepton pair, so we can use LEP $ee \rightarrow ll + \gamma$ searches for that region. The cross-section for $ee \rightarrow ll\gamma$ can be readily obtained by multiplying the cross-section in Eq. 29 by the appropriate branching fraction for the mass range. The constraints are shown in brown in Fig. (9).

## F. Constraints from Solar Emission, HB Stars and SN1987A

Supernova 1987 (SN1987A), HB stars and solar emission constraints can impose stringent constraints on ALPs for masses $\lesssim 1$ GeV (see [42, 88] for instance). However, such constraints are only relevant if the ALP’s dominant coupling is to photons. Since the coupling to photons in the type of models we are considering in this paper is favored to be much less than that to leptons (at least for the scalar case), those limits need to be revisited.

If $\phi$ is produced in the sun or in HB stars, then it could affect the measured energy loss rate of the star when it streams out, in addition to affecting the star’s evolution. Therefore one can obtain a bound on $m_{\phi}$ and $g_{\phi l l}$ by requiring that the amount of energy carried away by the scalar be less than the observed limits.

In this model where the scalar’s couplings to leptons are dominant, $\phi$ will be mainly produced in stars via its interaction with electrons. Assuming that the electrons are nonrelativistic and nondegenerate, which is a good assumption for the solar and HB stars’ mediums, the Compton-like scattering $e^{-}\gamma \rightarrow e^{-}\phi$ dominates over both Bremsstrahlung $e^{+}e^{-} \rightarrow e^{+}e^{-}\phi$ and electron-positron annihilation. For $m_{\phi} \ll m_e$, we can neglect the recoil energy of the electron and to a good approximation the energy loss rate per unit volume is given by [89]:

$$Q = \frac{n_e}{\pi^2} \int_{m_{\phi}}^{\infty} d\omega \frac{\sigma(\omega)\omega^3}{e^{\omega/T} - 1}, \qquad (31)$$

where $\omega$ is the energy of the photon and $n_e$ is the number density of the electrons, which in terms of the electron fraction in the star $y_e$, the mass density of the star $\rho$, and the atomic mass unit $m_u$, is given by:

$$n_e = \frac{y_e \rho}{m_u}. \qquad (32)$$
This calculation was done by Grifols & Massó in [90] with $m_\phi$ neglected. Here, I will keep the mass of the scalar but assume it's less than $m_e$, and keep the assumption that the recoil energy of the electron is small. The cross-section is given by:

$$ \sigma(e^{-}\gamma \rightarrow e^{-}\phi) = \frac{e^{2}g_{\phi ll}^{2}m_{e}^{2}}{32\pi\omega^{2}}F(x, y), \quad (33) $$

where $x = \frac{m_\phi}{\omega}$, $y = \frac{m_e}{\omega}$ and $F(x, y)$ is some complicated function that I relegate to Appendix A. The bound on the solar (HB) emission rate is [89]

$$ \dot{\epsilon} \le 2(10) \text{ erg g}^{-1} \text{ s}^{-1}. \quad (34) $$

Given this bound, we can use Eq. 33 to solve Eq. 31 numerically in order to find the excluded region in the parameter space. Here I set the average temperature to be $10^7(10^8)$ K for the solar (HB) medium, the average density to be $10^2(10^4)$ g cm$^{-2}$ for the solar (HB) medium and use $y_e = 0.5$ for both. The yellow region in Fig. (9) shows the excluded part of the parameter space by the solar constraints, while the dark green region shows the excluded part by HB stars. As the plot shows, all masses below $\sim 350$ KeV are excluded. However, we must note that for $m_\phi$ close to the electron mass, the bound is less rigorous as our assumption of a small electron's recoil energy becomes less valid. However, I checked numerically that including the recoil energy of the final state enellectron does not significantly impact the bound. Also, notice that the HB constraints are much more stringent than the solar ones in spite of the weaker bound in Eq. 34. This is due to the higher temperature of HB stars compared with the sun, which yields a larger Boltzmann factor in Eq. 31.

Of course, this analysis is valid only if $\phi$ streams freely out of the star, i.e. if it does not get trapped inside the medium of the star. In order to verify the validity of this assumption, we can calculate the mean free path of $\phi$ and compare it with the radius of the star. The mean free path is given by:

$$ \lambda \sim \frac{1}{\sigma n_e} = \left[ \frac{e^2 g_{\phi ll}^2}{32 \pi m_u} m_e^2 \omega^2 y_e \rho F(x, y) \right]^{-1}. \quad (35) $$

Assuming that the average photons energy is given by:

$$ \langle \omega \rangle = \frac{\pi^4}{30 \zeta(3)} T \approx 2.701 T, \quad (36) $$

one can easily check that for the entire mass range of interest, the mean free path is orders of magnitude larger than the radius of the sun or the typical radius of an HB star, thereby justifying the free streaming assumption.

A similar argument applies for scalars produced in supernovas. However, it was shown in [40] that for a scalar produced in the supernova core, the mean free path is given by:

$$ \lambda_{\text{mpf}} \sim 10 \text{ m} \left( \frac{g_{\phi ll}^{-1}}{10^6 \text{ GeV}} \right)^2, \quad (37) $$

which means that for $g_{\phi ll}$ larger than $10^{-6}$ GeV$^{-1}$, the scalar gets trapped in the core and never streams out. Thus, SN1987A does not costrain the favored region and therefore we ignore it.

<figure><img src="image_10.png" /><figcaption>FIG. 9. The exclusion plot of the $m_\phi - g_{\phi ll}$ parameter space corresponding to the $g_{\phi\gamma\gamma} = 10^{-11}$ GeV$^{-1}$ (top) and $g_{\phi\gamma\gamma} = 10^{-6}$ GeV$^{-1}$ (bottom) benchmark points. The plot shows the $2\sigma$ allowed region for $\Delta a_{\mu(e)}$ (beige); the excluded region by the NLO correction to the leptonic Z decay (red), the excluded region due to the NLO corrections to $ee \rightarrow \tau\tau$ (light green), Orsay (magenta), E137 (orange), LEP mono-$\gamma$ and $ee \rightarrow ll + \gamma$ searches (brown), solar emission (yellow) and HB star emission (dark green); and the projected region for NA64 (dashed). The plots also show the favored points that correspond to the minimum $\chi^2$ in the two scenarios.</figcaption></figure>

## G. Discussion

We have shown that the constraints from BaBar, beam dump experiments, the NLO correction to the Z decay, LEP mono-$\gamma$ searches, $ee \rightarrow \tau\tau(\gamma)$ searches, and solar and HB star constraint, exclude a significant part of the parameter space. Fig. (9) shows that all masses below $\sim 350$ KeV are excluded by solar and HB constraints; beam dump experiments exclude the region above 1 MeV up to $\sim 30$ MeV, while masses above $\sim 10$ GeV are in significant tension with the NLO corrections to the Z decay to $\tau\tau$ and the NLO corrections to $ee \rightarrow \tau\tau$, although this region is not entirely excluded. This leaves a large open window between $\sim 30$ MeV and $\sim 10$ GeV that is
most favored to be explored.

Part of this region is projected to be explored by the proposed NA64 experiment, which is projected to cover the mass range from 1 MeV and up to the di-muon mass. The first benchmark point corresponding to $g_{\phi\gamma\gamma} = 10^{-6}$ GeV$^{-1}$ lies in this region. On the other hand, the second benchmark point lies outside this region but could be explored by the Belle II experiment [91]. The Belle II experiment has recently started collecting data and it is expected that by 2025, it would have reached a total integrated luminosity of 50 ab$^{-1}$. This projects it to be more sensitive than the BaBar experiment, and therefore might help explore more of the parameter space above the dimuon threshold which is not covered by NA64 and up to ~5 GeV. This leaves two windows, one above ~5 and the other is between ~350 KeV and 1 MeV.

The International Linear Collider (ILC) [92] (if built) might help explore both regions. For the former, the ILC's ultra-precision measurement of the Z production and subsequent decay can help improve the bounds from the Z decay by lowering the measured uncertainties, while for the latter region, the $e^{+}e^{-} \to h + \mathcal{E}_{T}$ channel can help investigate the hypothetical process $e^{+}e^{-} \to h + \phi$.

A final point to mention is that due to the assumption made in Eq. 4, the dominant coupling will be to the tau lepton, and therefore collider and beam dump experiments limits are more stringent. If we assume a different type of coupling to leptons, say by assuming a suppressed coupling to the tau compared to the muon, then the constraints will be alleviated, and more of the parameter space will open. We will very briefly discuss one model where this can be achieved in the next section.

## IV. A RADION SOLUTION FOR THE $g-2$ ANOMALY?

In the Randall-Sundrum (RS) model [93], the radion is the scalar field that parametrizes the fluctuations of the extra dimension around its potential minimum. The radion could pose an interesting possibility for solving the $g-2$ anomaly due to its unique couplings to matter. More specifically, the radion's coupling to matter is highly model-dependent and varies according to the localization of the matter fields on either of the branes or in the bulk. The coupling to brane-localized matter is given by [94]:

$$ \frac{\phi(x)}{\Lambda_{UV,IR}} \mathrm{Tr} T_{\mu\nu}, \qquad (38) $$

where $\phi(x)$ is the 4D radion field, $\Lambda_{UV,IR}$ is the radion constant on the UV and IR branes respectively and $T_{\mu\nu}$ is the stress-energy tensor. Since the UV scale is typically many orders of magnitude larger than the IR scale, it is possible to suppress the coupling to the tau lepton compared to the muon by assuming that the former

is localized on the UV brane, while assuming that the latter is localized on the IR brane. This way, one could alleviate all of the constraints (except for beam dump experiments) as the couplings will be rescaled by $m_\mu/m_\tau$.

Although the typical mass of the radion is comparable to the Electroweak (EW) scale (~100 GeV), much lighter masses can be achieved through the Contino-Pomarol-Rattzzi (CPR) mechanism [95] as was demonstrated in [96, 97].

Another interesting aspect of a radion solution is that the radion could couple to nucleons and pions through quarks and gluons [96], which presents additional experimental probes. Focusing on the coupling to pions, we can write the effective Lagrangian as:

$$ \mathcal{L}_{\phi\pi\pi} = g_{\phi\pi\pi} m_{\pi}^2 \phi\pi^+\pi^-, \qquad (39) $$

where $g_{\phi\pi\pi}$ is the effective radion coupling to pions with dimension $(mass)^{-1}$. If the radion is heavy enough, it could decay to $\pi^+\pi^-$:

$$ \Gamma(\phi \to \pi^+\pi^-) = \frac{g_{\phi\pi\pi}^2 m_{\pi}^4}{16\pi m_{\phi}^2} \left(1 - \frac{4m_{\pi}^2}{m_{\phi}^2}\right)^{\frac{1}{2}}. \qquad (40) $$

The decay width in Eq. 40 could be small for typical values of $g_{\phi\pi\pi}$ but could still be measurable. For instance, the decay width for a 400 MeV radion with $g_{\phi\pi\pi} = 0.3$ GeV$^{-1}$ would be $\simeq 1.2$ KeV.

For lighter masses, searches for the rare pion decay $\pi^- \to \mu^- \bar{\nu}_\mu \phi$ could provide an interesting search option. If we assume that $g_{\phi\pi\pi}$ dominates over $g_{\phi\bar{\nu}}$, then the branching fraction of this hypothetical decay would be given by:

$$ \mathrm{Br}(\pi \to \mu\nu\phi) \simeq 1.35 \times 10^{-2} \left( \frac{g_{\phi\pi\pi}}{\mathrm{GeV}^{-1}} \right)^2 \%. \qquad (41) $$

For instance, $g_{\phi\pi\pi} \sim 0.03$ GeV$^{-1}$ would yield a branching fraction comparable to the observed rare decay $\pi^+ \to e^+e^-e^+\nu_e$.

## V. CONCLUSIONS

The $g-2$ anomaly remains one of the best ways to search for physics BSM. In this paper, I investigated a class of models with a scalar/pseudoscalar that has a coupling to leptons proportional to the lepton's mass.

We saw in this paper that for the case of a pseudoscalar solution, there is no region in the parameter space that could simultaneously solve both the electron and the muon anomalies with Yukawa couplings of the form in Eq. 4. However, if no assumption is made regarding the form of the Yukawas, it is possible to have a pseudoscalar solution for both of the anomalies. Nevertheless, this solution is not very attractive since it would require tuning
the scalar's coupling to photons to be somewhat large in order for both of the Yukawa couplings to have the same sign. Such a large coupling to photons would be disfavored by cosmological observations.

On the other hand, a scalar can simultaneously provide a solution for both anomalies while having the required form of the Yukawa couplings. In such a case, we demonstrated that such a solution favors smaller couplings to photons ~$10^{-6}$ GeV$^{-1}$, and we established the corresponding favored region in the $m_\phi - g_{\phi ll}$ parameter space that corresponds to two representative benchmark points. We investigated the experimental constraints from the BaBar experiment, beam dump experiments, the NLO corrections to the Z decay and to ee → ττ, LEP mono-γ and ee → ll + γ searches, and from solar and HB emission bounds, and we saw that a significant part of the parameter space is excluded. In particular, we found only two open windows for $m_\phi$ between ~350 KeV and 1 MeV, and between ~30 MeV and ~50 GeV, with the region above ~10 GeV being in significant tension with experiment.

The tools used in this paper can be used to constrain other solutions to the g-2 anomaly, such as solutions that adopt the Z' or the dark photon to explain it. I expect that the limits on these solutions would not be too different from the scalar case for the same range of masses and couplings, they are nonetheless worthwhile investigating.

Future experiments, such as NA64, Belle II, and the ILC can help explore significant regions of the parameter space, and one hopes that in the near future, enough data would be collected to shed more light on the remaining open windows, and thus help explore the viability of this solution.

## ACKNOWLEDGMENTS

I would like to thank John Terning for his valuable insight. I would also like to thank John Conway, Lloyd

Knox and Max Chertok for answering my questions.

## Appendix A: Explicit Expression for $F(x, y)$

$$ F(x, y) = \frac{2y}{(1 + 2y + 2y^2)^4 (1 + 4y^4 + 4x^2y^2)} \times \left[ (1 + 2y + 2y^2)(1 + 2y + 2y^2 + 4y^3)^2 (1 + 4y^4 + 4x^2y^2) \times \log \left( \frac{1 + 2y^2 + 2y\sqrt{1-x^2}}{1 + 2y^2 - 2y\sqrt{1-x^2}} \right) - 4y\sqrt{1-x^2} (16y^8(2x^2+9) + 16y^7(4x^2+7) + 16y^6(5x^2+4) + 40y^5(2x^2+1) + 24y^4(2x^2+1) + 4y^3(6x^2+5) + 4y^2(x^2+3) + 64y^{10} + 128y^9 + 6y+1) \right] \quad (A1) $$

**Added note:** After finalizing this paper, a new results [98] was published claiming a 2.4σ discrepancy in the measured electron anomaly:

$$ \Delta a_e = -88 \pm 36 \times 10^{-14}. \quad (A2) $$

If this yet uncorroborated result is true, then our results cannot explain both anomalies simultaneously at the 2σ level, although the results remain valid at the 3σ level. This is because the electron anomaly will be negative at the 2σ level and a novel explanation would be needed to explain the sign discrepancy between the two anomalies. In this case, the electron anomaly in this work would serve more as a constraint on the allowed region of muon anomaly, and the result in Figure 9 would be relevant for the muon anomaly only. Either way, a corroborating result is still needed to confirm the value and the sign of the electron anomaly.

[1] P. J. Mohr, B. N. Taylor and D. B. Newell, Rev. Mod. Phys. 84, 1527 (2012), arXiv: 1203.5425

[2] G. W. Bennett et al. [Muon g-2 Collaboration], Phys. Rev. Lett. 89, 101804 (2002), arXiv: hep-ex/020800.

[3] G. W. Bennett et al. [Muon g-2 Collaboration], Phys. Rev. Lett. 92, 161802 (2004), arXiv: hep-ex/040100.

[4] G. W. Bennett et al. [Muon g-2 Collaboration], Phys. Rev. D 73, 072003 (2006), arXiv: hep-ex/060203.

[5] C. Y. Chen, H. Davoudiasl, W. J. Marciano and C. Zhang, Phys. Rev. D 93, no. 3, 035006 (2016), arXiv: 1511.04715 [hep-ph].

[6] W. J. Marciano, A. Masiero, P. Paradisi and M. Passera, Phys. Rev. D 94, no. 11, 115033 (2016), arXiv: 1607.01022 [hep-ph].

[7] P. deNiverville, H. S. Lee and M. S. Seo, arXiv: 1806.00757 [hep-ph].

[8] T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. D 91, no. 3, 033006 (2015), arXiv: 1412.8284 [hep-ph].

[9] J. Grange et al. [Muon g-2 Collaboration], arXiv: 1501.06858

[10] G. Venanzoni [Fermilab E989 Collaboration], Nucl. Part. Phys. Proc. 273-275, 584 (2016), arXiv: 1411.2555.

[11] N. Saito [J-PARC g-2/EDM Collaboration], AIP Conf. Proc. 1467, 45 (2012).

[12] T. Mibe [J-PARC g-2 Collaboration], Nucl. Phys. Proc. Suppl. 218, 242 (2011).

[13] M. Lindner, M. Platscher and F. S. Queiroz, "A Call for New Physics: The Muon Anomalous Magnetic Moment and Lepton Flavor Violation," Phys. Rept. 731, 1 (2018) [arXiv:1610.06587 [hep-ph]].

[14] D. Stockinger, J. Phys. G 34, R45 (2007), arXiv: hep-
ph/060916.

[15] P. Langacker, Rev. Mod. Phys. **81**, 1199 (2009), arXiv: 0801.1345 [hep-ph].

[16] S. Baek, N. G. Deshpande, X. G. He and P. Ko, Phys. Rev. D **64**, 055006 (2001), arXiv: hep-ph/010414.

[17] E. Ma, D. P. Roy and S. Roy, Phys. Lett. B **525**, 101 (2002), arXiv: hep-ph/011014.

[18] S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B **513**, 119 (2001), arXiv: hep-ph/010222.

[19] M. Pospelov, Phys. Rev. D **80**, 095002 (2009), arXiv: 0811.1030 [hep-ph].

[20] J. Heeck and W. Rodejohann, Phys. Rev. D **84**, 075007 (2011), arXiv: 1107.5238 [hep-ph].

[21] K. Harigaya, T. Igari, M. M. Nojiri, M. Takeuchi and K. Tobe, JHEP **1403**, 105 (2014), arXiv: 1311.0870 [hep-ph].

[22] W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. Lett. **113**, 091801 (2014), arXiv: 1406.2332 [hep-ph].

[23] W. Altmannshofer, C. Y. Chen, P. S. Bhupal Dev and A. Soni, “Lepton flavor violating Z explanation of the muon anomalous magnetic moment,” Phys. Lett. B **762**, 389 (2016), arXiv: 1607.06832 [hep-ph].

[24] N. V. Krasnikov, arXiv: 1702.04596 [hep-ph].

[25] E. O. Iltan and H. Sundu, arXiv: hep-ph/010310.

[26] Y. Omura, E. Senaha and K. Tobe, JHEP **1505**, 028 (2015) [arXiv: 1502.07824 [hep-ph]].

[27] A. Broggio, E. J. Chun, M. Passera, K. M. Patel and S. K. Vempati, JHEP **1411**, 058 (2014), arXiv: 1409.3199 [hep-ph].

[28] L. Wang and X. F. Han, JHEP **1505**, 039 (2015), arXiv: 1412.4874 [hep-ph].

[29] T. Abe, R. Sato and K. Yagyu, JHEP **1507**, 064 (2015), arXiv: 1504.07059 [hep-ph].

[30] J. Cao, P. Wan, L. Wu and J. M. Yang, “Lepton-Specific Two-Higgs Doublet Model: Experimental Constraints and Implication on Higgs Phenomenology,” Phys. Rev. D **80**, 071701 (2009) [arXiv:0909.5148 [hep-ph]].

[31] B. Batell, M. Pospelov and A. Ritz, Phys. Rev. D **83**, 054005 (2011) [arXiv:0911.4938 [hep-ph]].

[32] S. Dutta, A. Goyal and M. P. Singh, JHEP **07**, 076 (2019) [arXiv:1809.07877 [hep-ph]].

[33] A. Freitas, J. Lykken, S. Kell and S. Westhoff, JHEP **1405**, 145 (2014), arXiv: 1402.7065 [hep-ph].

[34] D. Chakraverty, D. Choudhury and A. Datta, Phys. Lett. B **506**, 103 (2001), arXiv: hep-ph/010218.

[35] K. m. Cheung, Phys. Rev. D **64**, 033001 (2001), arXiv: hep-ph/010223.

[36] H. Davoudiasl, H. S. Lee and W. J. Marciano, Phys. Rev. D **89**, no. 9, 095006 (2014), arXiv: 1402.3620 [hep-ph].

[37] C. Hagmann, H. Murayama, G. G. Raffelt, L. J. Rosenberg, and K. V. Bibber (2008).

[38] G. L. Liu and Q. G. Zeng, Eur. Phys. J. C **79**, no.7, 612 (2019) [arXiv:1811.04777 [hep-ph]].

[39] Y. Nomura and J. Thaler, Phys. Rev. D **79**, 075008 (2009) [arXiv:0810.5397 [hep-ph]].

[40] R. Essig, R. Harnik, J. Kaplan and N. Toro, Phys. Rev. D **82**, 113008 (2010) [arXiv:1008.0636 [hep-ph]].

[41] C. Y. Chen, M. Pospelov and Y. M. Zhong, Phys. Rev. D **95**, no. 11, 115005 (2017), arXiv: 1701.07437 [hep-ph].

[42] J. Jaeckel and M. Spannowsky, Phys. Lett. B **753**, 482 (2016), arXiv: 1509.00476 [hep-ph].

[43] B. Batell, N. Lange, D. McKeen, M. Pospelov and A. Ritz, “Muon anomalous magnetic moment through

the leptonic Higgs portal,” Phys. Rev. D **95**, no. 7, 075003 (2017) [arXiv:1606.04943 [hep-ph]].

[44] J. P. Lees et al. [BaBar Collaboration], Phys. Rev. D **94**, no. 1, 011102 (2016), arXiv: 1606.03501 [hep-ex].

[45] M. Bauer, M. Neubert and A. Thamm, “Collider Probes of Axion-Like Particles,” JHEP **1712**, 044 (2017) [arXiv: 1708.00443 [hep-ph]].

[46] O. Lebedev, W. Loinaz and T. Takeuchi, “Constraints on R-parity violating couplings from lepton universality,” Phys. Rev. D **61**, 115005 (2000)

[47] O. Lebedev, W. Loinaz and T. Takeuchi, “Constraints on two Higgs doublet models at large tan beta from W and Z decays,” Phys. Rev. D **62**, 055014 (2000)

[48] C. Patrignani et al. [Particle Data Group], Chin. Phys. C **40**, no. 10, 100001 (2016).

[49] G. Abbiendi et al. [OPAL Collaboration], Eur. Phys. J. C **33**, 173 (2004)

[50] G. Abbiendi et al. [OPAL Collaboration], Eur. Phys. J. C **19**, 587 (2001)

[51] G. Abbiendi et al. [OPAL Collaboration], Eur. Phys. J. C **13**, 553 (2000)

[52] M. Acciarri et al. [L3 Collaboration], Z. Phys. C **62**, 551 (1994).

[53] M. Z. Akrawy et al. [OPAL Collaboration], Phys. Lett. B **247**, 458 (1990).

[54] P. D. Acton et al. [OPAL Collaboration], Z. Phys. C **58**, 219 (1993).

[55] M. Z. Akrawy et al. [OPAL Collaboration], Phys. Lett. B **235**, 379 (1990).

[56] P. Abreu et al. [DELPHI Collaboration], Eur. Phys. J. C **16**, 371 (2000).

[57] M. Z. Akrawy et al. [OPAL Collaboration], Phys. Lett. B **240**, 497 (1990).

[58] P. Aarnio et al. [DELPHI Collaboration], Phys. Lett. B **241**, 425 (1990).

[59] P. Abreu et al. [DELPHI Collaboration], Nucl. Phys. B **367**, 511 (1991).

[60] D. Buskulic et al. [ALEPH Collaboration], Z. Phys. C **60**, 71 (1993).

[61] D. Decamp et al. [ALEPH Collaboration], Z. Phys. C **48**, 365 (1990).

[62] O. Adriani et al. [L3 Collaboration], Phys. Lett. B **306**, 187 (1993).

[63] D. Decamp et al. [ALEPH Collaboration], Z. Phys. C **53**, 1 (1992).

[64] D. Buskulic et al. [ALEPH Collaboration], Phys. Lett. B **378**, 373 (1996).

[65] I. Adachi et al. [TOPAZ Collaboration], Phys. Lett. B **208**, 319 (1988).

[66] H. R. Band et al. [MAC Collaboration], Phys. Lett. B **218**, 369 (1989).

[67] C. Velissaris et al. [AMY Collaboration], Phys. Lett. B **331**, 227 (1994).

[68] B. Howell et al. [TOPAZ Collaboration], Phys. Lett. B **291**, 206 (1992).

[69] H. J. Behrend et al. [CELLO Collaboration], Phys. Lett. B **114B**, 282 (1982).

[70] W. Braunschweig et al. [TASSO Collaboration], Z. Phys. C **43**, 549 (1989).

[71] W. Bartel et al. [JADE Collaboration], Phys. Lett. B **161B**, 188 (1985).

[72] S. Hegner et al. [JADE Collaboration], Z. Phys. C **46**, 547 (1990).

[73] H. J. Behrend et al. [CELLO Collaboration], Phys. Lett.
B **222**, 163 (1989).

[74] W. Bartel et al. [JADE Collaboration], Z. Phys. C **30**, 371 (1986).

[75] K. K. Gan et al., Phys. Lett. **153B**, 116 (1985).

[76] M. Davier and H. Nguyen Ngoc, Phys. Lett. B **229**, 150 (1989).

[77] J. D. Bjorken et al., Phys. Rev. D **38**, 3375 (1988).

[78] S. Andreas et al., arXiv: 1312.3309 [hep-ex].

[79] S. N. Gninenko, N. V. Krasnikov and V. A. Matveev, Phys. Rev. D **91**, 095015 (2015), arXiv: 1412.1400 [hep-ph].

[80] A. Chapelain [Muon g-2 Collaboration], EPJ Web Conf. **137**, 08001 (2017), arXiv: 1701.02807.

[81] A. Heister et al. [ALEPH Collaboration], “Single photon and multiphoton production in $e^{+}e^{-}$ collisions at $\sqrt{s}$ up to 209-GeV,” Eur. Phys. J. C **28**, 1 (2003).

[82] R. Barate et al. [ALEPH Collaboration], “Single photon and multiphoton production in $e^{+}e^{-}$ collisions at a center-of-mass energy of 183-GeV,” Phys. Lett. B **429**, 201 (1998).

[83] M. Acciarri et al. [L3 Collaboration], “Determination of the number of light neutrino species from single photon production at LEP,” Phys. Lett. B **431**, 199 (1998).

[84] M. Acciarri et al. [L3 Collaboration], “Single and multiphoton events with missing energy in $e^{+}e^{-}$ collisions at 161-GeV < $\sqrt{s}$ < 172-GeV,” Phys. Lett. B **415**, 299 (1997).

[85] P. Achard et al. [L3 Collaboration], “Single photon and multiphoton events with missing energy in $e^{+}e^{-}$ collisions at LEP,” Phys. Lett. B **587**, 16 (2004).

[86] K. Ackerstaff et al. [OPAL Collaboration], “Multiphoton production in $e^{+}e^{-}$ collisions at $\sqrt{s}$ = 183-GeV,” Phys. Lett. B **438**, 379 (1998).

[87] J. Abdallah et al. [DELPHI Collaboration], “Photon events with missing energy in $e^{+}e^{-}$ collisions at $\sqrt{s}$ = 130-GeV to 209-GeV,” Eur. Phys. J. C **38**, 395 (2005).

[88] D. Cadamuro and J. Redondo, JCAP **1202**, 032 (2012), arXiv: 1110.2895 [hep-ph].

[89] G. Raffelt,

[90] J. A. Grifols and E. Masso, “Constraints on Finite Range Baryonic and Leptonic Forces From Stellar Evolution,” Phys. Lett. B **173**, 237 (1986).

[91] T. Abe et al. [Belle-II Collaboration], “Belle II Technical Design Report,” arXiv: 1011.0352.

[92] T. Behnke et al., “The International Linear Collider Technical Design Report - Volume 1: Executive Summary,” arXiv: 1306.6327.

[93] L. Randall and R. Sundrum, Phys. Rev. Lett. **83** (1999).

[94] C. Csaki, M. L. Graesser and G. D. Kribs, Phys. Rev. D **63**, 065002 (2001).

[95] R. Contino, A. Pomarol and R. Rattazzi, talk by R. Rattazzi at Planck 2010, CERN [slides]; talk by A. Pomarol, 2010 Madrid Christmas Workshop.

[96] F. Abu-Ajamieh, J. S. Lee and J. Terning, JHEP **1810**, 050 (2018), arXiv: 1711.02697 [hep-ph].

[97] F. Abu-Ajamieh, Int. J. Mod. Phys. A **33**, no. 24, 1850144 (2018), arXiv: 1803.01249 [hep-ph].

[98] T. Aoyama, T. Kinoshita and M. Nio, Atoms **7**, no.1, 28 (2019).