Exploring muonphilic ALPs at muon colliders

Axion-like particles (ALPs) are predicted to exist in a wide range of models that extend beyond the standard model (SM). The QCD axion, introduced originally to solve the strong CP problem, is one such model [1−5]. ALPs can also be generated from different spontaneous symmetry breaking patterns of global symmetries [6−9] as well as in string theory [10−13] and models of extra dimensions [14,15]. The broad spectrum of possible ALP masses makes them an attractive candidate for a variety of astrophysical and cosmological phenomena [16]. Sub-eV ALPs have been proposed as potential candidates for dark matter [17]. ALPs at different mass scales can also serve other purposes, such as acting as mediators to the dark sector [18,19], influencing the structure of the electroweak phase transition [20,21], and offering solutions to the hierarchy problem of the Higgs boson mass [22]. Understanding the characteristics and roles of ALPs is essential for unraveling the mysteries of the universe and advancing our knowledge of particle physics.

Various methods have been developed to search for ALPs, including laboratory-based experiments [23], astrophysical observations [24], and searches in high-energy collisions [25]. The current constraints on ALPs rely on their coupling strength and mass. For example, astrophysical observations of the diffuse gamma-ray background provide tight constraints on the coupling strength of sub-eV ALPs to photons [26−29], while experiments based on the LEP and LHC can limit the coupling strength of high-mass ALPs to SM particles [30−35]. With the advancements of experimental techniques, these bounds are expected to become even more stringent in the future, offering exciting new prospects for investigating the properties of ALPs.

In this study, we focus on muonphilic ALPs (μALPs), a specific type of ALP that predominantly interacts with muons [36−43]. These ALPs can be considered in an effective field theory framework [25,44−48], allowing us to study their couplings with SM particles independently. Bounds onμALPs for $ m_a < 2m_{\mu} $ have already been obtained from searches in supernovae [36,37,40] and atmospheric air showers [41]. For $ 2m_{\mu} < m_a \lesssim {\cal O}(1) $ GeV,μALPs can be largely produced in fixed target experiments [18], low-energy $ e^{+}e^{-} $ colliders [49], and Tera Z factories [43]. However, searching for GeV-scaleμALPs at high-energy colliders is challenging owing to the smallμALP production rate, as the coupling is proportional to the muon mass. Therefore, proposing newμALP production channels with sufficiently large cross sections at high-energy colliders is crucial to search for GeV-scaleμALPs.

Recently, a four-point interaction (W- $ \ell $ -ν-a), which has a coupling that is independent of the charged lepton mass, was proposed for the search for leptophilic ALPs [50]. This interaction is expected to arise from decays of $ \pi^{\pm} $ and $ K^{\pm} $ mesons and theWboson, with the novel energy enhancement effect. Similarly, this type ofW- $ \ell $ -ν-ainteraction with the energy enhancement effect has also been proposed as a promising approach to the search for leptophilic ALPs viat-channel processes ( $ \ell^+\ell^-\rightarrow \overline{\nu_{\ell}}a\nu_{\ell} $ and $ \ell^- p\rightarrow \nu_{\ell}aj $ ) at high-energy colliders [51]. In this study, we investigated the production of GeV-scaleμALPs from the above t-channel processes and their decay into a pair of muons at muon colliders [52−54]. Notably, when a lightμALP is highly-boost produced, the resulting pair of muons from theμALP decay is too collimated to pass the muon isolation criteria and forms a novel object known as a muon-jet [55−66].

We investigate three major signal processes at muon colliders: $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ . These signal production modes mainly rely on a four-point interaction,W-μ- $ \nu_{\mu} $ -a, and/or interactions from the chiral anomaly, whose couplings are not dependent on the muon mass. Generally, $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ yields the largest cross section, followed by $ \mu^{+}\mu^{-}\rightarrow \gamma a $ and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ in the electroweak violating (EWV) scenario. However, there is noW-μ- $ \nu_{\mu} $ -ainteraction in the electroweak preserving (EWP) scenario, and therefore $ \mu^{+}\mu^{-}\rightarrow \gamma a $ yields the largest cross section. In theEWVscenario, we discover that the channel $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ with theW- $ \ell $ -ν-ainteraction is the most important among these channels because of its novel energy-enhancement behavior. Our findings suggest that searching for the signature of two isolated muons (or a muon-jet) plus the missing energy in the EWV scenario at muon colliders can provide considerably stronger bounds than existing ones. On the other hand, searching for the signatures of two isolated muons (or a muon-jet) plus a photon and four isolated muons (or a muon-jet plus two isolated muons) in theEWPscenario at muon colliders may only slightly exceed existing bounds. Therefore, the muon collider is an ideal machine to search forμALPs and can also explore aμALP belonging to the EWV or EWP scenario.

This paper is organized as follows. In Sec. II, we provide a brief review of ALP-muon interactions andμALP decay modes. The method of distinguishing between different ALP-muon interaction types using the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ processes is discussed in Sec. III. We present the results of a full signal-to-background analysis at muon colliders and compare them with existing bounds of theμALP in Sec. IV. Finally, we summarize our findings in Sec. V. Supplementary materials, including kinematic distributions for both signals and SM backgrounds, and other tables are provided in Appendix A.

We consider ALPs generated from the breaking of the global Peccei-Quinn (PQ) symmetry [1], $ U(1)_{\text{PQ}} $ . Based on the structure of PQ symmetry, $ a(x)\rightarrow a(x)+\text{const} $ , the Lagrangian can be written in the form $ {\cal L}_{\mu\text{ALP}} = \partial_{\nu}a \; J^{\nu}_{ \text{PQ},\mu} $ . The general muon current is in the form

where Λ is the new physics scale, and $ c^V_{\mu} $ , $ c^A_{\mu} $ , and $ c_{\nu_{\mu}} $ are dimensionless couplings. Without the assumption of electroweak invariance, the condition $ c_{\nu_{\mu}} = c^V_{\mu} -c^A_{\mu} $ in Eq. (1) can be released ¹ . After integrating by parts for this Lagrangian, $ {\cal L}_{\mu\text{ALP}} $ can be written as [50]

where the symbols $ W^{\pm}_{\mu\nu} $ , $ Z_{\mu\nu} $ , and $ F_{\mu\nu} $ represent the field strength tensors of the massive gauge bosons $ W^{\pm} $ andZand the massless photon, respectively, the dual field strength tensor is defined as $ \tilde{F}_{\mu\nu}=\dfrac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma} $ , $ \alpha_{\text{em}} $ is the fine structure constant, $ g_W $ is the weak coupling constant, and $ s_W $ and $ c_W $ are the sine and cosine of the weak mixing angle, respectively.

In Eq. (2), we label the first term as " $ \boldsymbol{a\mu\mu} $ ," which can generateμALPs through muon radiation. However, this term is suppressed by $ m_{\mu}/\Lambda $ , necessitating high-intensity experiments to search for lightμALPs. The second to fourth terms, labeled as " $ \boldsymbol{aVV'} $ ," arise from the chiral anomaly and can produce lightμALPs through flavor-changing processes in meson decays [48,67,68]. HeavierμALPs can also be produced from these terms through gauge boson fusion and associated gauge boson production processes, despite not being proportional to $ m_{\mu} $ , but having $ \alpha_{\text{em}}/4\pi $ suppression. The terms in the final line of Eq. (2), labeled as " $ \boldsymbol{aW\mu\nu} $ ," are often overlooked in literature [69]. However, they are critical to our study, particularly when searching forμALPs at the GeV scale. This four-point interaction,W-μ- $ \nu_{\mu} $ -a, vanishes when the general muon current in Eq. (1) respects the electroweak symmetry. Moreover, this interaction is not related to $ m_{\mu} $ and has an obvious $ \left(\text{energy}/\Lambda\right) $ enhancement in specific processes. This enhancement behavior is crucial in constraining lightμALPs through decays of theWboson and charged mesons [50], as well as in searching for heavierμALPs int-channel processes such as $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ at muon colliders under the EWV scenario, which is defined in the next section.

However, searching forμALPs in collider experiments depends on their decay modes. ForμALP masses below the electroweak scale ( $ m_a\lesssim M_W $ ), their dominant decay modes are to $ \mu^{+}\mu^{-} $ and $ \gamma\gamma $ [25,45,70]. The decay widths are given by

where the coupling constant $ g_{a\gamma\gamma} $ is determined by the chiral anomaly and one-loop triangle Feynman diagrams and can be expressed as

and the loop function $ {\cal F} (z > 1) = \dfrac{1}{z}\text{arctan}^2\left(\dfrac{1}{\sqrt{1/z -1}}\right) $ . Here, we only consider the contribution from the muon loop because the contribution from theWboson is highly suppressed and can be safely neglected.

Figure 1shows the branching ratios for $ a\rightarrow \mu^{+}\mu^{-} $ and $ a\rightarrow \gamma\gamma $ . When $ m_a\lesssim M_W $ , the dominant decay mode ofμALPs is $ a\rightarrow \mu^{+}\mu^{-} $ . Because the partial decay width of $ a\rightarrow \gamma\gamma $ depends slightly on the muon mass and scales with $ m_a^3 $ , we can expect the branching ratio of $ a\rightarrow \gamma\gamma $ to increase with theμALP mass. It is important to note that this result is opposite to that of the electrophilic ALPs in Ref. [51] because the muon mass is considerably larger than the electron mass.

