Systematically Testing Singlet Models for $(g-2)_\mu$

We comprehensively study all viable new-physics scenarios that resolve the muon $(g-2)_\mu$ anomaly with only Standard Model singlet particles coupled to muons via renormalizable interactions. Since such models are only viable in the MeV -- TeV mass range and require sizable muon couplings, they predict abundant accelerator production through the same interaction that resolves the anomaly. We find that a combination of fixed-target (NA64$\mu$, $M^3$), $B$-factory (BABAR, Belle II), and collider (LHC, muon collider) searches can cover nearly all viable singlets scenarios, independently of their decay modes. In particular, future muon collider searches offer the only certain test of singlets above the GeV scale, covering all higher masses up to the TeV-scale unitarity limit for these models. Intriguingly, we find that $\mathcal{O}(100~\mathrm{GeV})$ muon colliders may yield better coverage for GeV-scale singlets compared to TeV-scale concepts, which has important implications for the starting center-of-mass energy of a staged muon collider program.


Introduction
The Fermilab Muon g − 2 collaboration has recently released its first measurement of the muon's anomalous magnetic moment [1][2][3][4]. Their new results are consistent with previous Brookhaven measurements [5], so the world average for a µ ≡ 1 2 (g − 2) µ now deviates from its Standard Model (SM) phenomenological prediction  by ∆a µ = (251 ± 59) × 10 −11 , (1.1) which yields a statistically significant 4.2 σ discrepancy and may be the first laboratory evidence of physics beyond the SM (BSM). 1 If the anomaly is due to BSM states, their presence induces the effective operator where µ L and µ c are the muon Weyl spinors, v = 246 GeV is the Higgs vacuum expectation value (VEV), C eff is a dimensionless coefficient and Λ is a BSM mass scale. This effective operator generates a contribution a BSM µ to (g −2) µ that can in principle resolve the anomaly by setting a BSM µ = ∆a µ . If the BSM states are charged under the electroweak (EW) gauge group, there is a vast landscape of models that can resolve the anomaly (see Refs. [47][48][49] for a review), hereafter referred to as electroweak (EW) models. Such models tend to have multiple free parameters and can, in principle, span a mass range for Λ anywhere between ∼ 100 GeV−100 TeV [47][48][49][50][51][52][53][54][55][56].
However, if the BSM states that generate Eq. (1.2) are singlets under the SM gauge group (singlet models), the Higgs VEV and chirality flip in Eq. (1.2) must both arise from the muon mass, so Eq. (1.2) becomes proportional to m µ , where M is the mass of a SM singlet particle that has been integrated out to generate these operators. This feature greatly simplifies the model landscape down to two renormalizable interactions: corresponding to BSM scalars (S) and vectors (V ), respectively, and greatly compresses the parameter space compared to EW models because of the muon mass suppression. 2 In this paper, we systematically analyze all singlet models that resolve the (g − 2) µ anomaly and survey how numerous proposed experiments can cover much of the remaining viable parameter space. Every singlet model can be uniquely specified by a dimensionless coupling g S,V , a corresponding singlet mass m S,V , and a possible additional decay width of the singlet to other non-muonic particles. Thus, singlet models are generically more predictive than electroweak models and can be studied without specifying many modeldependent features (for instance, the particle spectra of supersymmetric models). However, because the scalar singlet interaction in Eq. (1.4) is not gauge-invariant under the SM, it must be UV-completed, and the extra particles required for this UV completion will play 1 Lattice calculations of the Hadronic light-by-light contribution to aµ [31][32][33] as well as previous lattice calculations of the Hadronic Vacuum Polarization (HVP) [34][35][36][37][38][39][40][41][42] are in agreement with the phenomenological values. However, a more recent calculation of the HVP from the BMW collaboration [43] implies a value of aµ that is consistent with the measured value. This calculation might be in tension with the SM electroweak fit [44][45][46], for which future lattice calculations and improved R-ratio data will likely be needed to clarify this situation. 2 Alternative possibilities involving pseudo-scalars or axial vectors contribute the wrong sign to a BSM µ and do not resolve the anomaly, and interactions whose leading contribution to (g − 2)µ involves two-loop diagrams are severely constrained -see Appendix A for a discussion. an important role in further restricting the parameter space of the scalar singlet model once constraints from e.g. electroweak precision observables are taken into account. This paper is organized as follows. In Section 2 we summarize the landscape of singlet models and identify the allowed parameter space in each scenario, which is bounded from below by cosmology and from above by unitarity. In Sections 3, 4, and 5 we discuss future prospects for fixed target experiments, B-factories, and collider experiments, respectively. We conclude in Section 6 with an outlook on the types of experiments needed to fill any remaining gaps in parameter space to conclusively test singlet explanations for (g − 2) µ . Our main results are summarized in Figure 3. The Appendices contain discussions of 2-loop models, searches for singlets with non-muonic visible decays, and UV completions of the scalar singlet models.

Singlet Models
Throughout this paper, Singlet Models refers to the family of models where one of the interactions in Eq. (1.4) generates ∆a µ in Eq. (1.1). In this section, we outline the general properties of these models and identify their viable parameter space for resolving the (g−2) µ anomaly.

