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 the W boson, with the novel energy enhancement effect. Similarly, this type of W-$ \ell $-ν-a interaction with the energy enhancement effect has also been proposed as a promising approach to the search for leptophilic ALPs via t-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 no W-μ-$ \nu_{\mu} $-a interaction in the electroweak preserving (EWP) scenario, and therefore $ \mu^{+}\mu^{-}\rightarrow \gamma a $ yields the largest cross section. In the EWV scenario, we discover that the channel $ \mu^+\mu^-\rightarrow \nu_{\mu}a\overline{\nu_{\mu}} $ with the W-$ \ell $-ν-a interaction 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 the EWP scenario 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} $ and Z and 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 the W boson and charged mesons [50], as well as in searching for heavier μALPs in t-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 the W boson is highly suppressed and can be safely neglected.

Figure 1 shows 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 in Fig. 2 and 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. 3 shows 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 the EWV scenario, 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 in Fig. 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 in Fig. 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 the t-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 in Fig. 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 in Table 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 the W-μ-$ \nu_{\mu} $-a interaction. 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 to Fig. 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 of Fig. 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 a b-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 in Table 2, and several kinematic distributions are shown in Fig. A1 of 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. In Fig. 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 the b jet 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 significance Z [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 in Table 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 of Fig. 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 from c, b jets 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 in Table 3, and several kinematic distributions are shown in Fig. A2 of 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 in Fig. 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 in Table 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 the i-th muon and the photon. The cut-flow table detailing the signal and background for each event selection is presented in Table 4, with several relevant kinematic distributions shown in Fig. A3 of Appendix A. First, we find that two isolated muons in the signal are predominantly located in the relatively low transverse momentum regions (as shown in Fig. 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 in Fig. 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 in Table 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 from c and b jets 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 in Table 5, with several relevant kinematic distributions shown in Fig. A4 of Appendix A. Employing all of the event selections listed in Table 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 to Table 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 in Table 6, and several kinematic distributions are shown in Fig. A5 of 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 to Table 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 in Table 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 in Table 7, and the corresponding kinematic distributions can be found in Fig. A6 of Appendix A. Upon applying all of the event selections outlined in Table 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 to Table 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 in Fig. 5 to 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 in Table 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 in Fig. 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 in Fig. 5.

Figure 5.  (color online) Future bounds on $ c^A_\mu/\Lambda $ of GeV-scale μALPs in the EWV scenario 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 in Fig. 5 with 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 in Fig. 5 for comparison. First, according to the interaction in Eq. (1), there is a new W boson exotic decay channel, $ W^{+}\rightarrow \mu^{+}\nu_{\mu}a $. The precision measurements of W boson 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 in Fig. 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 in Fig. 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 in Fig. 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 into Z and 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 in Fig. 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 in Fig. 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 the EWP scenario 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 in Fig. 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 in Table 8, owing to the energy enhancement behavior resulting from the W-μ-$ \nu_{\mu} $-a interaction. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Fig. 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 in Table 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 in Table 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 in Tables A1 and A2.

$ 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 to Table 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 to Table 8, but for the $ \mu^{+}\mu^{-}\rightarrow \mu^{+}\mu^{-} a $ channel.