Figure 1.(color online) Decay branching ratios ofμALPs below the electroweak scale ( $ m_a\lesssim M_W $ ).

In this section, we focus on distinguishing between different types of ALP-muon interactions at muon colliders. First, we consider the signal process $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ with the relevant Feynman diagrams shown inFig. 2and numerically investigate the energy enhancement behavior of this process at muon colliders. To implement $ {\cal L}_{\mu\text{ALP}} $ from Eq. (2), we use FeynRules [71] and calculate cross sections for this process using Madgraph5_aMC@NLO [72] while varying the center-of-mass energy. As we know, the condition $ c_{\nu_{\mu}} = c^V_{\mu} -c^A_{\mu} $ is a criterion to determine whether ALP effective field theory is electroweak invariant or not. Therefore, we set $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ and $ c^V_{\mu} = c_{\nu_{\mu}}=0 $ as a benchmark point for the EWV scenario. Similarly, we set $ c^A_{\mu}/ \Lambda = c^V_{\mu}/ \Lambda = 10 $ TeV $ ^{-1} $ and $ c_{\nu_{\mu}}=0 $ as a benchmark point for the EWP scenario. We vary the center-of-mass energy $ \sqrt{s} $ between $ 1-15 $ TeV with $ m_a = 10 $ GeV at muon colliders.Fig. 3shows the energy enhancement behavior of cross sections in $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ , where the full contributions from the EWV and EWP scenarios are depicted by solid lines, and the contributions from $ \boldsymbol{aW\mu\nu} $ , $ \boldsymbol{aVV'} $ , and $ \boldsymbol{aW\mu\nu} $ in the EWV scenario are depicted by dashed lines.

Figure 2.(color online) Feynman diagrams for $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ . Here, the red marker indicates the $ \boldsymbol{aW\mu\nu} $ interaction, green indicates the $ \boldsymbol{aVV'} $ interaction, and blue indicates the $ \boldsymbol{a\mu\mu} $ interaction.

Figure 3.(color online) Energy enhancement behavior of cross sections in $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ with $ m_a = 10 $ GeV, $ c^A_{\mu}/\Lambda = 10 $ TeV $ ^{-1} $ , $ c^V_{\mu}=c_{\nu_{\mu}}=0 $ (EWV: solid-blue line) and $ c^A_{\mu} / \Lambda = c^V_{\mu}/ \Lambda = 10 $ TeV $ ^{-1} $ , $ c_{\nu_{\mu}}=0 $ (EWP: solid-olive line). In theEWVscenario, the dashed-red, dashed-green, and dashed-black lines are labeled as contributions from the $ \boldsymbol{aW\mu\nu} $ , $ \boldsymbol{aVV'} $ , and $ \boldsymbol{a\mu\mu} $ interactions, respectively.

As shown inFig. 3, the leading contribution in the EWV scenario originates from the $ \boldsymbol{aW\mu\nu} $ interaction, with the subleading contribution from the $ \boldsymbol{aVV'} $ interaction. The contribution from the $ \boldsymbol{aVV'} $ interaction is approximately seven orders of magnitude smaller than that from $ \boldsymbol{aW\mu\nu} $ , as depicted by the dashed lines inFig. 3. For $ \sqrt{s} = 1-2 $ TeV, the energy enhancement behaviors from these two interactions are evident because the momentum transfer size becomes sufficiently large, making the $ \boldsymbol{aW\mu\nu} $ and $ \boldsymbol{aVV'} $ interactions important. However, as the energy continues to increase, the growth rate becomes gentler because these two leading contributions steadily increase with the center-of-mass energy as $ \left(\text{energy}/\Lambda\right) $ .

Our numerical analysis reveals that the contribution from the $ \boldsymbol{aW\mu\nu} $ interaction is considerably greater than those from the $ \boldsymbol{aVV'} $ and $ \boldsymbol{a\mu\mu} $ interactions because of the novel energy enhancement behavior. Therefore, we show the analytical form for the amplitude square with the average (sum) over initial (final) polarization for the $ \boldsymbol{aW\mu\nu} $ interaction in the process $\mu^+(p_1)\mu^-(p_2) \rightarrow$ $ \nu_{\mu}(q_1)a(q_2)\overline{\nu_{\mu}}(q_3)$ ,

where $ s = \left(p_1 +p_2\right) ^2 = \left(q_1 + q_2 + q_3\right) ^2 $ , $ k = p_2 -q_3 $ , and $ k' = p_1 -q_1 $ . This shows that the amplitude square can be enhanced when the momentum transfer in thet-channel process is sufficiently large.

In the EWV scenario, the contribution from the $ \boldsymbol{a\mu\mu} $ interaction is negligible, whereas in the EWP scenario, both the $ \boldsymbol{aVV'} $ and $ \boldsymbol{a\mu\mu} $ interactions exist in $ \mu^{+}\mu^{-}\rightarrow \nu_\mu a\overline{\nu_\mu} $ . The cross sections have no obvious change with the increase in the center-of-mass energy in the EWP scenario owing to the lack of the energy enhancement effect. Lastly, the cross sections in the EWV scenario are more than six orders of magnitude larger than those in the EWP scenario for the $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ process, as shown inFig. 3. This is because there is an $ \boldsymbol{aW\mu\nu} $ interaction in the EWV scenario but not in the EWP scenario, and this interaction contributes to almost the entire cross-section amount in the EWV scenario. Therefore, this process is powerful in distinguishingμALPs in the EWV scenario from those in the EWP scenario.

Next, we discuss four optimal channels for searching forμALPs at muon colliders. Among these channels, we specifically consider theμALP couplings that are independent of the muon mass. TheseμALP production channels are $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ , and $ \mu^{+}\mu^{-}\rightarrow Z a $ . To compare theseμALP production channels at a muon collider with $ \sqrt{s} = 3 $ TeV, we calculate their cross sections in both the EWV and EWP scenarios with the benchmark point $ m_a = 10 $ GeV and $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ , as shown inTable 1. In the EWV scenario, we find that the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ channel has the largest cross section owing to its energy-enhancing behavior caused by theW-μ- $ \nu_{\mu} $ -ainteraction. The cross section of this channel is approximately six to seven orders of magnitude larger than those of other channels. However, in the EWP scenario, there is no energy-enhancing behavior in the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ channel, such that the cross section of this channel becomes smaller than those of other channels. At this point, the cross section of the $ \mu^{+}\mu^{-}\rightarrow \gamma a $ channel is the largest, making it the most prospective search channel in the EWP scenario.

Production channel	Cross section/fb
	$ {EWV} $	$ {EWP} $
$ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a \overline{\nu}_{\mu} $	$ 3.13\times 10^{4} $	$ 9.69\times 10^{-3} $
$ \mu^{+}\mu^{-} \rightarrow \mu^{+}\mu^{-}a $	$ 1.45\times 10^{-2} $	$ 1.69\times 10^{-2} $
$ \mu^{+}\mu^{-} \rightarrow a \gamma $	$ 7.72\times 10^{-2} $	$ 8.18\times 10^{-2} $
$ \mu^{+}\mu^{-} \rightarrow Z a $	$ 3.58\times 10^{-3} $	$ 2.54\times 10^{-2} $

Table 1.Cross sections of differentμALP production channels at a muon collider with $ \sqrt{s} = 3 $ TeV and the benchmark point $ m_a = 10 $ GeV and $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ .

In this section, we conduct a signal and background analysis on three specific search channels: the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu}_{\mu} $ , $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ processes. Our goal is to predict the future bounds for GeV-scaleμALPs at muon colliders and compare them with existing bounds.

To illustrate, we analyze the process $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu}_{\mu} $ in the EWV scenario and its relevant SM backgrounds in the context of the popular muon collider proposal with $ \sqrt{s}=3 $ TeV [53,54]. According toFig. 1, theμALP mainly decays into $ \mu^+\mu^- $ when $ m_a\lesssim M_W $ . Hence, we focus on the $ a\rightarrow \mu^+\mu^- $ decay mode in our analysis. TheμALP becomes highly boosted at the muon collider when it is sufficiently light; hence, $ \mu^{+} \mu^{-} $ in the final state may be too collimated to pass the muon isolation criterion at detectors. Taking a cone size $ R=0.1 $ as the muon isolation criterion at the muon collider, we find that a pair of muons cannot be isolated from each other at detectors when $ m_a \lesssim 15 $ GeV (parton-level), as shown in the left panel ofFig. 4. We can group this type of collimated non-isolated muon as a special signature "muon-jet" ( $ J_{\mu} $ ), which is a non-QCD jet-like structure that deposits most of its energy in the the muon spectrometer and has a distinct signature from QCD jets. Therefore, we classify the signal signatures into two categories: (1) two isolated muons plus missing energy ( $ {\:/ E} $ ) for $ m_a\gtrsim 15 $ GeV, and (2) $ J_{\mu} $ plus $ {\:/ E} $ for $ m_a \lesssim 15 $ GeV.

Figure 4.(color online) Left panel: Distribution of the opening angle between two muons, $ \Delta R_{\mu^{+}\mu^{-}} $ , from $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ at a muon collider with $ \sqrt{s}=3 $ TeV. Four benchmark mass values ofμALP, $ m_{a}=5 $ , $ 20 $ , $ 50 $ , $ 80 $ GeV, are displayed. Right panel: Distribution of theμALP lab frame decay length from $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ at a muon collider with $ \sqrt{s}=3 $ TeV. The benchmark points $ m_a =5 $ , $ 20 $ GeV with $ c^A_{\mu} / \Lambda = 0.1 $ TeV $ ^{-1} $ are considered.

For the first signal signature, two relevant SM backgrounds, $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $ and $\mu^{+}\mu^{-}\rightarrow t\overline{t}\rightarrow (b\mu^+\nu_{\mu}) (\overline{b}\mu^-\overline{\nu_{\mu}})$ , are considered. We choose the benchmark point $ m_a = 50 $ GeV with $ c^A_{\mu} / \Lambda = 0.1 $ TeV $ ^{-1} $ to display the signal features. To generate Monte Carlo samples for both signal and background processes, we use Madgraph5_aMC@NLO [72] and pass them to Pythia8 [73] for QED and QCD showering and hadronization effects. We impose pre-selection cuts ( $ P^{\mu}_T > 5 $ GeV and $ \lvert\eta_{\mu}\rvert < 2.5 $ ) at the parton-level for both the signal and backgrounds. To simulate the detector effects, we use the muon collider template in Delphes3 [74], whose muon isolation criterion is consistent with Refs. [75,76]. We use the Cambridge/Aachen ( $ C/A $ ) jet clustering algorithm [77,78] and consider ab-jet tagging efficiency of $ \epsilon_b = 0.8 $ with charm-jet and light-jet fake rates of $ P_{c\rightarrow b}=0.1 $ and $ P_{j\rightarrow b}=10^{-3} $ , respectively. The following event selections are required to identify the signal signature and suppress background events:

● (1) $ N(\mu)\ge 2 $ with $ P_T^{\mu_1} > 200 $ GeV, $ P_T^{\mu_2} > 10 $ GeV, $ \lvert\eta_{\mu_{1,2}}\rvert < 1.5 $ ,

● (2) $ 1500 < {\not E} < 2800 $ GeV and $ \lvert\eta_{{\:/ E}}\rvert < 1.8 $ ,

● (3) Veto $ N(b)\ge1 $ GeV with $ P_T^{b} > 25 $ GeV,

● (4) $ {\not E}/M_{\mu_1\mu_2} > 32 $ ,

● (5) $ \lvert M_{\mu_1\mu_2}-m_a\rvert < 2 $ GeV,

● (6) $ 3.0 < \Delta\phi_{\mu_1, {\:/ E}} < 3.3 $ and $ 2.9 < \Delta\phi_{\mu_2, {\:/ E}} < 3.5 $ ,

where $ P_T^{\mu_1} $ , $ P_T^{\mu_2} $ ( $ \eta_{\mu_1} $ , $ \eta_{\mu_2} $ ) are the transverse momentum (pseudorapidity) of leading and subleading energetic muons, $ {\not E} $ is the missing energy, $ M_{\mu_1\mu_2} $ is the invariant mass of a muon pair, and $ \Delta\phi_{\mu_i, {\not E}} $ is the azimuthal angle between the i-th muon and $ {\not E} $ . The cut-flow table including the signal and backgrounds for each event selection is shown inTable 2, and several kinematic distributions are shown inFig. A1of Appendix A.

Cut flow inσ/fb	Signal	$ \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $	$ t\overline{t} $
Generator	$ 2.54 $	$ 162.70 $	$ 4.15\times 10^{-2} $
cut-(1)	$ 1.78 $	$ 18.60 $	$ 7.94\times 10^{-3} $
cut-(2)	$ 1.78 $	$ 11.28 $	$ 7.21\times 10^{-3} $
cut-(3)	$ 1.78 $	$ 11.27 $	$ 3.54\times 10^{-4} $
cut-(4)	$ 1.74 $	$ 0.15 $	$ 2.12\times 10^{-5} $
cut-(5)	$ 1.47 $	$ 1.23\times 10^{-2} $	$ 8.30\times 10^{-7} $
cut-(6)	$ 1.35 $	$ 6.17\times 10^{-3} $	$ 0 $

Table 2.Cut-flow table for $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell} (a\rightarrow \mu^{+}\mu^{-})\overline{\nu_{\ell}} $ and relevant SM backgrounds with the signature of two isolated muons plus $ {\not E} $ . The benchmark point $ m_a = 50 $ GeV with $ c^A_{\mu} / \Lambda = 0.1 $ TeV $ ^{-1} $ is chosen for the signal. Each event selection is mentioned in the main text. The "Generator" means the cross sections at the parton-level calculated by Madgraph5_aMC@NLO.

First, we find two isolated muons and $ \not E $ in the central region of the signal events. To select candidate events, we apply the following trigger criteria: $ P_T^{\mu_1} > 200 $ GeV, $ P_T^{\mu_2} > 10 $ GeV, and $ \not E > 1500 $ GeV. InFig. A1, the distributions of $ P_T^{\mu_1} $ and $\not E$ show two peaks that correspond to the $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $ process. The right peak of the $\not E$ distribution indicates that most of the energy is carried away by the neutrino pair, leaving minimal energy for the two muons, whereas the left peak indicates that each of the two muons and two neutrinos carries almost an equal share of the energy. For the $ \mu^{+}\mu^{-}\rightarrow t\overline{t}\rightarrow (b\mu^+\nu_{\mu})(\overline{b}\mu^-\overline{\nu_{\mu}}) $ process, the $\not E$ distribution peak is around $ 2600 $ GeV, indicating that the two neutrinos take away more energy. As the signal $\not E$ distribution peak is around $ 2000 $ GeV, we apply $\not E < 2800$ GeV to reduce these two background events. Moreover, to suppress thebjet background from the $ \mu^{+}\mu^{-}\rightarrow t\overline{t}\rightarrow (b\mu^+\nu_{\mu})(\overline{b}\mu^-\overline{\nu_{\mu}}) $ process, we veto events with $ N(b)\ge1 $ GeV and $ P_T^b > 25 $ GeV. We also apply the ratio $\not E/M_{\mu_1\mu_2}$ as a complementary selection for theμALP mass window, setting $\not E/M_{\mu_1\mu_2} > 32$ . This selection is based on the observation that the position of the average $ M_{\mu_1\mu_2} $ distribution of the $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $ process is larger than that of the signal, and the $\not E$ distribution from this background is relatively small in the range $ 1500 < \not E < 2800 $ GeV. TheμALP mass window selection effectively reduces these two backgrounds while keeping most of the signal events. By applying the cuts of $ \Delta\phi_{\mu_{1,2}, {\not E}} $ to reduce some events from $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $ , we observe that two isolated muons are well-separated from $\not E$ . In particular, the distribution of $ \Delta\phi_{\mu_2, {\not E}} $ is not so large in both $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}\mu^{+}\mu^{-} $ and $ \mu^{+}\mu^{-}\rightarrow t\overline{t}\rightarrow (b\mu^+\nu_{\mu})(\overline{b}\mu^-\overline{\nu_{\mu}}) $ compared to the signal. Finally, using a benchmark integrated luminosity $ {\cal L} = 120 $ fb $ ^{-1} $ of a muon collider, we define the signal significanceZ[79] as

where $ N_s $ and $ N_b $ are the relevant signal and background event numbers. Here, the systematic uncertainties are not considered in our simple analysis because the muon collider is still a future collider. After all of the event selections inTable 2, we find that the signal significance can reach $ Z=38 $ for our benchmark point of $ {\cal L} = 120 $ fb $ ^{-1} $ , which means $ c^A_{\mu} / \Lambda < 0.1 $ TeV $ ^{-1} $ is still detectable in the future.

In the above analysis, we consider promptμALP decay with the lab frame decay length $ \gamma\beta c\tau_a < 1 $ mm as a criterion at a muon collider. Here,γis the Lorentz factor,βis theμALP velocity, and $ \tau_a $ is the proper decay time of theμALP. However, as we can expect, theμALP lab frame decay length becomes longer when $ m_a $ , $ c^A_{\mu}/ \Lambda $ are small andβis large. In this situation,μALPs become long-lived particles (LLPs). We take two benchmark points, $ m_a =5, 20 $ GeV with $ c^A_{\mu} / \Lambda = 0.1 $ TeV $ ^{-1} $ , to display the $ \gamma\beta c\tau_a $ distribution from $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ at the muon collider in the right panel ofFig. 4. We discuss the situation ofμALPs as LLPs in Sec. IV.C.

For the second signal signature, possible SM backgrounds originate from $ \nu_{\ell}\overline{\nu_{\ell}}b\overline{b} $ and $ \nu_{\ell}\overline{\nu_{\ell}}c\overline{c} $ , where the heavy flavor mesons produced fromc,bjets can decay into a collimated muon pair and mimic $ J_{\mu} $ from the signal. The pre-selection cuts ( $ P_T^{\mu} > 5 $ GeV, $ \lvert\eta_{\mu}\rvert < 2.5 $ ) at the parton-level are used for signal and background processes. We take the signal benchmark point as $ m_a = 5 $ GeV and $ c^A_{\mu}/ \Lambda = 0.1 $ TeV $ ^{-1} $ . The $ C/A $ jet clustering algorithm for $ J_{\mu} $ with a cone size $ R = 0.1 $ , which corresponds to the muon isolation criterion at the muon collider, is applied. We set up event selections to identify the signal signature and suppress the background events below:

● (1) $ N(\mu)\ge 2 $ with $ P_T^{\mu_{1,2}} > 5 $ GeV, $ \lvert\eta_{\mu_{1,2}}\rvert < 2.5 $ ,

● (2) $ N(J_{\mu}) = 1 $ and $ P_T^{J_{\mu}} > 20 $ GeV, $ \lvert\eta_{J_{\mu}}\rvert < 2 $ ,

● (3) $ 1500 < {\not E} < 2800 $ GeV and $ \lvert\eta_{{\:/ E}}\rvert < 1.4 $ ,

● (4) Veto $ N(b)\ge1 $ GeV with $ P_T^{b} > 25 $ GeV,

● (5) $ \lvert M_{J_{\mu}}-m_a\rvert < 2 $ GeV.

A cut-flow table including the signal and backgrounds for each event selection is shown inTable 3, and several kinematic distributions are shown inFig. A2of Appendix A.

Cut flow inσ/fb	Signal	$ \nu_{\ell}\overline{\nu_{\ell}}c\overline{c} $	$ \nu_{\ell}\overline{\nu_{\ell}}b\overline{b} $
Generator	$ 2.73 $	$ 208.20 $	$ 633.60 $
$ \gamma\beta c\tau_a< 1 $ mm	$ 0.52 $	$ - $	$ - $
cut-(1)	$ 0.50 $	$ 4.86\times 10^{-3} $	$ 0.17 $
cut-(2)	$ 0.50 $	$ 1.39\times 10^{-3} $	$ 2.41\times 10^{-2} $
cut-(3)	$ 0.47 $	$ 0 $	$ 6.31\times 10^{-3} $
cut-(4)	$ 0.47 $	$ 0 $	$ 1.27\times 10^{-3} $
cut-(5)	$ 0.42 $	$ 0 $	$ 0 $

Table 3.Similar to Table.II, but for the signal benchmark point $ m_a = 5 $ GeV and $ c^A_{\mu} / \Lambda = 0.1 $ TeV $ ^{-1} $ as well as the signature of a $ J_{\mu} $ candidate plus $\not E$ .

For theμALP prompt decay, we first set $ \gamma\beta c\tau_a < 1 $ mm as a criterion. Then, two muons with $ P_T^{\mu} > 5 $ GeV and $ \lvert\eta_{\mu}\rvert < 2.5 $ are required to be detectable in the muon spectrometer. We consider a $ J_{\mu} $ candidate with $ P_T^{J_{\mu}} > 20 $ GeV and $\not E > 1500$ GeV as the trigger, which are mainly distributed in the central region. The $ J_{\mu} $ in the signal events arises from energeticμALPs, whereas in the background events, it arises from the decay of heavy flavor mesons. As shown inFig. A2, $ P_T^{J_{\mu}} $ of the signal is considerably larger than that of the backgrounds, and most of the background events are largely reduced after cut-(2). We further require the selection $\not E < 2800$ GeV, which retains most signal events while removing significant parts of the background events; in particular, the events from $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}c\overline{c} $ are entirely removed. To suppress $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}b\overline{b} $ , we veto $ N(b)\ge1 $ GeV with $ P_T^{b} > 25 $ GeV. We also require the jet mass of $ J_{\mu} $ to satisfy theμALP mass window selection, which can entirely remove events from $ \mu^{+}\mu^{-}\rightarrow \nu_{\ell}\overline{\nu_{\ell}}b\overline{b} $ . After all event selections inTable 3, we can take this signal benchmark point as background-free. The distribution of the peak of $ M_{J_{\mu}} $ is broader than that of $ M_{\mu_1\mu_2} $ because two muons within a $ J_{\mu} $ cannot pass the muon isolation criteria. The selection of the jet clustering method, in conjunction with the choice of cone size $ R=0.1 $ , can affect the four-momentum reconstruction of $ J_{\mu} $ . In some cases, one of the muons is outside the jet cone and cannot be reconstructed, leading to distortions in $ M_{J_{\mu}} $ compared to $ M_{\mu_1\mu_2} $ . With $ {\cal L} = 120 $ fb $ ^{-1} $ , there are 50 signal events remaining for this benchmark point after all event selections.