Vector Singlets
For a vector singlet with no axial couplings, the partial decay width to di-muons is for m V > 2m µ . The vector contribution to (g − 2) µ from Fig. 1 (left) can be written and the favored parameter space for resolving the (g − 2) µ anomaly is shown in the orange band of Fig. 2. In principle, there are many possible choices for abelian SM extensions that contain a vector-muon interaction. However, for anomaly-free U(1) gauge extensions to the SM, the favored parameter space for resolving (g − 2) µ is already excluded in most variants. Parameter space for which one SM singlet scalar or vector particle resolves the (g − 2) µ anomaly. The thickness of each band represents the ±2σ confidence interval and the vertical axis is the corresponding muon-singlet coupling from Eq. (1.4). The vertical shaded region represents the bound on light relativistic species present in equilibrium during big bang nucleosynthesis. The horizontal shaded region is the bound from partial wave unitarity. For a discussion of these bounds, see Section 2.3. Right: Minimum di-muon branching fraction for vector and scalar couplings that resolve the (g − 2) µ anomaly, from Eqs. (2.5) and (2.11).
• Kinetically-mixed dark photon The popular kinetically-mixed "dark photon" model features a new vector singlet corresponding to a secluded U(1) gauge group under which SM fermions carry no charges [57]. In this scenario the muon-dark photon coupling is induced through kinetic mixing between dark and visible photons, so the dark photon couples to all charged fermions with equal strength. Consequently, this model is excluded as a solution to (g − 2) µ by a variety of B-factory, beam dump, and rare meson decay searches -see [58,59] for a summary of constraints. This conclusion holds regardless of whether the dark photon decays visibly or invisibly. 3 • Anomaly-free U(1) extensions Minimal abelian extensions to the SM promote an anomaly-free linear combination of baryon (B), lepton (L), and lepton flavor (L i ) quantum numbers, where i = e, µ, τ , to a gauge symmetry. Some examples include B −L, B −3L i , L i −L j . Since nearly all of these options include some coupling to first-generation SM fermions, they are subject to the same constraints that apply to the dark photon and are, therefore, also excluded [61][62][63] (unless new fields are introduced [64][65][66][67][68][69][70]). A well known exception is gauged [71][72][73][74], which has no tree-level couplings to first-generation particles and evades many of the bounds that arise from couplings to protons and electrons. This model is viable for vector masses in the ∼ 15 -210 MeV range [72] in which V decays exclusively to neutrinos and can be tested with low-energy fixed-target experiments -see Sec. 3 for more details.
In the case of the most general U (1) X extension for (g − 2) µ including right-handed neutrinos (needed for anomaly cancellation) one can have a vast family of models [75]. These can be classified as dark photon-like or L µ −L τ -like in case they include (former) or not (latter) interactions with first generation leptons. Most of these models are either excluded or highly constrained [75].
• Anomalous U(1) extensions Since most bounds on anomaly-free U(1) models involve constraints on their nonmuonic interactions, in principle it is possible to open up parameter space by coupling the vector only to muons. However, such an approach introduces interactions that grow with energy via non-decoupling triangle diagrams that amplify otherwise rare SM processes [62,76,77]. While such models have not been studied exhaustively in the context of explaining (g − 2) µ , it is expected that they would be strongly constrained by existing bounds (e.g. rare meson decays involving muon final states). We note that a fully unitary theory with anomalous interactions at low energies must ultimately invoke additional SM charged states to cancel triangle diagrams at some scale; however, there is considerable model dependence in how this is achieved.
For the remainder of this paper, we will take an agnostic approach towards the models that realize the vector-muon coupling, while noting that specific realizations of vector singlets might have additional bounds beyond those that we present. Taking into account only the irreducible coupling to muons which generates (g − 2) µ , we will present a phenomenological treatment that studies only di-muon decays or invisible decays, where E also allows for a long-lived V that effectively gives rise to missing energy signatures. For m V > 2m µ in this agnostic approach, the total width is bounded by unitarity at [78]. Saturating this maximum width, even if the dominant V branching fraction is into other channels (visible or invisible), the minimum vector branching fraction to di-muons is where we have used Eq. (2.1) in the m V m µ limit and expressed Eq. (2.2) in terms of the central value in Eq. (1.1); this minimum value applies to any vector singlet that resolves the (g − 2) µ anomaly. In Fig. 2 (right panel) we show the minimum branching fraction as a function of singlet mass.
As an example of how this general analysis becomes constrained in the context of particular models, if V is the gauge boson of an anomaly-free group involving L µ , then the only viable parameter space is for m V < 2m µ with predominantly invisible V →νν decays, since m V > 2m µ is excluded by the irreducible visible decay channel. Other channels due to, for example, a small kinetic mixing between V and the photon, are strongly suppressed.

Scalar Singlets
For scalar singlets, the partial decay width to di-muons is for m S > 2m µ , and the corresponding contribution to (g − 2) µ from Fig. 1 (right) can be written (2.7) The favored parameter space for resolving the (g − 2) µ anomaly is shown in the green band of Fig. 2. This contribution is comparable to the corresponding expression for vector singlets in Eq. (2.2). Unlike in the vector case, the scalar Yukawa interaction Sµ L µ c in Eq. (1.4) is not gauge-invariant under SU(2) L × U(1) Y and must arise from a higher-dimension operator. In principle, this could be avoided if S mass-mixes with the SM Higgs. However, the SM Higgs-muon Yukawa coupling is y µ ≈ 6 × 10 −4 , so from Fig. 2, such a contribution is too small for all but the lightest S masses, even before taking into account strong suppression from the singlet mass or its tiny ( 1) mixing angle with the Higgs [79]. Thus, the lowestdimension operator that can generate a flavor-specific coupling to muons is where L is the second-generation lepton doublet and Λ is a BSM mass scale generated by heavier particles that are integrated out to generate the singlet-muon Yukawa coupling after electroweak symmetry breaking. We defer a detailed discussion of UV completions for the scalar singlet scenarios to Section 5.2 and Appendix C. Here we only point out that since scalar Yukawa interactions are not constrained by anomaly cancellation requirements, the flavor structure of these interactions is only constrained by experimental considerations: • Mass Proportional Couplings: As noted above, S cannot mass-mix with the SM Higgs because the muon Yukawa is too small to explain (g − 2) µ , but a larger muon coupling can be achieved if S is part of a two-Higgs doublet model with leptophilic couplings. In such a model, the scalar couplings are proportional to charged lepton masses and can be larger than their nominal SM values (see [80,81] for an example). However, the parameter space for resolving the (g − 2) µ anomaly within this model has been excluded by the BABAR experiment in a null search for BSM e + e − → τ + τ − µ + µ − events [82].
• Flavor Specific Couplings: The interaction in Eq. (2.8) can also arise by coupling the singlet to new electroweakcharged fermions that mix with the muon, followed by integrating them out at some scale Λ. In principle, such an approach allows any flavor structure (subject to empirical constraints) and many options have been considered in the literature [83][84][85][86][87][88][89][90].
Thus, for the remainder of this paper, we will not assume any particular flavor structure in the S couplings. For light scalars below the di-muon threshold, we will consider m S < 2m µ : S → γγ, E (BR = 1) . (2.9) As before, E represents invisible decays or a long-lived S. Note that for singlet scalars, the di-photon coupling arises at one-loop through the S coupling to muons, which is required for (g − 2) µ . 4 As we will see in Sec. 3.1, the decay length assuming only the γγ channel is macroscopic, such that it may manifest as missing energy E in a collider experiment but be visible at a beam-dump experiment. Note also that the γγ channel is forbidden for vectors by the Landau-Yang theorem. Since S → e + e − is not required by (g − 2) µ and depends on the flavor structure of the model, in our model-agnostic treatment we will not consider it in detail in this work, but we briefly comment on this channel in Appendix B. For scalars above the di-muon threshold, we will consider the irreducible decay channel As in the vector case discussed above, the minimum muon couplings required for (g − 2) µ and the unitarity upper limit on the singlet's total width Γ S,max ≈ m S /2 imply a minimum di-muon branching fraction independently of other decay channels where we have used Eq. (2.6) in the m S m µ limit and followed the same logic that leads to Eq. (2.5) in the vector case. In Fig. 2 (right panel), we show the minimum branching fraction as a function of singlet mass.