In the EWP scenario, we employ a different approach to search forμALPs compared to the EWV scenario. As explained toward the end of Sec. III, the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ process in the EWP scenario lacks energy-enhancement behavior, leading to a smaller production cross section. For this reason, we focus on the twoμALP production channels $ \mu^{+}\mu^{-}\rightarrow \gamma a $ and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ , which have larger production cross sections, for the signal-to-background analyses to obtain stronger future bounds. We analyze the process $ \mu^{+}\mu^{-}\rightarrow \gamma a $ ( $ a\rightarrow \mu^{+}\mu^{-} $ ) using the same method as in Sec. IV.A. The details and results are presented below. The signal signatures are first classified into two categories: (1) two isolated muons plus a photon (γ) for $ m_a\gtrsim 15 $ GeV, and (2) $ J_{\mu} $ plusγfor $ m_a \lesssim 15 $ GeV. To investigate the first signal signature, we consider the relevant SM background $ \mu^{+}\mu^{-}\rightarrow \gamma\mu^{+}\mu^{-} $ and choose the benchmark point $ m_a = 50 $ GeV with $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ to display the signal features. The following event selections are required to identify the signal signature and suppress background events:

● (1) $ N(\mu)\ge 2 $ with $ P_T^{\mu_1} > 100 $ GeV, $ 10 < P_T^{\mu_2} < 500 $ GeV, $ \lvert\eta_{\mu_{1,2}}\rvert < 1.5 $ ,

● (2) $ E_{\gamma} > 1450 $ GeV and $ \lvert\eta_{\gamma}\rvert < 1.6 $ ,

● (3) $ 2.9 < \Delta\phi_{\mu_1, \gamma} < 3.4 $ and $ 2.9 < \Delta\phi_{\mu_2,\gamma} < 3.3 $ ,

● (4) $ E_{\gamma}/M_{\mu_1\mu_2} > 29 $ ,

● (5) $ \lvert M_{\mu_1\mu_2}-m_a\rvert < 2.0 $ GeV.

Here, $ E_{\gamma} $ is the energy of the photon, and $ \Delta\phi_{\mu_i,\gamma} $ is the azimuthal angle between thei-th muon and the photon. The cut-flow table detailing the signal and background for each event selection is presented inTable 4, with several relevant kinematic distributions shown inFig. A3of Appendix A. First, we find that two isolated muons in the signal are predominantly located in the relatively low transverse momentum regions (as shown inFig. A3), whereas a number of signal events for $ E_{\gamma} $ are mainly located in the relatively higher energy regions. This is because two muons in the signal are produced from the decay of theμALPs, which are secondary particles. By contrast, two muons from the background events mainly originate from the initial particles. Similarly, the photon in the signal is produced from the initial muons and therefore becomes more energetic. However, for the background process, the energy of photons is more divided by $ P_{z} $ of muons; hence, the leading photon energy is smaller than that of the signal process, as shown inFig. A3. To select candidate events, we apply the following trigger criteria: $ P_T^{\mu_1} > 100 $ GeV, $ P_T^{\mu_2} > 10 $ GeV, and $ E_{\gamma} > 1450 $ GeV. To reduce background events, we apply the cut of $ \Delta\phi_{\mu_{1,2},\gamma} $ because we observe that two isolated muons are well-separated from the photon. Additionally, we implemented a complementary selection based on the ratio $ E_{\gamma}/M_{\mu_1\mu_2} $ to further reduce the contribution from $ \mu^+\mu^-\rightarrow \gamma\mu^+\mu^- $ . Specifically, we set the ratio $ E_{\gamma}/M_{\mu_1\mu_2} > 29 $ , which effectively reduces background events while retaining most of the signal events. Since the average position of the $ M_{\mu_1\mu_2} $ distribution for the $ \mu^{+}\mu^{-}\rightarrow \gamma\mu^{+}\mu^{-} $ process is considerably wider than that of the signal, theμALP mass window can further reduce the number of background events. Finally, we consider a benchmark integrated luminosity of a muon collider of $ {\cal L} = 1000 $ fb $ ^{-1} $ to our analysis. After all of the event selections inTable 4, we find that the signal significance can reach $ Z=6.379 $ .