Viable Mass Range
Although the couplings required to explain (g−2) µ for light singlets are feeble from a collider perspective, they nonetheless suffice to thermalize these new states with the SM radiation bath in the early universe, provided the initial temperature exceeds a given choice of singlet mass. For light singlets (m S,V T ), the early universe production and annihilation rates scales as Γ S,V ∼ g 2 S,V T , where T is the temperature of the radiation bath. Comparing this rate to Hubble expansion implies thermalization for temperatures where g is the number of relativistic species in the bath and M Pl = 1.22 × 10 19 GeV is the Planck mass. Thus, for singlet masses below the MeV scale with couplings that resolve the (g − 2) µ anomaly in Fig. 2 (left panel), there will be a thermal ∝ T 3 number density Projections are indicated with colored lines. The M 3 [92], NA64µ [93], and ATLAS fixed-target [94] experiments probe invisibly-decaying singlets; projections here assume a 100% invisible branching fraction (see Sec. 3). The BABAR limits and Belle II projections are computed following the procedure described in Sec. 4. The LHC limits and HL-LHC projections in the 3µ/4µ channels, along with the mass range disfavored by UV completions for scalar singlets, are discussed in Sec. 5 and Appendix C. The purple muon collider projections are based on proposed analyses [48] reviewed in Sec. 5.3. For scalar singlets whose width is determined entirely by the muon coupling (top right), we also show projections for a S → γγ beam dump search [81] on the minimal assumption that the scalar-photon coupling arises from integrating out the muon as discussed in Sec. 3.1. Bottom: Same as the top row, only here we assume that for m S,V > 2m µ , the singlets have the minimum dimuon branching fraction consistent with unitarity using Eqs. (2.5) and (2.11). The curves which are unaffected by this change of muonic branching fraction correspond to searches that are insensitive to the singlet's decay modes. Projections for M 3 , NA64µ, and ATLAS fixed-target experiments assume a 100% invisible branching fraction for m S/V > 2m µ , which is model-dependent.
of singlets in the early universe at temperatures T ∼ MeV, which increases the expansion rate during Big Bang Nucleosynthesis (BBN) and spoils the successful SM prediction of light element yields. Avoiding the BBN bound requires [72,95] m V 10 MeV , m S 1 MeV , (2.13) which are represented in the shaded regions labeled "cosmology" in Figs. 2 and 3. For singlets above the muon mass, the BSM contribution to (g − 2) µ scales as a BSM µ ∝ g 2 S,V /m S,V so, in principle, any choice of mass can resolve the anomaly. However, for sufficiently large couplings, singlet interactions violate partial wave unitarity. As shown in [48], imposing tree-level unitarity constraints on muonic Bhabha scattering µ 12π. This implies an upper bound of m S < 2.7 TeV and m V < 1.1 TeV, or a UV completion involving additional states at the same mass scale as the singlet itself. The combination of cosmology and unitarity constraints provides a robust two-sided boundary to the viable mass range for singlet solutions to (g − 2) µ : (2.14) In Fig. 3, we show the viable parameter space for both visibly-and invisibly-decaying S and V , along with constraints and projections from fixed-target experiments, B-factories, and colliders, which we discuss in detail in the following sections.

Beam Dump Searches for S → γγ
In a beam dump experiment, a high-intensity muon beam impinges on a large, dense target which is surrounded by shielding material to block SM particles. If a long-lived, visiblydecaying singlet is produced in the dump, it can travel unimpeded through the shielding and decay to visible SM particles in a downstream detector. Such a setup can cover nearly all of the remaining parameter space for scalar singlets with m S < 2m µ that decay via S → γγ, induced at one-loop through the minimal scalarmuon Yukawa coupling required to address the (g − 2) µ anomaly [81]: where c is an order-one number and in the second expression the muon has been integrated out. If the S-µ Yukawa coupling is the only contribution to this operator, then for g S that Figure 1. Dark bremsstrahlung signal process for simplified models with invisibly decaying scalar (left) and vector (right) forces that couple predominantly to muons. In both cases, a relativistic muon beam is incident on a fixed target and scatters coherently o↵ a nucleus to produce the new particle as initial-or final-state radiation.
length in the nominal LDMX setup allows for a larger signal production rate while exploiting the fact that the muons will lose much less energy than electrons in a similarly-sized target. In analogy with similar processes involving electron beams, one can take advantage of the distinctive kinematics of the radiated massive scalar or vector particle S, V to distinguish signal from background (see Fig. 1). The Fermilab muon beam option provides several advantages over existing proposals for new physics searches with either electron beams or high-energy muon beams: • Bremsstrahlung backgrounds suppressed. The principal reducible backgrounds for LDMX are dominated by hadronic processes initiated by a real bremsstrahlung photon. Relative to electron beams, the M 3 bremsstrahlung rate is suppressed by (m e /m µ ) 2 ⇡ 2 ⇥ 10 5 , so background rejection becomes much simpler for muon beams for an equivalent target thickness.
• Compact experimental design. For m S,V ⌧ E beam , the signal production cross section is largely independent of beam energy. However, compared to the CERN/SPS option [11], with ⇠ 100 200 GeV beam muons, a lower-energy, e.g. 15 GeV, muon beam allows for greater muon track curvature and, therefore, a more compact experimental design. In particular, percentlevel momentum resolution is possible in M 3 with the target placed in the magnetic field region, reducing acceptance losses from having the magnet downstream of the target.
We propose a two-phase experiment, each covering a well-motivated region of parameter space: • Phase 1: (g 2) µ search. With 10 10 muons on target (MOT) and existing detector technology, we will show that our setup can probe the entire (g 2) µ region not currently excluded by experiments, for vectors with m V . 500 MeV and scalars with m S . 100 MeV which couple exclusively to muons and decay invisibly. 2 Here we are agnostic as to the UV completion of such a model, and we are simply aiming for an apples-to-apples comparison between a virtual S or V contributing to (g 2) µ and a real S or V emitted from an initial-or final-state muon.
• Phase 2: Thermal muon-philic DM search. With a larger flux of 10 13 MOT and upgraded detector performance to reject backgrounds at the level of 10 13 , our setup can probe a significant portion of parameter space for which DM is thermally produced through U (1) Lµ L⌧ gauge 2 Models with a more complicated dark sector can fail our search criteria, for example an inelastic DM model V ! 1 2 where the decay 2 ! 1e + e is prompt and proceeds through a di↵erent mediator which couples to electrons [18][19][20].
resolve (g − 2) µ anomaly, the decay length is [81] where E S is the energy of the S particle and the fiducial values for g S and m S are chosen to match the favored (g − 2) µ parameter space in Fig. 2

(left panel).
In the top-right plot of Fig. 3 we show the projections from Ref. [81] for a proposed beam dump search for long lived particles with a 3 GeV beam and 3×10 14 muons on target, assuming only the minimal muon-induced coupling to di-photons (Model B from Ref. [81]). However, note that if other model-dependent decay processes are introduced, the singlet S might not be a long-lived particle and this search strategy would lose sensitivity. Similar sensitivity could be achieved for the S → e + e − channel if the singlet-electron coupling is chosen to reproduce the same decay length as in Eq. (3.2), but unlike the minimal muon-induced di-photon decay, this displacement scale does not automatically follow from the same coupling that resolves the (g − 2) µ anomaly (also see Appendix B for further discussion).

Missing Momentum Searches for S/V → E
The proposed NA64µ [105,106] and M 3 [92] experiments aim to test the remaining (g − 2) µ parameter space for invisibly-decaying singlet particles. In these experiments, a muon beam of moderate intensity impinges on an active target, which monitors the beam energy as it passes through. Signal events are those where the outgoing muon energy is reduced by O(50%) and no other SM particles are observed in downstream veto detectors; the limiting factor for the luminosity is the requirement to only have a single muon per bunch to mitigate pileup. In Figure 3, we show the NA64µ and M 3 projections for 10 12 and 10 13 muons on target, respectively; these reference values match the ultimate reach projections from Refs. [105] and [92], where Ref. [105] only considered invisibly-decaying vectors. Note that these projections cover the remaining parameter space that resolves the (g − 2) µ discrepancy within the anomaly-free L µ − L τ gauge extension discussed in Sec. 2.1, where there is an order-1 invisible branching fraction to neutrinos for m V > 2m µ . Preliminary data is expected from NA64µ [111] in the very near future. Figure 5: Representative Feynman diagrams that yield singlet scalar (left) and vector (right) production at a B-factory via e + e − annihilation. The BABAR and Belle II search strategies discussed in Sec. 4 involve the 4µ channel in which S/V → µ + µ − decays yield 4µ final states with an opposite-sign di-muon resonance that reconstructs the singlet's mass.