Cut flow inσ/fb	Signal	$ \mu^{+}\mu^{-}\rightarrow \gamma\mu^{+}\mu^{-} $
Generator	$ 6.84\times 10^{-2} $	$ 179.80 $
cut-(1)	$ 2.70\times 10^{-2} $	$ 4.72 $
cut-(2)	$ 2.69\times 10^{-2} $	$ 0.98 $
cut-(3)	$ 2.48\times 10^{-2} $	$ 0.56 $
cut-(4)	$ 2.03\times 10^{-2} $	$ 2.70\times 10^{-2} $
cut-(5)	$ 1.68\times 10^{-2} $	$ 3.06\times 10^{-3} $

Table 4.Cut-flow table for $ \mu^{+}\mu^{-}\rightarrow \gamma a $ and the relevant SM background with the signature of two isolated muons plus a photon. The benchmark point $ m_a = 50 $ GeV with $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ is chosen for the signal.

For the second signal signature, possible SM backgrounds originate from $ \gamma b\overline{b} $ and $ \gamma c\overline{c} $ , where the heavy flavor mesons generated fromcandbjets can decay into a collimated muon pair and mimic $ J_{\mu} $ from the signal. The pre-selection cuts ( $ P_T^{\mu} > 5 $ GeV, $ \lvert\eta_{\mu}\rvert < 2.5 $ ) at the parton-level are used for the signal and background processes. We still choose the signal benchmark point as $ m_a = 5 $ GeV and $ c^A_{\mu}/ \Lambda = 10 $ TeV $ ^{-1} $ . Applying the $ C/A $ jet clustering algorithm for $ J_{\mu} $ with a cone size $ R = 0.1 $ , we implement event selection criteria to isolate the signal and suppress background events as specified below:

● (1) $ N(\mu)\ge 2 $ with $ P_T^{\mu_{1,2}} > 5 $ GeV, $ \lvert\eta_{\mu_{1,2}}\rvert < 2.5 $ ,

● (2) $ N(J_{\mu}) = 1 $ and $ P_T^{J_{\mu}} > 750 $ GeV, $ \lvert\eta_{J_{\mu}}\rvert < 1.5 $ ,

● (3) $ E_{\gamma} > 1400 $ GeV and $ \lvert\eta_{\gamma}\rvert < 1.0 $ ,

● (4) $ \lvert M_{J_\mu}-m_a\rvert < 3.0 $ GeV,

● (5) $ 100 < E_{\gamma}/M_{J_\mu} < 400 $ .

The reasons for implementing these event selection criteria are similar to those described earlier and will not be repeated here. A cut-flow table detailing the signal and backgrounds for each event selection is presented inTable 5, with several relevant kinematic distributions shown inFig. A4of Appendix A. Employing all of the event selections listed inTable 3, we find that the signal significance can achieve a value of $ Z=20.10 $ with $ {\cal L} = 1000 $ fb $ ^{-1} $ .

Cut flow inσ/fb	Signal	$ \gamma c\overline{c} $	$ \gamma b\overline{b} $
Generator	$ 8.03\times 10^{-2} $	$ 5.96 $	$ 8.48 $
cut-(1)	$ 7.78\times 10^{-2} $	$ 4.70\times 10^{-3} $	$ 1.78\times 10^{-2} $
cut-(2)	$ 6.36\times 10^{-2} $	$ 7.63\times 10^{-4} $	$ 5.49\times 10^{-4} $
cut-(3)	$ 5.29\times 10^{-2} $	$ 1.19\times 10^{-5} $	$ 1.50\times 10^{-4} $
cut-(4)	$ 4.43\times 10^{-2} $	$ 0 $	$ 1.33\times 10^{-4} $
cut-(5)	$ 3.61\times 10^{-2} $	$ 0 $	$ 4.99\times 10^{-5} $

Table 5.Similar toTable 4, but for the signal benchmark point $ m_a = 5 $ GeV and $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ as well as the signature of a $ J_{\mu} $ candidate plus a photon.

In the subsequent section, we analyze the process $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ (including $ \mu^{+}\mu^{-}\rightarrow Z a\rightarrow (\mu^+\mu^-) a $ ) using the same method as in Sec. IV.A. The signal signatures are classified into two categories: (1) four isolated muons for $ m_a \gtrsim 15 $ GeV, and (2) $ J_{\mu} $ plus two isolated muons for $ m_a \lesssim 15 $ GeV. To investigate the first signal signature, we consider the relevant SM background $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-}\mu^{+}\mu^{-} $ and choose the same signal benchmark point $ m_a = $ 50 GeV with $ c_{\mu}^A / \Lambda = 10 $ TeV $ ^{-1} $ to display the signal features. The following event selections are required to identify the signal signature and suppress background events:

● (1) $ N(\mu)\ge 4 $ with $ P_T^{\mu_{1,2,3,4}} > 5 $ GeV, $ \lvert\eta_{\mu_{1,2,3,4}}\rvert < 2.5 $ ,

● (2) $ P_T^{\mu_1} > 200 $ GeV, $ P_T^{\mu_{2,3}} > 100 $ GeV, $ \lvert\eta_{\mu_{1}}\rvert < 2.0 $ and $ \lvert\eta_{\mu_{4}}\rvert < 1.5 $ ,

● (3) $ \Delta\phi_{\mu_2,\mu_4} > 0.5 $ and $ \Delta\phi_{\mu_3,\mu_4} > 1 $ ,

● (4) $ P_T^{\mu_4}/M_{\mu_1 \mu_4} > 0.05 $ ,

● (5) $ \lvert M_{\mu_1\mu_4}-m_a\rvert < 5.0 $ GeV.

A cut-flow table including the signal and background for each event selection is shown inTable 6, and several kinematic distributions are shown inFig. A5of Appendix A.

Cut flow inσ/fb	Signal	$ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-}\mu^{+}\mu^{-} $
Generator	$ 1.06\times 10^{-3} $	$ 2.82 $
cut-(1)	$ 1.03\times 10^{-3} $	$ 1.26 $
cut-(2)	$ 6.53\times 10^{-4} $	$ 3.26\times 10^{-1} $
cut-(3)	$ 4.74\times 10^{-4} $	$ 1.59\times 10^{-1} $
cut-(4)	$ 4.35\times 10^{-4} $	$ 5.80\times 10^{-2} $
cut-(5)	$ 3.43\times 10^{-4} $	$ 3.67\times 10^{-3} $

Table 6.Similar toTable 4, but for the $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ channel and relevant SM backgrounds with the signature of four isolated muons.

Four isolated muons with $ P^{\mu}_T > 5 $ GeV and $ \lvert\eta_{\mu}\rvert < 2.5 $ are applied as trigger criteria. We observed that $ P^{\mu_1}_T $ , $ P^{\mu_2}_T $ , and $ P^{\mu_3}_T $ of the signal are more energetic than those of the background. However, $ \mu_1 $ and $ \mu_4 $ from the signal are distributed in the central regions relative to those from the background, whose $ \mu_1 $ and $ \mu_4 $ are mainly generated by the initial muons with forward and backward directions. Therefore, we choose cut-(2) in the above to select candidate events. We check all combinations of four muons in the final state to reconstruct a pair of muons that comes from theμALP decay and find that the pair of $ \mu_1 $ and $ \mu_4 $ is most likely to reconstruct the mass ofμALP. To reduce background events, we apply the event selections based on $ \Delta\phi_{\mu_{2,4}} $ and $ \Delta\phi_{\mu_{3,4}} $ . We observe that the muons produced by the decay ofμALPs in the signal are well-separated from $ \mu_2 $ and $ \mu_3 $ . Furthermore, we incorporate an additional selection criterion involving the ratio $ P_T^{\mu_4}/M_{\mu_1 \mu_4} $ to further suppress the contribution from $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-}\mu^{+}\mu^{-} $ . Specifically, we set the ratio $ P_T^{\mu_4}/M_{\mu_1 \mu_4} > 0.05 $ , which effectively reduces background events while retaining the majority of signal events. Since the $ M_{\mu_{1}\mu_{4}} $ distribution for the $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-}\mu^{+}\mu^{-} $ process is almost concentrated in the region greater than $ 85 $ GeV, theμALP mass window can further reduce the number of background events. Finally, after all of the event selections inTable 6, we find that the signal significance can reach $ Z=10.44 $ for $ {\cal L} = 1000 $ fb $ ^{-1} $ in our analysis.

For the second signal signature, possible SM backgrounds arise from $ \mu^{+}\mu^{-}c \overline{c} $ and $ \mu^{+}\mu^{-}b \overline{b} $ . The pre-selection cuts ( $ P_T^{\mu} > 5 $ GeV, $ \lvert\eta_{\mu_{}}\rvert < 2.5 $ ) at the parton-level are used for the signal and background processes. We still choose the signal benchmark point $ m_a = 5 $ GeV and $ c^A_{\mu} / \Lambda = 10 $ TeV $ ^{-1} $ . Applying the same $ C/A $ jet clustering algorithm for $ J_{\mu} $ with a cone size $ R = 0.1 $ , we set up event selections to pick up the signal and suppress background events as specified below.

● (1) $ N(\mu)\ge 4 $ with $ P_T^{\mu_{1,2,3,4}} > 5 $ GeV, $ \lvert\eta_{\mu_{1,2,3,4}}\rvert < 2.5 $ ,

● (2) $ P_T^{\mu_{2,3}} > 100 $ GeV, $ \lvert\eta_{\mu_{2,3}}\rvert < 2.0 $ ,

● (3) $ N(J_{\mu}) = 1 $ with $ 300 < P_T^{J_{\mu}} < 1400\; {\rm GeV}$ , $ \lvert\eta_{J_\mu}\rvert < 1.6 $ ,

● (4) $ \Delta\phi_{J_{\mu},\mu_2} > 4.5 $ and $ 1.0 < \Delta\phi_{\mu_{2},\mu_3} < 5.0 $ ,

● (5) $ P_{T}^{\mu_2}/M_{J_{\mu}} < 500 $ ,

● (6) $ \lvert M_{J_{\mu}}-m_{a}\rvert < 4 $ GeV.

The event selection criteria are implemented for reasons similar to those previously described and will not be reiterated here. The cut-flow table, which includes both signal and background events for each selection, is presented inTable 7, and the corresponding kinematic distributions can be found inFig. A6of Appendix A. Upon applying all of the event selections outlined inTable 7, we observe a significant signal significance of $ Z=19.86 $ for our benchmark point with $ {\cal L} = 1000 $ fb $ ^{-1} $ .

Cut flow inσ/fb	Signal	$ \mu^{+}\mu^{-}c \overline{c} $	$ \mu^{+}\mu^{-}b \overline{b} $
Generator	$ 1.26\times 10^{-3} $	$ 52.94 $	$ 90.18 $
cut-(1)	$ 1.13\times 10^{-3} $	$ 1.98\times 10^{-3} $	$ 7.01\times 10^{-2} $
cut-(2)	$ 7.19\times 10^{-4} $	$ 5.82\times 10^{-4} $	$ 2.12\times 10^{-2} $
cut-(3)	$ 6.63\times 10^{-4} $	$ 2.65\times 10^{-4} $	$ 1.52\times 10^{-2} $
cut-(4)	$ 6.61\times 10^{-4} $	$ 1.06\times 10^{-4} $	$ 1.25\times 10^{-2} $
cut-(5)	$ 6.03\times 10^{-4} $	$ 0 $	$ 8.21\times 10^{-3} $
cut-(6)	$ 5.39\times 10^{-4} $	$ 0 $	$ 5.41\times 10^{-4} $

Table 7.Similar toTable 5, but for the $ \mu^+\mu^-\rightarrow \mu^+\mu^- a $ channel and relevant SM backgrounds with the signature of a $ J_{\mu} $ candidate plus two isolated muons.

The study of signal benchmark points is extended to a wide range of $ m_a $ by employing the search strategies outlined in Sec. IV.A and Sec. IV.B, resulting in the identification of possible future bounds within a $ 95$ % confidence level (CL) ( $ Z=1.96 $ ). To conservatively demonstrate the signal significance of the case without the survival background event after all event selections or the case of a background-free assumption, a minimum of $ 10 $ signal events is required to be present, and only signal efficiencies larger than $ 10$ % are considered in the analysis. Our study is restricted to $ 1 $ GeV $ \le m_a \lesssim M_W $ forμALPs below the electroweak scale. We first summarize our results inFig. 5to search the $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a \overline{\nu}_{\mu} $ channel in the EWV scenario. The dotted lines are used for the case of promptμALP decay ( $ \gamma\beta c\tau_a < 1 $ mm) at a muon collider with $ {\cal L} = 120 $ fb $ ^{-1} $ . Note that the lower bound of $ m_a $ originates from technical issues in the $ J_{\mu} $ analysis. When $ m_a < 1 $ GeV, theμALP mass window selection is no longer sufficiently powerful to distinguish the signal from backgrounds. However, since the event selections for two isolated muons plus $ {\not E} $ at the muon collider are sensitive to the values of $ m_a $ , event selections are dynamically optimized for different $ m_a $ to suit each case, as shown inTable 8. The case of prompt decay ofμALPs with two isolated muons at the muon collider yields a background-free scenario when $ m_a\lesssim 30 $ GeV, causing a cusp point at $ m_a = 30 $ GeV inFig. 5. Similarly, for the prompt decay ofμALPs with $ J_{\mu} $ at the muon collider, SM background events can only survive after all event selections when $ m_a = 10 $ GeV, resulting in a cusp point at $ m_a = 10 $ GeV inFig. 5.

Figure 5.(color online) Future bounds on $ c^A_\mu/\Lambda $ of GeV-scaleμALPs in theEWVscenario from the muon collider with $ {\cal L} = 120 $ fb $ ^{-1} $ and $ {\cal L} = 840 $ fb $ ^{-1} $ within a $ 95$ % CL or 10 survival events for background-free cases (dotted lines for theμALP prompt decay and dashed lines for theμALP as an LLP) as well as existing bounds (bulk regions). Here, we label " $ 2\mu $ " and " $ J_{\mu} $ " to identify two types of signatures at the muon collider. $ {\cal B}_W\lesssim 10^{-5} $ represents $ {\cal B} (W^{\pm}\rightarrow \mu^{\pm}\nu_{\mu}a) < 10^{-5} $ [50] (blue bulk). For lightμALPs, BaBar [49] (hotpink bulk), $ B^{\pm}\rightarrow K^{\pm}a\rightarrow K^{\pm}(\gamma\gamma) $ [80] (lightskyblue bulk) are considered. Some other collider bounds are in order: ATLAS $ 2\gamma $ [32,33,35] (magenta bulk), ATLAS $ 3\gamma $ [34,35] (red bulk), OPAL $ 3\gamma $ [30,35] (green bulk), and ATLAS/CMS (PbPb) [81] (yellow bulk). Finally, the bound from $ (g-2)_{\mu} $ [82] is labeled as orange bulk.

$ m_a $ /GeV	$ \lvert\eta_{\mu_{1,2}}\rvert $	$ \lvert\eta_{{\:/ E}}\rvert $	$ {\:/ E}/M_{\mu_1\mu_2} $	$ \Delta M_{\mu_1\mu_2} $	$ \Delta\phi_{\mu_2, {\not E}} $
$ 10 $	$< 3.0 $	$< 3.0 $	$> 140 $	same	$ (2.5,3.6) $
$ 20 $	$< 2.0 $	$< 2.0 $	$> 70 $	same	$ (2.7,3.6) $
$ 30 $	$< 1.8 $	$<1.9 $	$> 50 $	same	same
$ 40 $	$< 1.6 $	same	$> 40 $	same	same
$ 60 $	same	same	$> 22 $	$< 1.5 $	same
$ 70 $	same	same	$> 22 $	$< 1.5 $	same
$ 80 $	same	same	$> 20 $	$< 2.3 $	$ (2.8,3.5) $
$ 85 $	same	same	$> 18 $	$< 2.3 $	$ (2.8,3.5) $

Table 8.Changes in several event selections in the $ {EWV}$ scenario with some $ m_a $ benchmark points for two isolated muons plus $ {\not E} $ at a muon collider, where $ \Delta M_{\mu_1\mu_2}\equiv \lvert M_{\mu_1\mu_2}-m_a\rvert $ , and "same" means the same event selection as the benchmark point $ m_a = 50 $ GeV in the main text.

ForμALPs as LLPs, we first consider the physical size of the radius of the proposed detectors for muon colliders [83]. Several relevant detector parameters for the inner and outer radii of the vertex detector, ECAL, HCAL, and muon system are summarized as follows: (1) $ 3.0 \leq R_{\text{vertex}} \leq 10.4 $ cm, (2) $ 150.0 \leq R_{\text{ECAL}} \leq 170.2 $ cm, (3) $ 174.0 \leq R_{\text{HCAL}} \leq 333.0 $ cm, and (4) $ 446.1 \leq R_{\text{muon}} \leq 645.0 $ cm. Therefore, we simply consider theμALP lab frame decay length within $ 10^{-3} \leq\gamma\beta c \tau_a\leq 6.4 $ m as a detectable LLP with a muon pair displaced vertex and displaced $ J_{\mu} $ signatures at a muon collider. We assume that both a muon pair displaced vertex and displaced $ J_{\mu} $ signatures at muon colliders are background-free after the trigger andμALP mass window selection implementation, as described previously. The analysis of a muon pair displaced vertex and displaced $ J_{\mu} $ signatures at the muon collider is performed using an integrated luminosity of $ {\cal L} = 840 $ fb $ ^{-1} $ . The results of the LLP study to search the $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a \overline{\nu}_{\mu} $ channel in the EWV scenario are summarized inFig. 5with dashed lines. The signal efficiency of the two isolated muon signatures decreases when $ m_a\lesssim 30 $ GeV because these two muons become too close to each other and cannot pass the muon isolation criterion. Similarly, grouping two muons inside a $ J_{\mu} $ candidate is challenging for $ m_a\gtrsim 10 $ GeV at the muon collider. Therefore, the analysis of signatures with two isolated muons and $ J_{\mu} $ is ideal forμALP searches in the middle $ m_a $ range.