B-Factories
We now discuss B-factory searches for singlets decaying to muons. Singlet models generate new contributions to e + e − → 4µ events from e + e − → µ + µ − +S/V → 4µ processes in which the singlet S or V is radiated from a final-state muon line and decays via S/V → µ + µ − , as depicted in Fig. 5. This search strategy can test the (g − 2) µ favored parameter space for light singlets with appreciable branching fractions to di-muons, but other search strategies are necessary for invisibly-decaying singlets. We briefly comment on other visible final-state searches at B-factories in Appendix B.
The BABAR collaboration has performed a search for muon-philic vectors with L µ −L τ gauge interactions and excluded the (g − 2) µ favored parameter space for masses between 2m µ and 10 GeV with 514 fb −1 of data [112]. Since an L µ − L τ gauge boson in this mass range can decay to muons, taus, and their corresponding neutrino flavors, in Fig. 3 (left panel) we rescale their limit for our vector V for which BR(V → µ + µ − ) = 1; for models in which V has a smaller branching fraction, the limit on g V gets weaker by a factor of BR(V → µ + µ − ).
In the right panel of Fig. 3, we reinterpret the BABAR limits from Ref. [112] in the context of our scalar singlet model. 5 Having already obtained the limit for vectors V , we can obtain the limit for scalars by applying the rescaling where σ is the singlet production cross section, S/V is the geometric acceptance, and each of the quantities here is mass-dependent. To obtain the necessary quantities for this rescaling, we simulated e + e − → µ + µ − S/V production events using MadGraph5 [113,114] which enable us to extract the cross sections and geometric acceptances assuming θ ∈ [18°, 150°] Figure 6: Singlet production at a hadron collider via the neutral (W -mediated) and charged (γ/Z-mediated) Drell-Yan processes. angular coverage in the BABAR lab frame, where θ is the polar angle of the final-state muon defined with respect to the incident electron beamline. We have also have verified that our simulated events adequately reproduce the L µ − L τ limits reported in Ref. [112] when adapted to that scenario.
In both panels of Fig. 3 we also show future projections for a 4µ search at Belle II with 50 ab −1 of luminosity. Here we assume that the BABAR background model and statistical uncertainties from Ref. [112] can be adapted to the Belle II experimental setup with a rescaled luminosity. To place a limit on the singlet couplings, we demand that the signal in each bin of di-muon reduced mass m R ≡ m 2 µ + µ − − 4m 2 µ not exceed the expected SM background by more than 2σ, where m µ + µ − is the reduced mass of an opposite state di-muon pair in the final state. A more rigorous, Belle II-specific systematic analysis is beyond the scope of this paper, but would be interesting for future work. Since the Belle II luminosity improves upon the BABAR data set by a factor of ∼ 100, this search alone will cover all remaining scalar singlets with mass between 2m µ and 10 GeV unless their branching fraction is to di-muons is less than ∼ 1%.

Collider Searches
We now discuss direct singlet searches at the LHC in the 3-and 4-muon final states, as well as inclusive singlet searches at future muon colliders. Since LHC searches have limited mass reach for scalars, it is important to take electroweak precision constraints from all possible types of UV completions for singlet scalar models into account, which sets a more restrictive upper mass limit than unitarity alone. Even so, we find that future muon colliders are required to exhaustively probe singlet explanations of the (g − 2) µ anomaly.

Singlet Production at the LHC
In order to produce S or V at a hadron collider, one needs to produce either one or two muons via charged (pp → ± ν ) or neutral (pp → + − ) Drell-Yan processes and then let one of the muons radiate a singlet. These processes are depicted in Fig. 6. Because the singlets can decay promptly back to a muon pair, the overall reactions look like threeand four-muon production, respectively. Backgrounds for these reactions include similar processes to those in Fig. 6, replacing the singlets by a γ/Z boson that further decays into a pair of muons, as well as Z + Z/W production where the vector bosons decay leptonically. These irreducible backgrounds combined can be sizeable compared to the signals, but with proper cuts on the invariant mass of opposite-charge di-muon pairs, one can isolate the signals with great sensitivity in the case where the singlets decay 100% to muons.
For our analysis, we implemented the Lagrangian Eq. (1.4) in FeynRules2 [115,116] and generated charged and neutral Drell-Yan events with MadGraph5 [113,114]. (For this clean muon signal, detector resolution and efficiency effects play a subdominant role and can be neglected in our estimates.) The cross sections for three-(left) and four-muon (right) production are presented in Fig. 7 where the SM background is shown in dashed gray and the vector and scalar singlets are shown in orange and green, respectively. These cross sections were obtained using the same loose cuts implemented in the four-lepton search by ATLAS in [117]: p µ T > 20 GeV (leading muon), p µ T > 10 GeV (sub-leading muon), p µ T > 5 GeV (the other muons), |η µ | < 2.7, ∆R µµ > 0.05, and m µµ > 5 GeV, where the small angular acceptance is required to allow topologies where a boosted particle decays into a pair of muons, and the cut in the invariant mass m µµ of any pair of opposite-sign muons is imposed to exclude leptons from quarkonia. These cross sections show that even without further cuts, the high-luminosity LHC (HL-LHC) should be able to probe singlet scalar masses up to m S m Z /2, but as we show below, the mass reach can be significantly increased by making use of the resonant nature of on-shell singlet production.
Note that in Fig. 7, the vector and the scalar production cross sections differ by about two orders of magnitude. This is for two reasons. First, the coupling g S/V fixed by (g − 2) µ is larger for vector singlets than for scalars by about a factor of three in the region above ∼ 1 GeV as shown in Fig. 2 (left). Second, due to the chiral flip at the scalar vertex in Eq. (1.4), the overall scalar production is a helicity-1 to a helicity-0 transition, so it is p-wave suppressed. This does not happen in vector models where an s-wave transition is possible, because the vector does not flip the chirality of the fermion lines. Also note that the charged Drell-Yan process (3µ) has a larger cross section than the neutral one (4µ). This is because W couplings to fermions are larger than Z couplings, and because the 3µ process includes two channels (µ + µ − µ ± ) whereas 4µ includes only one channel (µ + µ − µ + µ − ) and an extra suppression by statistical factors for identical particles in the final states.
Amongst recent LHC searches, the most relevant ones for this signal are two 139 fb −1 analyses conducted by ATLAS searching for new physics in 3µ and 4µ final states [117,118], as well as a 77 fb −1 Z → 4µ search for L µ −L τ gauge bosons by CMS [119]. The CMS search can be applied almost verbatim to our singlet vector scenario. Assuming the acceptances do not change significantly, we can also recast this search for the singlet scalar using the ratio between the vector and scalar singlet 4µ cross sections in Figure 7. The resulting limits are shown in Figure 3 and already exclude singlet scalar (vector) explanations of the (g − 2) µ anomaly in the mass range m S ∼ 6 − 30 GeV (m V ∼ 5 − 70 GeV), if these singlets decay exclusively to muons.
The mass reach can be significantly extended by dropping the requirement of an intermediate on-shell Z boson, and also considering 3µ final states. This means that recent multi-muon ATLAS analyses [117,118] have great exclusion power for these singlet scenarios. Unfortunately, these analyses do not explicitly look for an intermediate di-muon resonance. However, if the di-muon invariant mass distribution we describe below were made public for this or future such analyses, then these searches could supply the best limits on singlet scenarios that resolve the (g − 2) µ anomaly.
To demonstrate this, we derive sensitivity projections by considering the combined distribution of invariant masses of all possible opposite-sign di-muon pairs m µ + µ − in each event. Because there are multiple muons in the final states, the combinatorial background tends to spread the invariant mass spectrum over a wide range of values. Even so, the singlet-and Z-resonances remain clearly visible in the signal and background. For each singlet mass, we define a signal region in this combined m µ + µ − distribution that is centered on m S,V . The signal region has a variable width, ranging from 10 GeV for the lowest m S,V up to 30% of the singlet mass for larger values of m S,V . 6 Fig. 8 shows the LHC luminosity required using this strategy to exclude singlet scenarios that resolve the (g − 2) µ anomaly for masses above 5 GeV. This corresponds to a combined analysis of bump hunts in both charged and neutral Drell-Yan events, which we refer to from here on as a 3 + 4µ search. The 3µ events dominate the constraints for low singlet masses, whereas 4µ events dominate for large masses. The solid curve corresponds to the optimistic scenario where the singlets decay with a 100% branching ratio to muons, while the dashed curve corresponds to the minimum irreducible branching fraction to muons as shown in Fig. 2 (right).
In the optimistic scenario, one can see that 3+4µ searches at the LHC can probe singlet vector masses up to 700 GeV with current luminosity 140 ab −1 , whereas with the expected luminosity at HL-LHC the 3 + 4µ search could probe the entire parameter space up to the perturbativity limit at about 1 TeV. As discussed above, the signal cross section for scalars is much smaller than for vectors for a given singlet mass. For this reason the reach of the 3 + 4µ search is not as promising in the case of scalar singlets. With current (future) luminosity, the LHC could probe masses up to 50 (60) GeV. If the singlets only have the minimum irreducible branching ratio to muons, direct searches at the LHC have almost no power to probe this scenario. Even so, it is clear that this strategy greatly increases the mass reach compared to a search that is restricted to exotic Z decays. Figure 3 shows these HL-LHC projections as bounds on g S,V .

Electroweak Precision Constraints and Singlet Scalar UV Completions
As we can see from Fig. 8, singlet scalars can only be directly probed with 3 + 4µ searches at the HL-LHC if m S 60 GeV. It seems very difficult to directly probe heavier scalars coupling only to muons at a proton collider. Fortunately, we can exploit the fact that the scalar singlet model must have a UV completion featuring new electroweak charged states.
As discussed in Section 2.2, the scalar Yukawa interaction Sµ L µ c in Eq. (1.4) is not gauge-invariant and must arise from the dimension-5 operator where Λ is some BSM mass scale. We will assume that this operator is generated by treelevel exchange of some heavy mediator; if this interaction is somehow generated at higher loop level, the resulting electroweak states will be much lighter, resulting in even tighter constraints. There are only three possible types of diagrams that can generate this operator There are three possibilities to generate this operator via tree-level exchange of a mediator field, shown on the right. Assuming a minimal particle content for the additional fields, this corresponds to integrating out aut a fermion singlet (t-channel), a fermion doublet (t-channel), or a scalar doublet (s-channel).
for vectors for a given singlet mass. For this reason the reach of the 3 + 4µ search is not as promising in the case of scalar singlets. With current (future) luminosity, the LHC could probe masses up to 50 (60) GeV even in the optimistic scenario BR = 1.
We also show how the power of the 3 + 4µ search decreases once we assume a minimal branching ratio to muons as shown in Fig. 2 (right). Following the dashed lines, we can see that the luminosity required to see any deviation in 3 + 4µ searches is many orders of magnitude above HL-LHC for scalar singlets. For vector singlets, even assuming the minimal BR one can see deviations in 3 + 4µ at the HL-LHC except for a small region near the Z mass, and for singlet masses below 12 GeV. These projected constraints are also shown in Figure 3, translated into bounds on g S,V .

Electroweak precision constraints and UV completions
[DC: I rearranged this subsection a bit] As we can see from Fig. 7, singlet scalars can only be directly probed with 3 + 4µ searches at the HL-LHC if m S . 70 GeV. It seems very difficult to directly probe heavier singlets coupling only to muons at a proton collider. Fortunately, we can exploit the fact that the scalar singlet model must have a UV completion featuring new electroweak charged states.  There are three possibilities to generate this operator via tree-level exchange of a mediator field, shown on the right. Assuming a minimal particle content for the additional fields, this corresponds to integrating out a fermion singlet (t-channel), a fermion doublet (t-channel), or a scalar doublet (s-channel). at tree-level, shown in Figure 9: exchange of a fermion mediator that is an SU (2) L singlet, a doublet fermion mediator, or a doublet scalar mediator. The top two diagrams to the right correspond to what we call fermion UV completions and the bottom-right diagram is the scalar UV completion for the singlet scalar model. After electroweak symmetry breaking, the former class features heavy fermions that mix with the muon, whereas the latter features a heavy scalar doublet that mixes with the scalar singlet S.
It is easy to understand why electroweak precision tests (EWPT) [120] can constrain these UV completions. If the singlet accounts for the (g − 2) µ anomaly via an effective coupling to the muon g S ∼ v/Λ as in Figure 2, then heavier singlets require lower masses and/or larger couplings for the new mediator states. To derive indirect constraints on scalar singlet scenarios, we make the minimal assumption that the scenario is UV-completed by a single one of these three mediator types.
The main EWPT we focus on in our analysis is the modified lepton universality ratio In the SM, this ratio is close to unity except for small phase space differences and small loop corrections. Under the assumption that the singlet scalar resolves the (g−2) µ anomaly, both scalar and fermion UV completions generate significant Zµµ vertex corrections that can be probed by LEP measurements [120]. For scalar UV completions, the presence of additional scalars with electroweak charges also gives additional contributions to 3 + 4µ production at the LHC. The details of this analysis are provided in Appendix C. The important point is that we always consider the minimal size of R µe allowed by adjusting the various mediator couplings within the bounds set by unitarity, and that the simple singlet scalar model of Eq. (1.4) is a good IR effective theory in the non-excluded regions of parameter space. Our main results are the following: • Fermion UV completions of the singlet scalar violate EWPT for m S > 167 GeV.
• Scalar UV completions of the singlet scalar violate EWPT for m S > 615 GeV.
• If the singlet scalar only couples to the muons, then additional contributions to 3 + 4µ production at the LHC, together with EWPT, exclude m S < 63 GeV.
Since these constraints were derived by directly relating the required size of the dimension-5 operator Eq. (2.8) (to generate ∆a µ ) to the corresponding minimal Zµµ vertex correction, it is unlikely that these constraints would be greatly affected by considering more nonminimal UV completions, since additional couplings generally lead to additional observable effects unless deliberate cancellations are enforced. Therefore, it appears that with current LEP constraints, singlet scalar scenarios cannot be consistent with S masses above roughly 615 GeV, which has implications for future colliders that could aim to directly produce these singlets. We indicate this upper mass bound as the vertical dashed gray line in Figure 3.