Several existing bounds are also shown inFig. 5for comparison. First, according to the interaction in Eq. (1), there is a newWboson exotic decay channel, $ W^{+}\rightarrow \mu^{+}\nu_{\mu}a $ . The precision measurements ofWboson width ( $ \Gamma_W = 2.085\pm 0.042 $ GeV [84]) can indirectly testμALPs with $ m_a < M_W $ in the EWV scenario. Here, we conservatively require the branching ratio of $ W^{+}\rightarrow \mu^{+}\nu_{\mu}a $ to be less than $ 10^{-5} $ [50] and mark it as the blue bulk inFig. 5. For lighterμALPs ( $ m_a\lesssim 5 $ GeV), searching for four muons in the final state [80] (hotpink bulk) and $ B^{\pm}\rightarrow K^{\pm}a $ (light skyblue bulk) via BaBar experiments can already constrain some parameter space in the upper-left corner. For heavierμALPs ( $ m_a > 5 $ GeV), the ATLAS $ 2\gamma $ [32,33,35] (magenta bulk), ATLAS $ 3\gamma $ [34,35] (red bulk), OPAL $ 3\gamma $ [30,35] (green bulk), and ATLAS/CMS (PbPb) [81] (yellow bulk) can already exclude some parameter space with $ c^A_{\mu} / \Lambda\gtrsim 10^2 $ TeV $ ^{-1} $ . However, the precision measurements of muon magnetic moment can also provide constraints forμALPs. The combined measurement from Fermilab and Brookhaven is reported as $ a^{\text{EXP}}_{\mu}=116,592,061(41)\times 10^{-11} $ [85], and if we consider the lattice calculation for hadronic vacuum polarization (HVP), the SM prediction value changes to $ a^{\text{SM}}_{\mu}=116,591,954(55)\times 10^{-11} $ [86]. In this situation, the deviation of $ (g-2)_{\mu} $ is reported as $ \Delta a_{\mu} = a^{\text{EXP}}_{\mu} -a^{\text{SM}}_{\mu} = 107(69)\times 10^{-11} $ , and we consider the $ \Delta a_{\mu} $ observation within $ 2\sigma $ forμALPs in this study. The one-loop contributions from lightμALPs to $ (g-2)_{\mu} $ are negative and can be written as ²

where the first term originates from theμ-a-μloop, and the second term arises from theμ-a-γloop, as shown in Ref. [82], in the following form:

Here, $ r=m^2_{\mu}/m^2_a $ , $ \Delta^2 = m^2_{\mu}x^2 +m^2_a (1-x) $ , and $ \Lambda_{\rm loop} $ is the cut-off scale of the loop integration, which is taken to be $ 1 $ TeV here. The strongest constraint among all of the above is from $ (g-2)_{\mu} $ [82] (orange bulk), with $ c^A_{\mu} / \Lambda\gtrsim 10 $ TeV $ ^{-1} $ and extending to a wide range of $ m_a $ . It is important to note that all of the above bounds are rescaled according to our definition of ALP-muon interactions in Eq. (2) andμALP decay branching ratios inFig. 1. However, several other bounds, such as OPAL $ 2\gamma $ [30,35], Belle II [87], and LHCb [88], are so weak that we do not include them here. In comparison to these existing bounds, our proposal to search forμALPs via $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ at muon colliders remains attractive. Furthermore, the possible future bounds of $ c^A_\mu/\Lambda $ can reach less than $ 0.01-0.1 $ TeV $ ^{-1} $ , which opens new doors to explore $ m_a $ in the EWV scenario below the electroweak scale.