Future Muon Colliders
Muon colliders have recently attracted considerable interest as successor machines to the LHC . Representative muon colliders concepts at the energy frontier may feature multi-TeV scale center-of-mass energies and luminosities of order ∼ ab −1 . Such facilities would be ideal laboratories for studying heavy, muon-philic interactions that cannot be probed robustly with existing proton or electron collider concepts [48]. Indeed, it has been shown that such machines can perform model-independent tests of the (g − 2) µ anomaly via muonic Bhabha scattering [48] and the µ + µ − → γh channel [159,160]. We have also shown in previous work that heavy singlets responsible for the (g − 2) µ anomaly can be directly tested through µ + µ − → γS/V production independently of how the singlets decay [48]. Since 2-body particle production yields back-to-back final-state recoils, this strategy triggers on events where the final-state photon is the only object in a detector hemisphere surrounding its momentum vector. Due to the large couplings required for heavy singlets (see Fig. 2), such production events can exceed the SM prediction for hemispherically-isolated single photons. Furthermore, since the observable is the isolated photon, the sensitivity of this strategy is independent of how the singlet decays as it recoils away from the photon. In Fig. 3 we show muon collider projections for singlet sensitivity assuming a 215 GeV (3 TeV) center-of-mass energy and 0.4 ab −1 (1 ab −1 ) of integrated luminosity as purple solid (dashed) curves. For the 215 GeV muon collider, singlets below ∼ 100 GeV are probed via direct production, while heavier singlets are probed through their virtual contributions to muonic Bhabha scattering (µ + µ − → µ + µ − ). At a 3 TeV muon collider, direct production is always more sensitive as long as the singlet can be produced on-shell, so the analogous Bhabha projections do not contribute to the exclusion curve we show.

Conclusions
In this paper we have systematically investigated solutions to (g − 2) µ in which a SM singlet generates the dominant contribution to the observed discrepancy via renormalizable interactions with muons. There are only two such classes of models, involving a singlet scalar S or a vector V , and both are highly restricted by general principles: scalars by the requirement of UV completion, and vectors by anomaly cancellation or non-decoupling triangle diagrams. Using only robust arguments from cosmology and unitarity, we constrain the viable mass range for singlets to lie between a few MeV and a few TeV [47,48]. We have then surveyed a program of existing and future experiments to cover this parameter space, illustrated in Figure 3 and summarized here: • Model-Independent Probes: A future 215 GeV muon collider could cover the entire (g − 2) µ favored region above 2m µ for vectors and about 1 GeV for scalars. A 3 TeV muon collider is less sensitive to light singlets, excluding all favored masses above about 1 GeV for vectors and 10 GeV for scalars. Notably, this is independent of how the singlet decays, as it only relies on the photon kinematics in µ + µ − → γS/V events [48], or singlet contributions to Bhabha scattering.
• Di-muon Decays: If the singlet decays predominantly to di-muons (when kinematically allowed), the (g − 2) µ parameter space can be probed with a combination of B-factory and highenergy collider searches. Current BABAR and LHC data, together with future Belle II and HL-LHC searches, can cover nearly the entire (g − 2) µ parameter space for vectors above 2m µ ; scalar singlets will be robustly covered for masses between 2m µ and 60 GeV.
• Invisible Decays: Missing energy/momentum experiments (including NA64µ, which will be taking preliminary data imminently) will cover the entire (g − 2) µ region below 2m µ , for both scalars and vectors which decay invisibly on the scale of the experiment, and have additional reach into the GeV scale if the dominant decay is invisible, as is true for viable anomaly-free vector models such as an L µ − L τ gauge boson.
• Other Visible Decays: The only remaining gap in the coverage presented here involves prompt visible decays to e + e − , γγ, and π + π − (or other visible exotics) for singlets in the ∼ MeV-10 GeV mass window; higher masses will be robustly covered with muon colliders, and displaced decays will be covered by missing-energy/momentum experiments at which long-lived particles yield missing energy. However, it is generically difficult to engineer prompt and visible decays for singlets that resolve (g − 2) µ in this mass range. All anomaly-free vector models with visible decays in this range have already been ruled out [59], and anomalous vectors are likely subject to severe constraints from non-decoupling triangle diagrams. Scalar singlets can be viable, but the electron coupling is generically required to be smaller than the muon coupling, which would yield macroscopic decay lengths for much of the allowed mass range. Similarly, long decay lengths are expected if the scalar decays through the minimal di-photon coupling induced by integrating out the muon [81]. We leave a more detailed study of this remaining parameter space with non-muonic decay channels for future work.
Thus, a broad program of energy and intensity frontier probes can robustly cover all the currently viable (g − 2) µ parameter space with only a few narrow exceptions. We have also studied how the singlet scalar model may be UV completed, and there are just three classes of possibilities if the singlet-muon coupling is generated by tree-level exchange of heavy mediators. Electroweak precision measurements then suggest that the singlet scalar mass has to be less than about 600 GeV, a factor of five lower than the upper bound provided by naive unitarity arguments alone.
Our analyses show that even without a muon collider, the parameter space of singlets that solve the (g−2) µ anomaly can be significantly probed by beam dumps, B-factories, and HL-LHC searches. However, this leaves large gaps in coverage. For the optimistic scenario where the new singlet decays entirely to muons, scalars above 60 GeV are inaccessible. In the pessimistic scenario where the singlet only has the irreducible minimum branching ratio to muons required by unitarity, vectors between ∼ 20 − 200 GeV are inaccessible, while scalars represent something of a nightmare scenario, apparently impossible to see above ∼ 1 GeV with this suite of experimental probes.
On the other hand, it is clear that a muon collider is required to probe singlet scalar solutions to the (g − 2) µ anomaly. They are straightforward to UV complete with masses much larger than the HL-LHC reach, and it is similarly easy to imagine scenarios where they have other decays that hide them from muon searches. Furthermore, probing singlet solutions particularly motivates a staggered approach in center-of-mass energy for the muon collider program. A 3 TeV machine is obviously preferred to push the boundaries of the energy frontier, as well as to start probing solutions to (g − 2) µ that involve dominant contributions by new electroweak charged states [48]. However, such a high-energy machine actually has more difficulty discovering light singlet scalars below 10 GeV, while direct production at a 215 GeV muon collider is sensitive down to the upper mass limit of proposed muon fixed-target experiments.
It is interesting to consider whether other experimental probes could fill the gap up to scalar masses of 10 GeV, so that no gap in coverage remains if only a multi-TeV muon collider were built. For both the scalar and the vector, the limited number of visible final states available in this narrow mass range should allow a complete enumeration of all of the possible visible and invisible decays for total widths approaching the unitarity limit. The outcome of such a study could have major implications for the future muon collider program by setting the minimum mass scale of new physics that such a high-energy machine would have to probe in order to leave no gaps in coverage. We leave this for future investigation.
The (g − 2) µ anomaly is quite possibly our best hint of new BSM physics in the laboratory. Our work highlights the need for muon colliders to comprehensively probe the new physics that must exist if the anomaly is real, from very high masses [47,48] all the way down to the GeV scale. No other kind of experiment has such a wide dynamic range in energy and precision, and if SM explanations of the anomaly are ultimately excluded, a comprehensive muon collider program would be well motivated as a top priority for the international particle physics community.

A Appendix: Two-Loop Models
Throughout this paper, we have studied singlet particles that couple linearly to muon currents and resolve the (g − 2) µ anomaly at one loop through the Feynman diagrams shown in Fig. 1. However, it is logically possible that the leading contribution to a µ involves two-loop diagrams instead. In this appendix we assess the feasibility of two representative two-loop scenarios. 7

A.1 Bilinear Scalar Couplings
In principle, a quadratic (or higher polynomial) scalar coupling could yield appreciable contributions to (g − 2) µ involving two-loop diagrams with no corresponding one-loop contributions. By SM gauge invariance, the lowest-dimension operator that yields a two-loop correction to (g − 2) µ involves a scalar S with the interaction where Λ is a mass scale that arises from integrating out particles charged under SU(2) L × U(1) Y and we assume that the simpler Sµ L µ c Yukawa coupling is negligibly suppressed or forbidden due to a Z 2 or other symmetry under which the S is charged. By naive dimensional analysis, the contribution to (g − 2) µ from the interaction in Eq. (A.1) is estimated to be where the prefactor (16π 2 ) −2 arises from the two-loop phase space and we demand Λ 90 GeV to avoid model-independent LEP bounds on new electroweak charged particles [161] while maximizing their contribution to (g − 2) µ . Furthermore, this correction to (g − 2) µ saturates to a fixed value for m S m µ , so the normalization in Eq. (A.2) represents the maximum contribution to the anomaly due to this interaction. Despite these extremal choices, the magnitude of the effect in Eq. (A.2) appears barely large enough to account for the necessary ∆a µ ≈ 2.5 × 10 −9 from Eq. (1.1) even if C eff is sizeable, but is sufficiently close that a more careful calculation would be needed before fully excluding this possibility.
However, the interaction in Eq. (A.1) assumes that the two-loop interaction arises from integrating out a single particle with a coupling-to-mass ratio of order Λ, without introducing any intermediate steps. In practice, resolving this higher-dimension operator with renormalizable and gauge-invariant interactions at energies above the scale Λ is more challenging and further suppresses the BSM contribution to (g − 2) µ , effectively resulting in C eff 1. Indeed, the only single-particle UV completion of Eq. (A.1) involves coupling S to the Higgs doublet H and integrating out the Higgs boson after EWSB to obtain where m h = 125 GeV is the Higgs mass, y µ ∼ 10 −4 is the SM muon Yukawa coupling, and κ is a dimensionless parameter. Here Λ → m h and C eff → κy µ , and for any unitary choice of λ HS , the a BSM µ contribution of Eq. (A.2) is much too small to explain the anomaly. Avoiding this Yukawa suppression requires adding additional SM charged states that couple to both the muon and the Higgs (for example, chiral fermions which mass-mix with the muon). These states must also couple to the SM-singlet operator S 2 , which is dimension-2 and cannot couple to a dimension-3 fermion bilinear with a renormalizable interaction, so an additional scalar must also be added and then integrated out to generate S 2 µ L µ c at low energies. However, if this additional scalar couples linearly to the µ L µ c bilinear, it introduces quantitatively larger one-loop contributions to (g − 2) µ , which invalidates our motivation; if its coupling to this bilinear arises from some other interaction, even more fields must be added to the theory and integrated out. Furthermore, despite these difficulties, the net effect of performing all these steps must still yield C eff v 2 /Λ 2 ∼ O(1) in Eq. (A.1), which seems implausible. While we have not presented a rigorous theorem to eliminate such a possibility, even if it were possible to generate a large prefactor and suppress all one-loop corrections to (g − 2) µ , its contribution would still be too small based on Eq. (A.2), subject to the same caveats mentioned above.

A.2 Millicharged Particles
It has also been argued that a two-loop contribution to (g − 2) µ could arise from a large N χ 1 multiplicity of MeV-scale fermions χ with electromagnetic millicharges 1 [162]. In this scenario, loops of χ introduce new QED-like vacuum polarization diagrams whose contribution resolves the anomaly if N χ 2 ∼ 10 −3 . However, such particles are subject to stringent limits from the E137 electron beamdump experiment [109,163], which were not considered in Ref. [162]. At E137, millicharged particles could be radiatively produced in the beam dump via electron-nucleus scattering e − N → e − N χχ, pass unimpeded through the beam dump, and scatter in the downstream detector. Since E137 reported null results, there are nontrivial limits on similar models involving MeV-scale scalars, as studied in Ref. [109]. It is expected that this limit would also impose nontrivial bounds on the millicharge model that explains (g − 2) µ with twoloop diagrams. Such an analysis is beyond the scope of this paper, but is worth studying in future work.
B Appendix: B-factory Searches for S/V → e + e − , γγ For m S,V < 2m µ it may be possible to perform a B-factory search for S → γγ and S/V → e + e − decays, but understanding the relevant SM backgrounds calls for a dedicated analysis, which is beyond the scope of this paper. Since the singlet coupling to electrons is, in principle, unrelated to the muon coupling that resolves the (g − 2) µ anomaly, the decays can be appreciably displaced at B-factory energies. For example, if the singlet decays dominantly to electrons, its decay length is approximately where we have used the rest frame width Γ(S → e + e − ) = g 2 S,e m S /(8π) and g S,e is the singlet coupling to electrons. Thus, a major challenge of this search strategy is distinguishing the displaced decay signal from SM photon conversion backgrounds in which an e + e − pair is produced in photon-nucleus interactions in the detector. Similar considerations apply to S → γγ decays which proceed through a higher-dimension operator (see Sec. 3.1) and are expected to be similarly long-lived if the only couplings in the model are the muon coupling and the loop-induced γγ coupling. We emphasize again that the singlet lifetime is unrelated to the singlet-muon coupling and can vary considerably, so it is not currently known whether a direct search strategy is possible for these specific final states.
In the context of discovering singlets responsible for (g − 2) µ , these visible non-muonic decay searches require m S < 2m µ where it is still possible for singlets to decay to e + e − or γγ with a large branching fraction; viable singlets above the di-muon threshold will always have a larger branching fraction to di-muons 8 or to invisible final states, in which case the 4µ search described in Sec. 4 or the invisible missing energy/momentum searches from Sec. 3.2 are better strategies, respectively. Furthermore, singlets with very displaced decays to e + e − , γγ final states can also be discovered or falsified at missing energy/momentum experiments if they decay downstream of the detector to fake a missing energy signature (see Sec. 3.2). Thus, the B-factory searches for e + e − or γγ final states described here are only sensible for singlets that decay promptly to these particles in the m S < 2m µ regime.

C Appendix: Scalar Singlet UV Completions
As discussed in Sec. 5.2, we will investigate the constraints on singlet scalars in the context of the fermion and scalar UV completions illustrated in Fig. 9. Some of these models were studied in detail in [89]. The Lagrangians for these interactions are for the fermion singlet mediator, fermion doublet mediator and scalar doublet mediator respectively (for simplicity we are suppressing the "+h.c." in all Lagrangians in this section).