In addition, as discussed in Sec. III, in the EWV scenario, cross sections are more than six orders of magnitude larger than those in the EWP scenario for $ \mu^{+}\mu^{-}\rightarrow \nu_{\mu} a\overline{\nu_{\mu}} $ processes inFig. 3. Therefore, future bounds from this channel in the EWP scenario are as small as the existing bounds. Additionally, almost the entire cross-section originates from the $ \boldsymbol{aVV'} $ interaction in the EWP scenario at the muon collider. Comparing the EWP scenario with the EWV scenario, the longitudinal momentum ( $ P_z $ ) becomes larger than the transverse momentum ( $ P_T $ ) for a pair of two isolated muons because the dominant contribution in the signal process is $ \mu^+\mu^-\rightarrow Za\rightarrow (\nu\overline{\nu})(\mu^+\mu^-) $ instead of that from the four-point interaction. When $ m_a \gtrsim 30 $ GeV, the total energy is roughly equally divided intoZand the ALP, resulting in large changes in some kinematic distributions. For comparison, we use the same event selections for both the EWV and EWP scenarios. Most of the signal efficiencies are below $ 10$ % in the EWP scenario because the condition $ P_T^{\mu_1} > 200 $ GeV is too stringent in this situation. Meanwhile, the efficiency of the signal is also highly sensitive to $ {\not E}/M_{\mu_1\mu_2} $ . Thus, as we can expect, the distributions of two isolated muons in the EWP scenario are distinct from those in the EWV scenario. Eventually, the EWP signal efficiency is approximately $ 10$ % to $ 40$ % less than that of the EWV efficiency. At the same time, we explore the potential results from the search forμALPs in different channels with $ {\cal L} = 1000 $ fb $ ^{-1} $ . Except for $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , we also include $ \mu^{+}\mu^{-}\rightarrow \gamma a $ and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ in the EWP scenario. Owing to different generation mechanisms among these channels, the coverage range of the interval of $ m_a $ may vary. The case of the prompt decay ofμALPs with two isolated muons in the $ \mu^{+}\mu^{-}\rightarrow \gamma a $ channel yields a background-free scenario when $ m_a\gtrsim 10 $ GeV, causing a cusp point at $ m_a = 10 $ GeV inFig. 6. The case of the prompt decay ofμALPs with four isolated muons in the $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ channel causes a cusp point at $ m_a = 5 $ GeV inFig. 6. Finally, we find that a photon plus aμALP channel has the best potential when searching forμALPs in the EWP scenario. The possible future bounds on $ c^A_\mu/\Lambda $ can reach values lower than $ 1-10 $ TeV $ ^{-1} $ , which is only slightly greater than existing bounds.

Figure 6.(color online) Future bounds on $ c^A_\mu/\Lambda $ of GeV-scaleμALPs from the muon collider in theEWPscenario with $ {\cal L} = 1000 $ fb $ ^{-1} $ within a $ 95$ % CL or 10 survival events for background-free cases. Only the bounds from theμALP prompt decay are considered (dashed lines for $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , solid lines for $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , and dotted lines for $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ ). The labels " $ 2\mu $ " and " $ J_{\mu} $ " are the same as inFig. 5.

ALPs are pseudo-Nambu Goldstone bosons that exist beyond the SM. In the effective field theory framework, ALPs are allowed to have masses ranging from nearly massless to the electroweak scale or higher, and their couplings with SM particles can be investigated independently. Therefore, it is crucial to search for ALPs with various mass ranges and interaction types. This study focuses on the search for GeV-scaleμALPs, a specific type of ALP that interacts predominantly with muons, at a muon collider.

Producing GeV-scaleμALPs is challenging owing to their suppressed production cross sections, which are proportional to the square of the muon mass. Hence, a new proposal is necessary to produce them effectively at high-energy colliders. This study proposes four production channels that can be used to search forμALPs at muon colliders, $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ , $ \mu^{+}\mu^{-}\rightarrow \gamma a $ , $ \mu^{+}\mu^{-}\rightarrow Z a $ , and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ , which rely on a four-point interaction,W-μ- $ \nu_{\mu} $ -a, and/or interactions arising from the chiral anomaly, which do not depend on the muon mass. It is noteworthy that in the EWV scenario, the cross section of the $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ process is six to seven orders of magnitude larger than that of other channels, as shown inTable 8, owing to the energy enhancement behavior resulting from theW-μ- $ \nu_{\mu} $ -ainteraction. However, in the EWP scenario, the four-point interaction disappears, and the $ \mu^{+}\mu^{-}\rightarrow \gamma a $ channel has the largest cross section.

In the search for GeV-scaleμALPs at a muon collider, different search strategies are employed for the EWV and EWP scenarios. We focus on the production channels $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ in the EWV scenario and $ \mu^{+}\mu^{-}\rightarrow \gamma a $ and $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ in the EWP scenario. However, GeV-scaleμALPs mainly decay into a pair of muons. When a lightμALP is highly boosted and produced at a muon collider, these two muons are too collimated to pass standard muon isolation criteria and form a novel object called a muon-jet, $ J_{\mu} $ . Therefore, this study explores two types of signatures: (1) two isolated muons plus other parts, and (2) $ J_{\mu} $ plus other parts. These two signature types are complementary in the search for GeV-scaleμALPs. The signature of $ J_{\mu} $ can cover the low-massμALP detection range well, and the signature of two isolated muons can cover the high-massμALP detection range. After a comprehensive signal-to-background analysis for these two types of signatures at a muon collider, future bounds for $ c^A_\mu/\Lambda $ are shown to be more than three orders of magnitude stronger than existing bounds forμALPs with $ 1 $ GeV $ \leq m_a\lesssim M_W $ at an integrated luminosity of $ {\cal L} = 120 $ fb $ ^{-1} $ for promptμALP decay and $ {\cal L} = 840 $ fb $ ^{-1} $ for theμALP as a long-lived particle in the EWV scenario, as illustrated inFig. 5. However, future bounds for $ c^A_\mu/\Lambda $ ofμALPs with $ 1 $ GeV $ \leq m_a\lesssim M_W $ are shown to barely exceed existing bounds in the EWP scenario, even with an integrated luminosity of $ {\cal L} = 1000 $ fb $ ^{-1} $ , as illustrated inFig. 6. Overall, this study provides important insights into the potential of exploring GeV-scaleμALPs. Such efforts will motivate experimentalists to pursueμALP searches at future muon colliders.

In this appendix, we choose several representative kinematic distributions for both signals and backgrounds at a muon collider:

● For the signature of two isolated muons plus $\not E $ at a $ \mu^+\mu^- $ collider, $ P_T^{\mu_1} $ , $ \eta_{\mu_1} $ , $ \not E $ , $ M_{\mu_1\mu_2} $ , $ \Delta\phi_{\mu_1,\not E} $ , and $ \Delta\phi_{\mu_2,\not E} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 0.1 $ TeV $ ^{-1} $ are shown inFig. A1.

Figure A1.(color online) Several signal and background kinematic distributions for the signature of two isolated muons plus $ {\:/ E} $ at a $ \mu^+\mu^- $ collider. $ P_T^{\mu_1} $ , $ \eta_{\mu_1} $ , $ {\:/ E} $ , $ M_{\mu_1\mu_2} $ , $ \Delta\phi_{\mu_1,{\:/ E}} $ , and $ \Delta\phi_{\mu_2,{\:/ E}} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 0.1 $ TeV $ ^{-1} $ .

● For the signature of $ J_{\mu} $ plus $ \not E $ at a $ \mu^+\mu^- $ collider, $ P_T^{J_{\mu}} $ , $ \eta_{J_{\mu}} $ , $ \not E $ , and $ M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 0.1 $ TeV $ ^{-1} $ are shown inFig. A2.

Figure A2.(color online) Several signal and background kinematic distributions for the signature of $ J_{\mu} $ plus $\not E$ at a $ \mu^+\mu^- $ collider. $ P_T^{J_{\mu}} $ , $ \eta_{J_{\mu}} $ , $\not E$ , and $ M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 0.1 $ TeV $ ^{-1} $ .

● For the signature of two isolated muons plusγat $ \mu^+\mu^- $ colliders, $ P_T^{\mu_1} $ , $ \eta_{\mu_1} $ , $ E_{\gamma} $ , $ M_{\mu_1\mu_2} $ , $ \Delta\phi_{\mu_1,\gamma} $ , and $ E_{\gamma}/M_{\mu_1\mu_2} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ are shown inFig. A3.

Figure A3.(color online) Several signal and background kinematic distributions for the signature of two isolated muons plusγat a $ \mu^+\mu^- $ collider. $ P_T^{\mu_1} $ , $ \eta_{\mu_1} $ , $ E_{\gamma} $ , $ M_{\mu_1\mu_2} $ , $ \Delta\phi_{\mu_1,\gamma} $ and $ E_{\gamma}/M_{\mu_1\mu_2} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ .

● For the signature of $ J_{\mu} $ plusγat a $ \mu^+\mu^- $ collider, $ P_T^{J_{\mu}} $ , $ \eta_{J_{\mu}} $ , $ E_{\gamma} $ , $ \eta_{\gamma} $ , $ M_{J_{\mu}} $ , and $ E_{\gamma}/M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ are shown inFig. A4.

Figure A4.(color online) Several signal and background kinematic distributions for the signature of $ J_{\mu} $ plusγat a $ \mu^+\mu^- $ collider. $ P_T^{J_{\mu}} $ , $ \eta_{J_{\mu}} $ , $ E_{\gamma} $ , $ \eta_{\gamma} $ , $ M_{J_{\mu}} $ , and $ E_{\gamma}/M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ .

● For the signature of four isolated muons at a $ \mu^+\mu^- $ collider, $ P_T^{\mu_3} $ , $ \eta_{\mu_2} $ , $ \eta_{\mu_4} \Delta\phi_{\mu_3,\mu_4} $ , $ P_T^{\mu_4}/M_{\mu_{1}\mu_4} $ , and $ M_{\mu_{1}\mu_{4}} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ are shown inFig. A5.

Figure A5.(color online) Several signal and background kinematic distributions for the signature of four isolated muons at a $ \mu^+\mu^- $ collider. $ P_T^{\mu_3} $ , $ \eta_{\mu_2} $ , $ \eta_{\mu_4} \Delta\phi_{\mu_3\mu_4} $ , $ P_T^{\mu_4}/M_{\mu_{1}\mu_4} $ , and $ M_{\mu_{1}\mu_{4}} $ distributions for $ m_a = 50 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ .

● For the signature of two isolated muons plus $ J_{\mu} $ at a $ \mu^+\mu^- $ collider, $ \eta_{\mu_2} $ , $ \eta_{J_{\mu}} $ , $ P_T^{J_{\mu}} $ , $ \Delta \phi_{\mu_{2}\mu_3} $ , $ M_{\mu_{1}\mu_4}/(P_T^{\mu_1}+P_T^{\mu_4}) $ , and $ M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ are shown inFig. A6.

Figure A6.(color online) Several signal and background kinematic distributions for the signature of $ J_{\mu} $ plus two isolated muons at a $ \mu^+\mu^- $ collider. $ \eta_{\mu_2} $ , $ \eta_{J_{\mu}} $ , $ P_T^{J_{\mu}} $ , $ \Delta \phi_{\mu_{2}\mu_3} $ , $ M_{\mu_{1}\mu_4}/(P_T^{\mu_1}+P_T^{\mu_4}) $ , and $ M_{J_{\mu}} $ distributions for $ m_a = 5 $ GeV with $ c^A_\mu / \Lambda = 10 $ TeV $ ^{-1} $ .

Conversely, we modify the event selections for the detection of two isolated muons plus $ \not E $ with varying $ m_a $ at a muon collider (as shown inTable 8) to optimize the signal efficiency. Specifically, we adjust the ranges of $ \eta_{\mu_{1,2}} $ and $ \eta_{\not E} $ for small values of $ m_a $ because loosening these criteria can improve signal detection while still eliminating all background events with the current selection criteria. In addition, we adjust $ \not E/M_{\mu_1\mu_2} $ based on signal and background distributions because it decreases as $ m_a $ increases. Conversely, we do not optimize event selections when detecting $ J_{\mu} $ plus $ \not E $ with varying $ m_a $ at a muon collider since the relevant backgrounds are already unlikely to satisfy the conditions of two detectable muons in the muon spectrometer and forming an energetic $ J_{\mu} $ in the central region. As a result, nearly all of these signals are free from background events after the cut-(3) selection inTable 3. Similarly, we fine-tune event selections for the detection of two isolated muons plus a photon and four isolated muons at a muon collider for different $ m_a $ in the EWP scenario, as listed inTables A1andA2.

$ m_a $ /GeV	$ P_T^{\mu_2} $	$ \lvert\eta_{\mu_{1,2}}\rvert $	$ \lvert\eta_{E_{\gamma}}\rvert $	$ E_{\gamma}/M_{\mu_1\mu_2} $	$ \Delta M_{\mu_1\mu_2} $
$ 30 $	$> 300 $	same	$< 1.8 $	$> 46 $	same
$ 40 $	$> 400 $	same	same	$> 35 $	$< 1.8 $
$ 60 $	same	$< 1.2 $	$< 1.2 $	$> 24 $	$< 2.2 $
$ 70 $	same	$< 1.2 $	$< 1.2 $	$ (20,23) $	$< 2.8 $
$ 80 $	same	$< 1.2 $	$< 1.2 $	$ (17,20) $	$< 3.0 $
$ 85 $	same	$< 1.2 $	$< 1.2 $	$ (16,19) $	$< 3.0 $

Table A1.Similar toTable 8, but for the $ \mu^{+}\mu^{-}\rightarrow \gamma a $ channel.

$ m_a $ /GeV	$ \lvert\eta_{\mu_{1}}\rvert $	$ \lvert\eta_{\mu_{4}}\rvert $	$ \Delta \phi_{\mu_{3},\mu_4} $	$ P_{T}^{\mu_4}/M_{\mu_{1}\mu_{4}} $	$ \Delta M_{\mu_1\mu_4} $
$ 30 $	same	$<1.6 $	$>0.8 $	same	$<4.0 $
$ 40 $	same	same	same	same	$<4.0 $
$ 60 $	same	same	same	$>0.06 $	same
$ 70 $	$<1.8 $	$< 1.4 $	$>1.2 $	$>0.06 $	same
$ 80 $	$<1.8 $	$< 1.4 $	$>1.2 $	$>0.08 $	same
$ 85 $	$<1.8 $	$< 1.4 $	$>1.2 $	$>0.08 $	$< 4.0 $

Table A2.Similar toTable 8, but for the $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ channel.