C.1 Fermion UV Completions
The first two Lagrangians above can be expanded into where the second lines are written in four-component spinors. After integrating out χ and where M represents either the mass of χ or ψ d . Once the Higgs gets a VEV, the scalar singlet Yukawa interaction L ⊃ g S Sµ L µ c is generated, and we can identify the (g − 2) µ coupling as

Fermion Singlet UV Completion
The mass matrix for L I is given by whereỹ e = y e v/ √ 2 andỹ 1 = y 1 v/ √ 2. We can diagonalize this matrix with two real rotations M diagonal = L T M R, where L/R = cos θ L/R sin θ L/R − sin θ L/R cos θ L/R (C.7) and we can expand to obtain The couplings between the Z boson and the muon are Due to the mixing with χ the coupling g L gets shifted: This shift is what modifies the ratio R µe defined in Eq. (5.2) with respect to its SM value.

Fermion Doublet UV Completion
The mass matrix for L II is given by where this time the rotations are given by (C.14) In this case the right-handed Zµµ coupling gets shifted: Constraints Fig. 10 shows the SM expectation for the ratio R µe and the corresponding minimum deviation from the fermion UV completions L I (solid) and L II (dashed). These constraints were derived in the following way. For the L I UV completion, the relevant mixing angle of the new fermion with the muon is θ ∼ y 1 v/M , while g S = y 1 y 2 v √ 2M . To resolve the (g − 2) µ anomaly, g S = g S (m S ) is fixed as a function of the singlet mass, see Figure 2 (left). Therefore, for a given m S , the mixing angle is fully determined θ ∼ g S (m S )/y 2 , and is minimized by choosing y 2 to be as large as possible, in our case at the unitarity limit of √ 4π. For L II the argument is identical up to y 1 ↔ y 2 . The green lines in Fig. 10 translate into 2σ deviations in the ratio R µe for singlet masses m S above 167 GeV for the fermion UV completion L I , and above 181 GeV for the fermion UV completion L II .

C.2 Scalar UV Completion
The third Lagrangian in Eq. (C.1) can be expanded as follows: where φ u,d are the up and down components of the scalar doublet Φ. The parameter space of this UV completion is defined by the mass of the singlet m S , the mass of the doublet m Φ , the trilinear parameter κ, and the doublet-muon coupling y. For simplicity, we trade κ for the mixing parameter where ξ = 0 represents no mixing and ξ = 1 represents maximal mixing so that one of the eigenstates in the theory becomes massless. (We do not consider values of ξ > 1 for which the new scalars acquire VEVs and contribute to electroweak symmetry breaking.) For a given choice of m S (the mass parameter, not the mass eigenvalue) one can scan the plane (ξ, m Φ ) by fixing the coupling y so that a BSM µ = ∆a µ at each point in that plane. The summary of constraints on the parameter space of the scalar doublet UV completion is shown in Fig. 11. Grey regions at the bottom and top of the figures are excluded by perturbative unitarity in the y and κ couplings respectively. The region below the solid black line is excluded by EWPT. After mixing, this model contains four scalar eigenstates ϕ 1 (mostly S), ϕ 2 (mostly φ d CP-even), ϕ 3 (φ d CP-odd), and φ u (the charged component of Φ). The blue contours represent the mass of the ϕ 1 eigenstate and the orange contours show the coupling of this eigenstate to muons, corresponding to m S and g S in the singlet scalar effective theory (1.4). The regions to the left of the red lines would be excluded from The plots in Fig. 11 show that for m S < 81 GeV and m S > 830 GeV, the multiple constraints do not leave any viable parameter space for the model to generate (g − 2) µ , provided the new scalars do not have additional couplings to hidden sector states that reduce the branching fraction to muons. This is the reason why the scalar UV completion is only valid in the interval m 1 ∈ (63, 615) GeV, where m 1 is the mass of the eigenstate ϕ 1 . (Note the blue contours in Fig. 11 for the m S < 81 GeV and m S > 830 GeV panels.) If the singlet has additional couplings to invisible states, then the 3+4µ constraints would be weaker. We do not study this in detail, and therefore conservatively drop the lower bound on m 1 in this case. The upper bound is set by EWPT and is unaffected.
We now provide more details on each of the constraints shown in Fig. 11.
• Perturbative unitarity in y: In principle, the y interaction in Eq. (C.16) can produce the required contribution to (g − 2) µ by itself. This interaction requires a non-perturbative y coupling for masses above 560 GeV. By increasing the mixing between S and φ d via increasing ξ, one reduces the contribution to (g − 2) µ from the doublet and increases the contribution from the singlet S. This means that the larger the mixing, the larger the m Φ allowed by perturbativity of the coupling y. This explains the ascending behavior of the "y unitarity" constraint in the plots.
• EWPT: When the mixing ξ is small and the coupling y is large, the presence of the extra scalars ϕ i gives sizeable one-loop vertex corrections to the Zµµ coupling. Such a loop effect modifies the ratio R µe in a significant way for large enough values of the coupling y. This implies that EWPT constraints push the viable parameter space into the large mixing region, meaning large values of ξ to allow small values of y, except in two limiting regions: when m φ is small (hence y is also small), which suppresses the loop effect, and when m φ ∼ 1750 GeV where there is a dip due to an accidental cancellation between the different loop contributions to R µe .
• Perturbative unitarity in κ: For a given m S , increasing the mixing ξ implies decreasing the mass of the lightest eigenstate ϕ 1 while increasing κ. When a trilinear coupling is large compared to the mass scales in the theory, it can be constrained by perturbative unitarity. Our analysis is similar to the one in [47,48]. We consider the scattering of the lightest eigenstate ϕ 1 ϕ 1 → ϕ 1 ϕ 1 via Higgs exchange. From the scattering amplitude we calculate the partial wave expansion coefficient where p, k, s are the initial momentum, final momentum, and center of mass energy of the process, and θ is the scattering angle. All these quantities are defined in the center of mass frame. To find the constraints from unitarity on κ, for a given set of masses m S and m φ , we find the ξ value that saturates the unitarity condition |Re(a 0 )| < 1/2.
• LHC: Due to the four new scalars ϕ i one gets contributions to 3 + 4µ production in proton collisions. This happens through radiating some of the scalars ϕ i off muon lines in charged and neutral Drell-Yan production (through similar diagrams to those in Fig. 6) or by pair producing scalars that subsequently decay down to muons and/or neutrinos. Similar to the singlet-strahlung production analysis in Section 5.1, we derive the projected sensitivity of the LHC and HL-LHC. Unlike the singlet-only analysis we only require that the total BSM contribution to the 3+4µ production rate is smaller than a 2σ upward fluctuation of the SM contribution. This conservative choice reflects the fact that the BSM production of 3 and 4 muons includes both resonant and non-resonant processes, and a more comprehensive analysis of this scenario may yield stronger constraints. The BSM 3 + 4µ production cross section is dominated by the charged component of the doublet Φ, which is why the constraint vanishes for large m Φ . The zigzag behavior comes from the fact that the neutral scalars contribute and can in principle interfere with the contributions from the charged scalar in non-trivial ways. In any case, the trend holds that the constraints vanish for large m Φ . If the singlet has significant couplings to particles other than muons, e.g. invisible sector states, then these bounds would be weakened.