Long-lived sterile neutrinos from axionlike particles at the Super Tau–Charm Facility

In recent years, the search for new physics in the form of long-lived particles (LLPs) has become an important direction in particle physics [1−6]. This development is motivated partly by the absence of discoveries of heavy, promptly decaying new fundamental particles and partly by the rich motivations often associated with models predicting such LLPs, including explanations of non-vanishing neutrino masses or dark matter.

Among the many possible theoretical scenarios in which LLPs naturally arise, a class of portal-physics models is often studied, perhaps because these models serve as representative toy models of ultraviolet completions and predict rich phenomenology at various terrestrial experiments. Such portals include heavy neutral leptons (HNLs) [7−9] (see also Ref. [10] for a review), a pseudoscalar boson such as an axion or axionlike particle (ALP) [11−17], a dark scalar [18−23], and a dark photon [24−29]. Examples of scenarios involving more than one such portal include the dark axion portal [30], which predicts both an ALP and a dark photon; the hidden Abelian Higgs model [19, 31], which contains both an additional scalar and a dark photon; and an ALP-HNL portal [32−34].

HNLs are hypothetical spin-$ \dfrac{1}{2} $ particles that are singlets under the Standard-Model (SM) gauge group. They mix with active neutrinos and are highly motivated as an explanation of nonzero neutrino masses through the type-I seesaw mechanism [35−39] or various low-scale seesaw mechanisms [40−43], and as dark-matter candidates [44−51] if their masses are in the keV range. On the other hand, the QCD axion was originally proposed as a pseudo-Nambu-Goldstone boson arising from the spontaneous breaking of a global U(1)_PQ symmetry to solve the strong CP problem [11]. More generally, ALPs appear ubiquitously in extensions of the SM, where their masses and couplings are largely independent parameters. They are closely connected to the strong CP problem [11], dark matter [52−55], and the hierarchy problem [56]. UV completions that predict ALPs include string compactifications [57, 58], Froggatt-Nielsen models of flavor [59, 60], and supersymmetry models [61].

In this work, we focus on the ALP-HNL portal and perform a phenomenological study. In the ALP-HNL portal model, ALPs and HNLs can couple to each other, leading to phenomenology more complex than that arising from simple, single-portal scenarios. Previous phenomenological analyses of this model scenario have studied ALPs produced in heavy-scalar decays [34], electroweak processes [62], gluon fusion [62−64], meson mixing [65], and meson decays [66, 67], in collider or beam-dump experiment settings where the ALP further decays into a pair of HNLs; they have also studied the same new particles in the contexts of cosmology [68, 69] and dark matter [70].

In particular, for meson decays, one may assume tree-level quark-flavor-violating (QFV) couplings that lead to ALP production in association with another meson. Such QFV couplings are predicted in UV completions such as DFSZ- and KSVZ-type axion models [71−78]; see also [79−84] for additional possibilities. The quark-flavor-diagonal couplings can be suppressed in, for example, the astrophobic axion model [77] and several Froggatt-Nielsen model variants [72, 75, 85]. In this phenomenological analysis, we remain agnostic about the origin of the particular quark-flavor structure of the ALP that we assume and treat the ALP QFV couplings as independent parameters. For previous phenomenological studies in this direction, see, e.g., Refs. [81, 86−92].

We study ALPs that are produced in rare decays of charm mesons and subsequently decay promptly into a pair of long-lived HNLs, and estimate the search prospects for such long-lived HNLs at the proposed Super Tau–Charm Facility (STCF) [93, 94] in Hefei, China. STCF is planned to collide an electron beam and a positron beam at center-of-mass (COM) energies $ \sqrt{s} $ between 2 GeV and 7 GeV with symmetric beam energies. We restrict our analysis to the mode with $ \sqrt{s}=3.773 $ GeV, leading to resonant production of $ D^+ D^- $ pairs. The ALPs can originate from the decays of charged D mesons in association with a charged pion and then decay promptly into a pair of long-lived HNLs. To estimate the reconstruction efficiencies of displaced vertices (DVs) arising from the decays of the long-lived HNLs, we perform both fast Monte Carlo (MC) simulations at the truth level and dedicated MC simulations at the detector level with the OSCAR framework [95−98], which is based on GEANT4 [99, 100] and dedicated to realistic STCF-detector simulations. Taking into account the present bounds on the model parameters, we determine the sensitivity reach of STCF to long-lived HNLs in terms of the mass and mixing parameters.

This paper is organized as follows. We first introduce our model setup in Sec. II and then provide the details of our truth-level analysis, along with a background discussion, in Sec. III. In Sec. IV, we explain the detector-level analysis within the OSCAR framework. The numerical results are provided in Sec. V, and we summarize the paper in Sec. VI.

We work within an effective framework for an axion-like particle, defined as follows [34, 62, 70, 81, 88, 90, 101, 102]:

where a is a CP-odd ALP with mass $ m_a $, N denotes the HNL, Λ is the effective cutoff scale, and q labels the quark fields $ u_L $, $ u_R $, $ d_L $, and $ d_R $, with generation indices $ i,j=1,2,3 $.

In Eq. (1), the couplings $ g^{q}_{i,j} $ and $ g_N $ are dimensionless. Meson decays $ P\to P'/V a $, where P and $ P' $ denote pseudoscalar mesons and V denotes a vector meson, can be induced by a non-vanishing coupling $ g^q_{i,j} $ with $ i\neq j $. In this work, we restrict our analysis to the flavor indices $ (2,1) $ in the up-type quark sector. With the coupling $ g^{u+}_{2,1}=g_{2,1}^{u_R}+g_{2,1}^{u_L} $ ($ g^{u-}_{2,1}=g_{2,1}^{u_R}-g_{2,1}^{u_L} $), charm mesons decay into an ALP and a pseudoscalar (vector) meson. Specifically, a nonzero $ g^{u+}_{2,1} $ coupling can induce the charm-meson decays $ D^0\to \pi^0/\eta/\eta'\,a $, $ D^+\to \pi^+ \,a $, and $ D_s^+\to K^+ \, a $. By contrast, the $ g^{u-}_{2,1} $ coupling can lead to $ D^0\to \rho^0/\omega \,a $, $ D^+\to \rho^+ \,a $, and $ D_s^+\to K^{*+} \,a $. We assume that either $ g_{2,1}^{u_L} $ or $ g_{2,1}^{u_R} $ vanishes, so that both $ P\to P' $ and $ P\to V $ transitions occur through a single coupling; in the following, we denote the single non-vanishing ALP-quark coupling by $ g^u_{2,1} $. For the numerical analysis, we focus on the decay $ D^+\to \pi^+ a $ and its charge-conjugate channel, as our theoretical evaluations indicate that these channels provide the dominant decay mode for ALP production. The sensitivities can be enhanced by including the other channels.

In the remainder of the paper, we denote the coupling that induces the decay $ D^+\to \pi^+ a $ by $ g^u_{2,1} $. All other ALP couplings to quarks are assumed to be zero.

The decay rate of $ D^+\to \pi^+ a $ can be computed using the following formula [81]:

where $ \lambda(x,y)\equiv 1+x^2+y^2-2\,x -2\,y -2xy $ and $ m_{D^+/\pi^+} $ denotes the mass of the $ D^+/\pi^+ $ meson. $ F_0^{D^+\to \pi^+} $ is the $ D\to \pi $ transition form factor [103], which we extract from Ref. [104]. To calculate the branching ratio of $ D^+ \to \pi^+ a $, Br($ D^+ \to \pi^+ a $), we divide $ \Gamma\left(D^+ \to \pi^+ a\right) $ by the measured total decay width of the $ D^+ $ meson [105].

We derive the current bound on $ g^u_{2,1}/\Lambda $ from the existing limits on the decay branching ratio of $ D^0\to \pi^0 \nu\bar{\nu} $, since in our theoretical setup the ALP decays exclusively into a pair of long-lived HNLs, which are typically unobserved in the main detector. For $ D^0\to \pi^0 \nu\bar{\nu} $, which is a three-body decay, we assume that the kinematics do not differ substantially from those of the two-body decay $ D^0\to \pi^0 a $ that would arise from the $ g^u_{2,1} $ coupling. Using Eq. (2), we derive $ g^u_{2,1}/\Lambda\lesssim 2\times 10^{-4} $ TeV⁻¹ (see also Ref. [90]) based on the current upper bound obtained by BESIII [106]. In our numerical analysis, we restrict ourselves to two representative values of $ m_a $, 1.4 GeV and 1.6 GeV, and set $ g^u_{2,1}/\Lambda $ to its upper bound, $ 2\times 10^{-4} $ TeV⁻¹. Since we study HNL pairs produced in ALP decays, the HNL mass can therefore be at most 0.8 GeV in this work.

In the $ m_a $ range of interest, the ALP decays dominantly into a pair of HNLs with the decay width

where $ m_N $ is the HNL mass.

We emphasize that, in this work, we assume that the HNL N is Majorana in nature, that only one HNL is kinematically relevant, and that it mixes exclusively with the electron neutrino.

For the numerical analysis, we fix $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹ as a benchmark value, satisfying the perturbativity requirement $ g_N m_N/\Lambda<1 $. Thus, in the $ m_a $ and $ m_N $ ranges of interest, the ALP decays promptly unless $ m_N\approx 0 $ or $ m_N\approx m_a/2 $ [66, 67]. If $ g_N/\Lambda $ is orders of magnitude smaller, the ALP would become long-lived, thereby weakening the sensitivities; see Ref. [66], where this effect was studied. In addition, for a different value of $ g_N/\Lambda $, our results remain unchanged as long as the ALP decays promptly and exclusively into a pair of long-lived Majorana HNLs.

The HNL N participates in the SM weak interactions through mixing with the active neutrinos

where g denotes the SU(2) gauge coupling, $ V_{\alpha N} $ is the mixing parameter between N and the active neutrino $ \nu_\alpha $, with $ \alpha=e,\mu,\tau $, $ \ell $ is an SM charged lepton, $ P_L $ is the usual left-chiral projector, W and Z are the W- and Z-bosons, $ \theta_W $ denotes the weak mixing angle, and $ V^L_{\alpha i} $ is the neutrino mixing matrix in the left-handed sector between the flavor eigenstate $ \nu_\alpha $ and the mass eigenstate $ \nu_i $, with $ i=1,2,3 $.

We consider only two parameters within the pure HNL sector: $ m_N $, the mass of N, and $ V_{e N} $, its mixing strength with the electron neutrino. In the $ m_N $ range of interest, HNL decays can be well described by two-body decays alone. Thus, for the HNL decay, we need only compute the two-body decay widths at tree level, following Ref. [107].

We focus on the range $ 0.2\; {\rm{GeV}} \leq m_N \leq 0.8\; {\rm{GeV}} $. In this mass range, the existing bounds are dominated by beam-dump experiments: PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114]. These bounds were obtained for HNLs in minimal scenarios in which the HNLs do not couple to any beyond-the-Standard-Model (BSM) particles. However, some of these searches are applicable to our theoretical scenario, namely the CHARM [110] and BEBC [114] searches. These two searches obtained similar bounds in the mass range of interest, and we reinterpret the bounds from the CHARM search [110] in our model. Rather than performing a full recast based on MC simulations, which would be difficult, we follow the simple and rapid reinterpretation method based on rescaling proposed in Ref. [92]. We note that this approach provides approximate rather than precise results.

In Fig. 1, we show both the HNL $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $ (left) and a Feynman diagram of our signal process, $ D^+\to \pi^+ a $ followed by $ a\to N N $ (right). We note that $ c\tau $ is proportional to $ 1/|V_{eN}|^{2} $. The left panel shows that the HNL can be long-lived if $ m_N $ or $ |V_{eN}|^2 $ is small. The right panel illustrates the signal process considered here, in which the ALP is on-shell and decays promptly.

Figure 1.  (color online) Left: $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $. Right: Feynman diagram illustrating the decays $ D^+ \to \pi^+ a $ and $ a\to N N $, mediated by the couplings $ g^u_{2,1}/\Lambda $ and $ g_N/\Lambda $, respectively.

We work within an effective framework for an axion-like particle, defined as follows [34, 62, 70, 81, 88, 90, 101, 102]:

where a is a CP-odd ALP with mass $ m_a $, N denotes the HNL, Λ is the effective cutoff scale, and q labels the quark fields $ u_L $, $ u_R $, $ d_L $, and $ d_R $, with generation indices $ i,j=1,2,3 $.

In Eq. (1), the couplings $ g^{q}_{i,j} $ and $ g_N $ are dimensionless. Meson decays $ P\to P'/V a $, where P and $ P' $ denote pseudoscalar mesons and V denotes a vector meson, can be induced by a non-vanishing coupling $ g^q_{i,j} $ with $ i\neq j $. In this work, we restrict our analysis to the flavor indices $ (2,1) $ in the up-type quark sector. With the coupling $ g^{u+}_{2,1}=g_{2,1}^{u_R}+g_{2,1}^{u_L} $ ($ g^{u-}_{2,1}=g_{2,1}^{u_R}-g_{2,1}^{u_L} $), charm mesons decay into an ALP and a pseudoscalar (vector) meson. Specifically, a nonzero $ g^{u+}_{2,1} $ coupling can induce the charm-meson decays $ D^0\to \pi^0/\eta/\eta'\,a $, $ D^+\to \pi^+ \,a $, and $ D_s^+\to K^+ \, a $. By contrast, the $ g^{u-}_{2,1} $ coupling can lead to $ D^0\to \rho^0/\omega \,a $, $ D^+\to \rho^+ \,a $, and $ D_s^+\to K^{*+} \,a $. We assume that either $ g_{2,1}^{u_L} $ or $ g_{2,1}^{u_R} $ vanishes, so that both $ P\to P' $ and $ P\to V $ transitions occur through a single coupling; in the following, we denote the single non-vanishing ALP-quark coupling by $ g^u_{2,1} $. For the numerical analysis, we focus on the decay $ D^+\to \pi^+ a $ and its charge-conjugate channel, as our theoretical evaluations indicate that these channels provide the dominant decay mode for ALP production. The sensitivities can be enhanced by including the other channels.

In the remainder of the paper, we denote the coupling that induces the decay $ D^+\to \pi^+ a $ by $ g^u_{2,1} $. All other ALP couplings to quarks are assumed to be zero.

The decay rate of $ D^+\to \pi^+ a $ can be computed using the following formula [81]:

where $ \lambda(x,y)\equiv 1+x^2+y^2-2\,x -2\,y -2xy $ and $ m_{D^+/\pi^+} $ denotes the mass of the $ D^+/\pi^+ $ meson. $ F_0^{D^+\to \pi^+} $ is the $ D\to \pi $ transition form factor [103], which we extract from Ref. [104]. To calculate the branching ratio of $ D^+ \to \pi^+ a $, Br($ D^+ \to \pi^+ a $), we divide $ \Gamma\left(D^+ \to \pi^+ a\right) $ by the measured total decay width of the $ D^+ $ meson [105].

We derive the current bound on $ g^u_{2,1}/\Lambda $ from the existing limits on the decay branching ratio of $ D^0\to \pi^0 \nu\bar{\nu} $, since in our theoretical setup the ALP decays exclusively into a pair of long-lived HNLs, which are typically unobserved in the main detector. For $ D^0\to \pi^0 \nu\bar{\nu} $, which is a three-body decay, we assume that the kinematics do not differ substantially from those of the two-body decay $ D^0\to \pi^0 a $ that would arise from the $ g^u_{2,1} $ coupling. Using Eq. (2), we derive $ g^u_{2,1}/\Lambda\lesssim 2\times 10^{-4} $ TeV⁻¹ (see also Ref. [90]) based on the current upper bound obtained by BESIII [106]. In our numerical analysis, we restrict ourselves to two representative values of $ m_a $, 1.4 GeV and 1.6 GeV, and set $ g^u_{2,1}/\Lambda $ to its upper bound, $ 2\times 10^{-4} $ TeV⁻¹. Since we study HNL pairs produced in ALP decays, the HNL mass can therefore be at most 0.8 GeV in this work.

In the $ m_a $ range of interest, the ALP decays dominantly into a pair of HNLs with the decay width

where $ m_N $ is the HNL mass.

We emphasize that, in this work, we assume that the HNL N is Majorana in nature, that only one HNL is kinematically relevant, and that it mixes exclusively with the electron neutrino.

For the numerical analysis, we fix $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹ as a benchmark value, satisfying the perturbativity requirement $ g_N m_N/\Lambda<1 $. Thus, in the $ m_a $ and $ m_N $ ranges of interest, the ALP decays promptly unless $ m_N\approx 0 $ or $ m_N\approx m_a/2 $ [66, 67]. If $ g_N/\Lambda $ is orders of magnitude smaller, the ALP would become long-lived, thereby weakening the sensitivities; see Ref. [66], where this effect was studied. In addition, for a different value of $ g_N/\Lambda $, our results remain unchanged as long as the ALP decays promptly and exclusively into a pair of long-lived Majorana HNLs.

The HNL N participates in the SM weak interactions through mixing with the active neutrinos

where g denotes the SU(2) gauge coupling, $ V_{\alpha N} $ is the mixing parameter between N and the active neutrino $ \nu_\alpha $, with $ \alpha=e,\mu,\tau $, $ \ell $ is an SM charged lepton, $ P_L $ is the usual left-chiral projector, W and Z are the W- and Z-bosons, $ \theta_W $ denotes the weak mixing angle, and $ V^L_{\alpha i} $ is the neutrino mixing matrix in the left-handed sector between the flavor eigenstate $ \nu_\alpha $ and the mass eigenstate $ \nu_i $, with $ i=1,2,3 $.

We consider only two parameters within the pure HNL sector: $ m_N $, the mass of N, and $ V_{e N} $, its mixing strength with the electron neutrino. In the $ m_N $ range of interest, HNL decays can be well described by two-body decays alone. Thus, for the HNL decay, we need only compute the two-body decay widths at tree level, following Ref. [107].

We focus on the range $ 0.2\; {\rm{GeV}} \leq m_N \leq 0.8\; {\rm{GeV}} $. In this mass range, the existing bounds are dominated by beam-dump experiments: PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114]. These bounds were obtained for HNLs in minimal scenarios in which the HNLs do not couple to any beyond-the-Standard-Model (BSM) particles. However, some of these searches are applicable to our theoretical scenario, namely the CHARM [110] and BEBC [114] searches. These two searches obtained similar bounds in the mass range of interest, and we reinterpret the bounds from the CHARM search [110] in our model. Rather than performing a full recast based on MC simulations, which would be difficult, we follow the simple and rapid reinterpretation method based on rescaling proposed in Ref. [92]. We note that this approach provides approximate rather than precise results.

In Fig. 1, we show both the HNL $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $ (left) and a Feynman diagram of our signal process, $ D^+\to \pi^+ a $ followed by $ a\to N N $ (right). We note that $ c\tau $ is proportional to $ 1/|V_{eN}|^{2} $. The left panel shows that the HNL can be long-lived if $ m_N $ or $ |V_{eN}|^2 $ is small. The right panel illustrates the signal process considered here, in which the ALP is on-shell and decays promptly.

Figure 1.  (color online) Left: $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $. Right: Feynman diagram illustrating the decays $ D^+ \to \pi^+ a $ and $ a\to N N $, mediated by the couplings $ g^u_{2,1}/\Lambda $ and $ g_N/\Lambda $, respectively.

We work within an effective framework for an axion-like particle, defined as follows [34, 62, 70, 81, 88, 90, 101, 102]:

where a is a CP-odd ALP with mass $ m_a $, N denotes the HNL, Λ is the effective cutoff scale, and q labels the quark fields $ u_L $, $ u_R $, $ d_L $, and $ d_R $, with generation indices $ i,j=1,2,3 $.

In Eq. (1), the couplings $ g^{q}_{i,j} $ and $ g_N $ are dimensionless. Meson decays $ P\to P'/V a $, where P and $ P' $ denote pseudoscalar mesons and V denotes a vector meson, can be induced by a non-vanishing coupling $ g^q_{i,j} $ with $ i\neq j $. In this work, we restrict our analysis to the flavor indices $ (2,1) $ in the up-type quark sector. With the coupling $ g^{u+}_{2,1}=g_{2,1}^{u_R}+g_{2,1}^{u_L} $ ($ g^{u-}_{2,1}=g_{2,1}^{u_R}-g_{2,1}^{u_L} $), charm mesons decay into an ALP and a pseudoscalar (vector) meson. Specifically, a nonzero $ g^{u+}_{2,1} $ coupling can induce the charm-meson decays $ D^0\to \pi^0/\eta/\eta'\,a $, $ D^+\to \pi^+ \,a $, and $ D_s^+\to K^+ \, a $. By contrast, the $ g^{u-}_{2,1} $ coupling can lead to $ D^0\to \rho^0/\omega \,a $, $ D^+\to \rho^+ \,a $, and $ D_s^+\to K^{*+} \,a $. We assume that either $ g_{2,1}^{u_L} $ or $ g_{2,1}^{u_R} $ vanishes, so that both $ P\to P' $ and $ P\to V $ transitions occur through a single coupling; in the following, we denote the single non-vanishing ALP-quark coupling by $ g^u_{2,1} $. For the numerical analysis, we focus on the decay $ D^+\to \pi^+ a $ and its charge-conjugate channel, as our theoretical evaluations indicate that these channels provide the dominant decay mode for ALP production. The sensitivities can be enhanced by including the other channels.

In the remainder of the paper, we denote the coupling that induces the decay $ D^+\to \pi^+ a $ by $ g^u_{2,1} $. All other ALP couplings to quarks are assumed to be zero.

The decay rate of $ D^+\to \pi^+ a $ can be computed using the following formula [81]:

where $ \lambda(x,y)\equiv 1+x^2+y^2-2\,x -2\,y -2xy $ and $ m_{D^+/\pi^+} $ denotes the mass of the $ D^+/\pi^+ $ meson. $ F_0^{D^+\to \pi^+} $ is the $ D\to \pi $ transition form factor [103], which we extract from Ref. [104]. To calculate the branching ratio of $ D^+ \to \pi^+ a $, Br($ D^+ \to \pi^+ a $), we divide $ \Gamma\left(D^+ \to \pi^+ a\right) $ by the measured total decay width of the $ D^+ $ meson [105].

We derive the current bound on $ g^u_{2,1}/\Lambda $ from the existing limits on the decay branching ratio of $ D^0\to \pi^0 \nu\bar{\nu} $, since in our theoretical setup the ALP decays exclusively into a pair of long-lived HNLs, which are typically unobserved in the main detector. For $ D^0\to \pi^0 \nu\bar{\nu} $, which is a three-body decay, we assume that the kinematics do not differ substantially from those of the two-body decay $ D^0\to \pi^0 a $ that would arise from the $ g^u_{2,1} $ coupling. Using Eq. (2), we derive $ g^u_{2,1}/\Lambda\lesssim 2\times 10^{-4} $ TeV⁻¹ (see also Ref. [90]) based on the current upper bound obtained by BESIII [106]. In our numerical analysis, we restrict ourselves to two representative values of $ m_a $, 1.4 GeV and 1.6 GeV, and set $ g^u_{2,1}/\Lambda $ to its upper bound, $ 2\times 10^{-4} $ TeV⁻¹. Since we study HNL pairs produced in ALP decays, the HNL mass can therefore be at most 0.8 GeV in this work.

In the $ m_a $ range of interest, the ALP decays dominantly into a pair of HNLs with the decay width

where $ m_N $ is the HNL mass.

We emphasize that, in this work, we assume that the HNL N is Majorana in nature, that only one HNL is kinematically relevant, and that it mixes exclusively with the electron neutrino.

For the numerical analysis, we fix $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹ as a benchmark value, satisfying the perturbativity requirement $ g_N m_N/\Lambda<1 $. Thus, in the $ m_a $ and $ m_N $ ranges of interest, the ALP decays promptly unless $ m_N\approx 0 $ or $ m_N\approx m_a/2 $ [66, 67]. If $ g_N/\Lambda $ is orders of magnitude smaller, the ALP would become long-lived, thereby weakening the sensitivities; see Ref. [66], where this effect was studied. In addition, for a different value of $ g_N/\Lambda $, our results remain unchanged as long as the ALP decays promptly and exclusively into a pair of long-lived Majorana HNLs.

The HNL N participates in the SM weak interactions through mixing with the active neutrinos

where g denotes the SU(2) gauge coupling, $ V_{\alpha N} $ is the mixing parameter between N and the active neutrino $ \nu_\alpha $, with $ \alpha=e,\mu,\tau $, $ \ell $ is an SM charged lepton, $ P_L $ is the usual left-chiral projector, W and Z are the W- and Z-bosons, $ \theta_W $ denotes the weak mixing angle, and $ V^L_{\alpha i} $ is the neutrino mixing matrix in the left-handed sector between the flavor eigenstate $ \nu_\alpha $ and the mass eigenstate $ \nu_i $, with $ i=1,2,3 $.

We consider only two parameters within the pure HNL sector: $ m_N $, the mass of N, and $ V_{e N} $, its mixing strength with the electron neutrino. In the $ m_N $ range of interest, HNL decays can be well described by two-body decays alone. Thus, for the HNL decay, we need only compute the two-body decay widths at tree level, following Ref. [107].

We focus on the range $ 0.2\; {\rm{GeV}} \leq m_N \leq 0.8\; {\rm{GeV}} $. In this mass range, the existing bounds are dominated by beam-dump experiments: PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114]. These bounds were obtained for HNLs in minimal scenarios in which the HNLs do not couple to any beyond-the-Standard-Model (BSM) particles. However, some of these searches are applicable to our theoretical scenario, namely the CHARM [110] and BEBC [114] searches. These two searches obtained similar bounds in the mass range of interest, and we reinterpret the bounds from the CHARM search [110] in our model. Rather than performing a full recast based on MC simulations, which would be difficult, we follow the simple and rapid reinterpretation method based on rescaling proposed in Ref. [92]. We note that this approach provides approximate rather than precise results.

In Fig. 1, we show both the HNL $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $ (left) and a Feynman diagram of our signal process, $ D^+\to \pi^+ a $ followed by $ a\to N N $ (right). We note that $ c\tau $ is proportional to $ 1/|V_{eN}|^{2} $. The left panel shows that the HNL can be long-lived if $ m_N $ or $ |V_{eN}|^2 $ is small. The right panel illustrates the signal process considered here, in which the ALP is on-shell and decays promptly.

Figure 1.  (color online) Left: $ c\tau $ as a function of $ m_N $ for $ |V_{eN}|^2=1 $. Right: Feynman diagram illustrating the decays $ D^+ \to \pi^+ a $ and $ a\to N N $, mediated by the couplings $ g^u_{2,1}/\Lambda $ and $ g_N/\Lambda $, respectively.

The Super Tau–Charm Facility is a next-generation $ e^+e^- $ collider proposed for construction in China. It is designed to achieve a peak luminosity of $ 0.5 \times 10^{35}\;{\rm{ cm}}^{-2}{\rm{s}}^{-1} $. A key feature of STCF is its capability to perform energy scans over the COM-energy range of 2–$ 7\;{\rm{ GeV}} $, providing extensive opportunities to search for BSM physics in the tau–charm energy region. In particular, as a dedicated charm-meson factory operating at the $ \psi(3770) $ resonance, STCF is expected to produce approximately $ 2.8 \times 10^{9} $ $ D^+D^- $ meson pairs with an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $ collected in a single year of operation. This yield is approximately two orders of magnitude larger than that of BESIII [115].

STCF comprises multiple subdetectors, including the Inner Tracking system (ITK), the Main Drift Chamber (MDC), the DIRC-like Time-of-flight detector for particle identification, and the Electromagnetic Calorimeter [93]. For charged particles in the momentum range relevant to this study (as discussed below), particle identification and track reconstruction are primarily performed by the ITK and MDC. Therefore, in the phenomenological analysis, we restrict the geometric acceptance to these two subdetectors by requiring the HNLs to decay within this region; in addition, further geometric constraints are imposed to suppress background contributions. Specifically, the geometric coverage of the ITK and MDC can be approximated as $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 6\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $, where z and r denote the longitudinal and radial decay positions of the HNL, respectively. Moreover, near the interaction point (IP), prompt tracks, soft particles, and their multiple interactions with the detector material can generate numerous detector responses, making particle identification and track reconstruction challenging in this region. To suppress these potential backgrounds, we require $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 10\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $ as our fiducial volume.

Through mixing with the electron neutrino, HNLs can decay into a variety of SM final states. In general, a higher fraction of charged particles in the HNL final states leads to more precise DV reconstruction and improved invariant-mass ($ M_{\rm{inv}} $) resolution, thereby enhancing background suppression. Therefore, we restrict our analysis to channels that contain at least two charged particles and do not involve significant missing energy: (1) leptonic $ \ ll^{\prime}\nu:\ e^{\pm}\mu^{\mp}\nu_{\mu},\ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and (2) semi-leptonic $ \ e^{\pm}M^{\mp}: e^{\pm}\pi^{\mp},\ e^{\pm}\rho^{\mp},\ e^{\pm}K^{\mp} $. We plot their branching ratios as functions of $ m_N $ in Fig. 2.

Figure 2.  (color online) Decay branching ratios of N into final states containing at least two charged particles as functions of $ m_N $.

We select the $ D^- $-meson produced in the decay of $ \psi(3770) $ as the tagging side and fully reconstruct it through SM decay channels (single tag) to suppress the $ q\bar{q} $ continuum background. Its kinematic properties provide constraints on the signal process on the opposite side. On the signal side, we search for the proposed process associated with the $ D^+ $-meson. A single ALP produced on the signal side can decay into two HNLs, which may, in principle, both decay visibly within the detector acceptance. In practice, the presence of a single DV is generally sufficient for effective background suppression, and the additional background rejection provided by the single tag further strengthens this suppression. Therefore, reconstructing only one HNL on the signal side is considered sufficient to identify the event as a signal process.

After applying the selections described above, we expect the dominant remaining background to arise from photon conversion, which is difficult to suppress using the same strategy as that employed for backgrounds originating near the IP. Despite imposing a stringent geometric requirement $ r \gt 10\;{\rm{ cm}} $, together with the DV constraint and the single-tag strategy, reconstructing only one HNL on the signal side remains susceptible to contamination from photon-conversion backgrounds. Therefore, we exclude all leptonic channels, including $ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and $ e^{\pm}\mu^{\mp}\nu_{\mu} $, because they can be directly affected by the misreconstruction of lepton pairs originating from photon-conversion processes ($ \gamma \to e^+e^- / \mu^+\mu^- $).

For the STCF COM energy considered here ($ \sqrt{s}= 3.773\;{\rm{ GeV}} $), and taking into account initial-state radiation, the produced $ D^+ $- and $ D^- $-mesons have energies approximately in the range $ 1.870 $–$ 1.885\;{\rm{ GeV}} $. The momentum of the $ \pi^+ $-meson from $ D^+ $-meson decays ($ D^+ \to \pi^+ a $) is required to be at least $ 0.1\;{\rm{ GeV}} $ to penetrate a sufficient number of layers in the MDC and produce reconstructable tracks. Therefore, we require the ALP mass to be below $ 1.6\;{\rm{ GeV}} $, which further implies that the HNL mass must be less than half of the ALP mass, specifically below $ 0.8\;{\rm{ GeV}} $. Furthermore, an ALP with a mass $ m_a $ below 1.0 GeV would make the HNLs pair-produced from ALP decays kinematically incapable of decaying into $ e^{\pm}K^{\mp} $, a scenario that we intend to avoid. Therefore, we constrain the $ m_a $ range to 1.0–1.6 GeV. Specifically, the upper bound of 1.6 GeV and a slightly lower value of 1.4 GeV are selected as benchmark points. In addition, to avoid excessively soft charged particles in the HNL decay final states, which would degrade the reconstruction efficiency, we impose a lower bound on the HNL mass of $ 0.2\;{\rm{ GeV}} $. Consequently, the range of $ m_N $ is set to 0.2–0.8 GeV.

In Fig. 2, we observe that near $ 0.8\;{\rm{ GeV}} $, the opening of the $ e^{\pm}\rho^{\mp} $ channel leads to a turning point in the total branching ratios, while the branching ratio of the $ e^{\pm}\pi^{\mp} $ channel begins to decrease more rapidly. Within the HNL mass range of interest, 0.2–$ 0.8\;{\rm{ GeV}} $, the $ e^{\pm}\pi^{\mp} $ channel clearly dominates among the semi-leptonic final states. Therefore, for semi-leptonic modes, we select only the $ e^{\pm}\pi^{\mp} $ channel as the signal final state of the HNL. Accordingly, the final-state particle combination considered in this work is $ \pi^+ e^{\pm}\pi^{\mp} $. Its physical background can arise only from $ D^{\pm} \to \pi^+\pi^- e^{\pm} \nu_e $, which can be readily rejected by the displaced-vertex requirement, as discussed above, rendering this a very clean final state. We anticipate that, in a realistic detector environment, additional backgrounds, such as combinatorial backgrounds or fake-track backgrounds, may still be present; these will be discussed in the next section.

In addition, the MDC of the STCF detector adopts a superlayer wire configuration similar to those of Belle II and BESIII. Such a structure causes the detector efficiency for reconstructing charged-particle tracks and performing particle identification (PID) to gradually decrease as the particle production point moves away from the IP, with this effect being particularly pronounced for increasing radial production positions r. In this phenomenological analysis, to account for this decrease in detector efficiency, we model it using a linear function as

This approach has been used in Refs. [67, 83, 116−119]. In the next section, we present an improved, dedicated detector-level analysis to estimate the relevant effects.

In the phenomenological analysis, we use a fast MC simulation program with phase-space decays to generate truth-level kinematic information and evaluate the signal acceptance. Specifically, the kinematic properties of the HNL are generated through a cascade decay chain and subsequent boosts into the laboratory frame. The decay probability is then determined from its proper decay length ($ c\tau $) to assess whether it decays within the fiducial volume. Using this procedure, we scan the HNL mass over the range 0.2–$ 0.8\;{\rm{ GeV}} $ and its proper decay length $ c\tau $ over the range $ 10^{-6} $–$ 10^{10}\;{\rm{ m}} $. The final expression for the expected number of signal events $ N_s $, incorporating both the geometric acceptance and the linear efficiency function $ L(r) $, is given by

where

Here, $ N_{D^+D^-} $ is the total number of $ D^+D^- $ pairs produced for an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $, and $ \sum_i {\rm{Br}}^{\,i}_{\rm tag} \cdot \epsilon^{\,i}_{\rm tag} $ denotes the sum over all SM decay channels suitable for fully reconstructing the tagging-side $ D^- $ meson, with each channel weighted by its corresponding single-tag efficiency. The factor of 2 accounts for the two D mesons in each signal event. Using the single-tag $ D^+D^- $ production events reported in Ref. [93], we estimate a total tagging efficiency of approximately 19.64%. $ {\rm{Br}}(a \to NN)= $ 100%. $ N_{{\rm{MC}}} $ in Eq. (7) is the total number of MC-simulated events. The final integral in Eq. (8) simultaneously incorporates the linear efficiency function $ L(r) $ and the decay probability within the fiducial volume. $ R_r^{N_{i,j}} = $ $ \beta_r^{N_{i,j}} \gamma c\tau =\dfrac{p_T^{N_{i,j}}}{m_N} c \tau $ is the HNL transverse decay length in the laboratory frame, where $ j=1,2 $ denotes the two HNLs produced from the ALP decay. In addition, $ z_{N_{i,j}}=r \cot \theta_{N_{i,j}} $ is the longitudinal decay position of the HNL, with r denoting the transverse decay position and θ the HNL polar angle. The function $ A(z) $ is introduced to impose the longitudinal fiducial-volume constraint, defined as

The Super Tau–Charm Facility is a next-generation $ e^+e^- $ collider proposed for construction in China. It is designed to achieve a peak luminosity of $ 0.5 \times 10^{35}\;{\rm{ cm}}^{-2}{\rm{s}}^{-1} $. A key feature of STCF is its capability to perform energy scans over the COM-energy range of 2–$ 7\;{\rm{ GeV}} $, providing extensive opportunities to search for BSM physics in the tau–charm energy region. In particular, as a dedicated charm-meson factory operating at the $ \psi(3770) $ resonance, STCF is expected to produce approximately $ 2.8 \times 10^{9} $ $ D^+D^- $ meson pairs with an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $ collected in a single year of operation. This yield is approximately two orders of magnitude larger than that of BESIII [115].

STCF comprises multiple subdetectors, including the Inner Tracking system (ITK), the Main Drift Chamber (MDC), the DIRC-like Time-of-flight detector for particle identification, and the Electromagnetic Calorimeter [93]. For charged particles in the momentum range relevant to this study (as discussed below), particle identification and track reconstruction are primarily performed by the ITK and MDC. Therefore, in the phenomenological analysis, we restrict the geometric acceptance to these two subdetectors by requiring the HNLs to decay within this region; in addition, further geometric constraints are imposed to suppress background contributions. Specifically, the geometric coverage of the ITK and MDC can be approximated as $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 6\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $, where z and r denote the longitudinal and radial decay positions of the HNL, respectively. Moreover, near the interaction point (IP), prompt tracks, soft particles, and their multiple interactions with the detector material can generate numerous detector responses, making particle identification and track reconstruction challenging in this region. To suppress these potential backgrounds, we require $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 10\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $ as our fiducial volume.

Through mixing with the electron neutrino, HNLs can decay into a variety of SM final states. In general, a higher fraction of charged particles in the HNL final states leads to more precise DV reconstruction and improved invariant-mass ($ M_{\rm{inv}} $) resolution, thereby enhancing background suppression. Therefore, we restrict our analysis to channels that contain at least two charged particles and do not involve significant missing energy: (1) leptonic $ \ ll^{\prime}\nu:\ e^{\pm}\mu^{\mp}\nu_{\mu},\ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and (2) semi-leptonic $ \ e^{\pm}M^{\mp}: e^{\pm}\pi^{\mp},\ e^{\pm}\rho^{\mp},\ e^{\pm}K^{\mp} $. We plot their branching ratios as functions of $ m_N $ in Fig. 2.

Figure 2.  (color online) Decay branching ratios of N into final states containing at least two charged particles as functions of $ m_N $.

We select the $ D^- $-meson produced in the decay of $ \psi(3770) $ as the tagging side and fully reconstruct it through SM decay channels (single tag) to suppress the $ q\bar{q} $ continuum background. Its kinematic properties provide constraints on the signal process on the opposite side. On the signal side, we search for the proposed process associated with the $ D^+ $-meson. A single ALP produced on the signal side can decay into two HNLs, which may, in principle, both decay visibly within the detector acceptance. In practice, the presence of a single DV is generally sufficient for effective background suppression, and the additional background rejection provided by the single tag further strengthens this suppression. Therefore, reconstructing only one HNL on the signal side is considered sufficient to identify the event as a signal process.

After applying the selections described above, we expect the dominant remaining background to arise from photon conversion, which is difficult to suppress using the same strategy as that employed for backgrounds originating near the IP. Despite imposing a stringent geometric requirement $ r \gt 10\;{\rm{ cm}} $, together with the DV constraint and the single-tag strategy, reconstructing only one HNL on the signal side remains susceptible to contamination from photon-conversion backgrounds. Therefore, we exclude all leptonic channels, including $ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and $ e^{\pm}\mu^{\mp}\nu_{\mu} $, because they can be directly affected by the misreconstruction of lepton pairs originating from photon-conversion processes ($ \gamma \to e^+e^- / \mu^+\mu^- $).

For the STCF COM energy considered here ($ \sqrt{s}= 3.773\;{\rm{ GeV}} $), and taking into account initial-state radiation, the produced $ D^+ $- and $ D^- $-mesons have energies approximately in the range $ 1.870 $–$ 1.885\;{\rm{ GeV}} $. The momentum of the $ \pi^+ $-meson from $ D^+ $-meson decays ($ D^+ \to \pi^+ a $) is required to be at least $ 0.1\;{\rm{ GeV}} $ to penetrate a sufficient number of layers in the MDC and produce reconstructable tracks. Therefore, we require the ALP mass to be below $ 1.6\;{\rm{ GeV}} $, which further implies that the HNL mass must be less than half of the ALP mass, specifically below $ 0.8\;{\rm{ GeV}} $. Furthermore, an ALP with a mass $ m_a $ below 1.0 GeV would make the HNLs pair-produced from ALP decays kinematically incapable of decaying into $ e^{\pm}K^{\mp} $, a scenario that we intend to avoid. Therefore, we constrain the $ m_a $ range to 1.0–1.6 GeV. Specifically, the upper bound of 1.6 GeV and a slightly lower value of 1.4 GeV are selected as benchmark points. In addition, to avoid excessively soft charged particles in the HNL decay final states, which would degrade the reconstruction efficiency, we impose a lower bound on the HNL mass of $ 0.2\;{\rm{ GeV}} $. Consequently, the range of $ m_N $ is set to 0.2–0.8 GeV.

In Fig. 2, we observe that near $ 0.8\;{\rm{ GeV}} $, the opening of the $ e^{\pm}\rho^{\mp} $ channel leads to a turning point in the total branching ratios, while the branching ratio of the $ e^{\pm}\pi^{\mp} $ channel begins to decrease more rapidly. Within the HNL mass range of interest, 0.2–$ 0.8\;{\rm{ GeV}} $, the $ e^{\pm}\pi^{\mp} $ channel clearly dominates among the semi-leptonic final states. Therefore, for semi-leptonic modes, we select only the $ e^{\pm}\pi^{\mp} $ channel as the signal final state of the HNL. Accordingly, the final-state particle combination considered in this work is $ \pi^+ e^{\pm}\pi^{\mp} $. Its physical background can arise only from $ D^{\pm} \to \pi^+\pi^- e^{\pm} \nu_e $, which can be readily rejected by the displaced-vertex requirement, as discussed above, rendering this a very clean final state. We anticipate that, in a realistic detector environment, additional backgrounds, such as combinatorial backgrounds or fake-track backgrounds, may still be present; these will be discussed in the next section.

In addition, the MDC of the STCF detector adopts a superlayer wire configuration similar to those of Belle II and BESIII. Such a structure causes the detector efficiency for reconstructing charged-particle tracks and performing particle identification (PID) to gradually decrease as the particle production point moves away from the IP, with this effect being particularly pronounced for increasing radial production positions r. In this phenomenological analysis, to account for this decrease in detector efficiency, we model it using a linear function as

This approach has been used in Refs. [67, 83, 116−119]. In the next section, we present an improved, dedicated detector-level analysis to estimate the relevant effects.

In the phenomenological analysis, we use a fast MC simulation program with phase-space decays to generate truth-level kinematic information and evaluate the signal acceptance. Specifically, the kinematic properties of the HNL are generated through a cascade decay chain and subsequent boosts into the laboratory frame. The decay probability is then determined from its proper decay length ($ c\tau $) to assess whether it decays within the fiducial volume. Using this procedure, we scan the HNL mass over the range 0.2–$ 0.8\;{\rm{ GeV}} $ and its proper decay length $ c\tau $ over the range $ 10^{-6} $–$ 10^{10}\;{\rm{ m}} $. The final expression for the expected number of signal events $ N_s $, incorporating both the geometric acceptance and the linear efficiency function $ L(r) $, is given by

where

Here, $ N_{D^+D^-} $ is the total number of $ D^+D^- $ pairs produced for an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $, and $ \sum_i {\rm{Br}}^{\,i}_{\rm tag} \cdot \epsilon^{\,i}_{\rm tag} $ denotes the sum over all SM decay channels suitable for fully reconstructing the tagging-side $ D^- $ meson, with each channel weighted by its corresponding single-tag efficiency. The factor of 2 accounts for the two D mesons in each signal event. Using the single-tag $ D^+D^- $ production events reported in Ref. [93], we estimate a total tagging efficiency of approximately 19.64%. $ {\rm{Br}}(a \to NN)= $ 100%. $ N_{{\rm{MC}}} $ in Eq. (7) is the total number of MC-simulated events. The final integral in Eq. (8) simultaneously incorporates the linear efficiency function $ L(r) $ and the decay probability within the fiducial volume. $ R_r^{N_{i,j}} = $ $ \beta_r^{N_{i,j}} \gamma c\tau =\dfrac{p_T^{N_{i,j}}}{m_N} c \tau $ is the HNL transverse decay length in the laboratory frame, where $ j=1,2 $ denotes the two HNLs produced from the ALP decay. In addition, $ z_{N_{i,j}}=r \cot \theta_{N_{i,j}} $ is the longitudinal decay position of the HNL, with r denoting the transverse decay position and θ the HNL polar angle. The function $ A(z) $ is introduced to impose the longitudinal fiducial-volume constraint, defined as

The Super Tau–Charm Facility is a next-generation $ e^+e^- $ collider proposed for construction in China. It is designed to achieve a peak luminosity of $ 0.5 \times 10^{35}\;{\rm{ cm}}^{-2}{\rm{s}}^{-1} $. A key feature of STCF is its capability to perform energy scans over the COM-energy range of 2–$ 7\;{\rm{ GeV}} $, providing extensive opportunities to search for BSM physics in the tau–charm energy region. In particular, as a dedicated charm-meson factory operating at the $ \psi(3770) $ resonance, STCF is expected to produce approximately $ 2.8 \times 10^{9} $ $ D^+D^- $ meson pairs with an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $ collected in a single year of operation. This yield is approximately two orders of magnitude larger than that of BESIII [115].

STCF comprises multiple subdetectors, including the Inner Tracking system (ITK), the Main Drift Chamber (MDC), the DIRC-like Time-of-flight detector for particle identification, and the Electromagnetic Calorimeter [93]. For charged particles in the momentum range relevant to this study (as discussed below), particle identification and track reconstruction are primarily performed by the ITK and MDC. Therefore, in the phenomenological analysis, we restrict the geometric acceptance to these two subdetectors by requiring the HNLs to decay within this region; in addition, further geometric constraints are imposed to suppress background contributions. Specifically, the geometric coverage of the ITK and MDC can be approximated as $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 6\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $, where z and r denote the longitudinal and radial decay positions of the HNL, respectively. Moreover, near the interaction point (IP), prompt tracks, soft particles, and their multiple interactions with the detector material can generate numerous detector responses, making particle identification and track reconstruction challenging in this region. To suppress these potential backgrounds, we require $ -140\; {\rm{cm}} \lt z \lt 140\; {\rm{cm}} $ and $ 10\; {\rm{cm}} \lt r \lt 85\; {\rm{cm}} $ as our fiducial volume.

Through mixing with the electron neutrino, HNLs can decay into a variety of SM final states. In general, a higher fraction of charged particles in the HNL final states leads to more precise DV reconstruction and improved invariant-mass ($ M_{\rm{inv}} $) resolution, thereby enhancing background suppression. Therefore, we restrict our analysis to channels that contain at least two charged particles and do not involve significant missing energy: (1) leptonic $ \ ll^{\prime}\nu:\ e^{\pm}\mu^{\mp}\nu_{\mu},\ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and (2) semi-leptonic $ \ e^{\pm}M^{\mp}: e^{\pm}\pi^{\mp},\ e^{\pm}\rho^{\mp},\ e^{\pm}K^{\mp} $. We plot their branching ratios as functions of $ m_N $ in Fig. 2.

Figure 2.  (color online) Decay branching ratios of N into final states containing at least two charged particles as functions of $ m_N $.

We select the $ D^- $-meson produced in the decay of $ \psi(3770) $ as the tagging side and fully reconstruct it through SM decay channels (single tag) to suppress the $ q\bar{q} $ continuum background. Its kinematic properties provide constraints on the signal process on the opposite side. On the signal side, we search for the proposed process associated with the $ D^+ $-meson. A single ALP produced on the signal side can decay into two HNLs, which may, in principle, both decay visibly within the detector acceptance. In practice, the presence of a single DV is generally sufficient for effective background suppression, and the additional background rejection provided by the single tag further strengthens this suppression. Therefore, reconstructing only one HNL on the signal side is considered sufficient to identify the event as a signal process.

After applying the selections described above, we expect the dominant remaining background to arise from photon conversion, which is difficult to suppress using the same strategy as that employed for backgrounds originating near the IP. Despite imposing a stringent geometric requirement $ r \gt 10\;{\rm{ cm}} $, together with the DV constraint and the single-tag strategy, reconstructing only one HNL on the signal side remains susceptible to contamination from photon-conversion backgrounds. Therefore, we exclude all leptonic channels, including $ e^{+}e^{-}\nu_{e},\ \mu^{+}\mu^{-}\nu_{e} $, and $ e^{\pm}\mu^{\mp}\nu_{\mu} $, because they can be directly affected by the misreconstruction of lepton pairs originating from photon-conversion processes ($ \gamma \to e^+e^- / \mu^+\mu^- $).

For the STCF COM energy considered here ($ \sqrt{s}= 3.773\;{\rm{ GeV}} $), and taking into account initial-state radiation, the produced $ D^+ $- and $ D^- $-mesons have energies approximately in the range $ 1.870 $–$ 1.885\;{\rm{ GeV}} $. The momentum of the $ \pi^+ $-meson from $ D^+ $-meson decays ($ D^+ \to \pi^+ a $) is required to be at least $ 0.1\;{\rm{ GeV}} $ to penetrate a sufficient number of layers in the MDC and produce reconstructable tracks. Therefore, we require the ALP mass to be below $ 1.6\;{\rm{ GeV}} $, which further implies that the HNL mass must be less than half of the ALP mass, specifically below $ 0.8\;{\rm{ GeV}} $. Furthermore, an ALP with a mass $ m_a $ below 1.0 GeV would make the HNLs pair-produced from ALP decays kinematically incapable of decaying into $ e^{\pm}K^{\mp} $, a scenario that we intend to avoid. Therefore, we constrain the $ m_a $ range to 1.0–1.6 GeV. Specifically, the upper bound of 1.6 GeV and a slightly lower value of 1.4 GeV are selected as benchmark points. In addition, to avoid excessively soft charged particles in the HNL decay final states, which would degrade the reconstruction efficiency, we impose a lower bound on the HNL mass of $ 0.2\;{\rm{ GeV}} $. Consequently, the range of $ m_N $ is set to 0.2–0.8 GeV.

In Fig. 2, we observe that near $ 0.8\;{\rm{ GeV}} $, the opening of the $ e^{\pm}\rho^{\mp} $ channel leads to a turning point in the total branching ratios, while the branching ratio of the $ e^{\pm}\pi^{\mp} $ channel begins to decrease more rapidly. Within the HNL mass range of interest, 0.2–$ 0.8\;{\rm{ GeV}} $, the $ e^{\pm}\pi^{\mp} $ channel clearly dominates among the semi-leptonic final states. Therefore, for semi-leptonic modes, we select only the $ e^{\pm}\pi^{\mp} $ channel as the signal final state of the HNL. Accordingly, the final-state particle combination considered in this work is $ \pi^+ e^{\pm}\pi^{\mp} $. Its physical background can arise only from $ D^{\pm} \to \pi^+\pi^- e^{\pm} \nu_e $, which can be readily rejected by the displaced-vertex requirement, as discussed above, rendering this a very clean final state. We anticipate that, in a realistic detector environment, additional backgrounds, such as combinatorial backgrounds or fake-track backgrounds, may still be present; these will be discussed in the next section.

In addition, the MDC of the STCF detector adopts a superlayer wire configuration similar to those of Belle II and BESIII. Such a structure causes the detector efficiency for reconstructing charged-particle tracks and performing particle identification (PID) to gradually decrease as the particle production point moves away from the IP, with this effect being particularly pronounced for increasing radial production positions r. In this phenomenological analysis, to account for this decrease in detector efficiency, we model it using a linear function as

This approach has been used in Refs. [67, 83, 116−119]. In the next section, we present an improved, dedicated detector-level analysis to estimate the relevant effects.

In the phenomenological analysis, we use a fast MC simulation program with phase-space decays to generate truth-level kinematic information and evaluate the signal acceptance. Specifically, the kinematic properties of the HNL are generated through a cascade decay chain and subsequent boosts into the laboratory frame. The decay probability is then determined from its proper decay length ($ c\tau $) to assess whether it decays within the fiducial volume. Using this procedure, we scan the HNL mass over the range 0.2–$ 0.8\;{\rm{ GeV}} $ and its proper decay length $ c\tau $ over the range $ 10^{-6} $–$ 10^{10}\;{\rm{ m}} $. The final expression for the expected number of signal events $ N_s $, incorporating both the geometric acceptance and the linear efficiency function $ L(r) $, is given by

where

Here, $ N_{D^+D^-} $ is the total number of $ D^+D^- $ pairs produced for an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $, and $ \sum_i {\rm{Br}}^{\,i}_{\rm tag} \cdot \epsilon^{\,i}_{\rm tag} $ denotes the sum over all SM decay channels suitable for fully reconstructing the tagging-side $ D^- $ meson, with each channel weighted by its corresponding single-tag efficiency. The factor of 2 accounts for the two D mesons in each signal event. Using the single-tag $ D^+D^- $ production events reported in Ref. [93], we estimate a total tagging efficiency of approximately 19.64%. $ {\rm{Br}}(a \to NN)= $ 100%. $ N_{{\rm{MC}}} $ in Eq. (7) is the total number of MC-simulated events. The final integral in Eq. (8) simultaneously incorporates the linear efficiency function $ L(r) $ and the decay probability within the fiducial volume. $ R_r^{N_{i,j}} = $ $ \beta_r^{N_{i,j}} \gamma c\tau =\dfrac{p_T^{N_{i,j}}}{m_N} c \tau $ is the HNL transverse decay length in the laboratory frame, where $ j=1,2 $ denotes the two HNLs produced from the ALP decay. In addition, $ z_{N_{i,j}}=r \cot \theta_{N_{i,j}} $ is the longitudinal decay position of the HNL, with r denoting the transverse decay position and θ the HNL polar angle. The function $ A(z) $ is introduced to impose the longitudinal fiducial-volume constraint, defined as

To simulate more accurately the dependence of the detector efficiency on the radial production position r presented in Sec. III, we employ the SNiPER-based OSCAR framework, which was specifically developed for STCF offline data analysis. OSCAR enables the generation of physical events, simulation of detector responses, event reconstruction, and analysis of the signal process.

During the event generation and simulation phase, we employ the KKMC generator [120, 121] to simulate $ e^+e^- $ collisions and produce the $ \psi(3770) $ resonance. We use the StcfEvtGen module to handle the production and decay processes of charm mesons in both the signal and tag channels.

In the phenomenological estimation of the total tagging efficiency discussed in the previous section, all known decay channels of the $ D^- $ mesons listed in Ref. [105] are taken into account. At the reconstructed level, we use the channel $ D^- \to K^+ \pi^- \pi^- $ as a representative benchmark for evaluating the reconstruction efficiencies of the signal channel, motivated by its relatively large branching ratio and clean charged final states. Assuming that the signal-event reconstruction efficiencies are largely independent of the specific tag channel chosen, we include all tag channels in the numerical computation of the signal-event yields and use the reconstruction efficiencies of the $ D^- \to K^+ \pi^- \pi^- $ channel for each of them.

Finally, we select candidate particles, extract relevant observables from kinematic properties, and apply selection criteria to identify signal events and compute the efficiencies.

In Table 1, we list the signal-event selection efficiency cut flows for a benchmark point with $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV, combining data samples with $ c\tau $ = 1 m, 5 m, and 10 m. On the tag side, in each event, we select a set of charged-particle candidates by using $ \Delta E_{{\rm{tag}}}= E_{K\pi\pi}-E_{{\rm{beam}}} $, the difference between the energy of the system of promptly produced charged particles and the beam energy $ E_{{\rm{beam}}}=1.8865\; {\rm{GeV}} $. Since the $ \Delta E_{{\rm{tag}}} $ distribution exhibits a sufficiently narrow peak, we impose a $ \pm 3\sigma $ requirement around its central value, where $ 1\sigma = 0.0045 $ GeV is obtained by fitting the $ \Delta E_{{\rm{tag}}} $ distribution with a Crystal Ball function. Subsequently, we compute the beam-constrained mass $ M_{{\rm{BC}}}=\sqrt{E^2_{{\rm{beam}}}-|{\boldsymbol{p}}_{D^-}|^2} $ for each candidate combination and apply a $ \pm 3\sigma $ window around the central value $ M_{{\rm{BC}}}=1.8697\; {\rm{GeV}} $ to the event yield, where $ 1\sigma = 0.0011 $ GeV is obtained in the same manner as for $ \Delta E_{{\rm{tag}}} $. Here, $ {\boldsymbol{p}}_{D^-} $ is the three-momentum of all candidate particles used to reconstruct the $ D^- $ meson in the COM system.

On the signal side, we first select charged particles originating from DVs by requiring $ d_0 \gt 0.4\; {\rm{cm}} $ and $ z_0 \gt 0.6\; {\rm{cm}} $, which are the transverse and longitudinal distances, respectively, of the point of closest approach to the IP. This requirement distinguishes displaced tracks from candidate particles on the tag side as well as from other prompt particles. We then require at least one displaced π candidate and at least one displaced e candidate, with the π and e carrying opposite electric charges. We also require at least one promptly produced π with a charge opposite to that of the tag-side $ D^- $ meson. Subsequently, a $ \Delta E_{{\rm{sig}}} $ variable, defined as the difference between the energy of the displaced charged particles and the beam energy, is used to further select these charged particles and obtain a set of signal-side candidates; however, no window requirement around the central value of $ \Delta E_{{\rm{sig}}} $ is imposed. Finally, we require the transverse momenta of both the displaced π and displaced e candidates to satisfy $ p_T \gt 0.1\; {\rm{GeV}} $ to improve the $ \pi/e $ separation capability and track-reconstruction precision. The final signal is observed through the distribution of the invariant mass $ M^2_{{\rm{inv}}}(\pi e) $ reconstructed from the $ \pi e $ pair, as illustrated in Fig. 3 for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV. We choose not to apply additional selections on $ M^2_{{\rm{inv}}}(\pi e) $, since doing so could potentially lead to a loss of signal peaks in a real experimental analysis.

Figure 3.  (color online) Distributions of the squared invariant mass of the $ \pi e $ system for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV at $ \sqrt{s}=3.773 $ GeV, obtained from OSCAR simulations. The black, blue, and red curves correspond to the raw MC samples, single-tagged samples, and full-selection samples, respectively. Samples with $ c\tau=1 $ m, 5 m, and 10 m are combined.

As shown in Table 1, after applying the PID selection criteria on the signal side, the corresponding selection efficiency decreases significantly. This reduction is mainly driven by the requirement that displaced pions and electrons be identified in front of the MDC. Specifically, this effect arises from both the decrease in detector reconstruction efficiency with increasing HNL radial decay position r and the fact that some long-lived HNLs in our event samples do not decay inside the fiducial volume. Meanwhile, in Fig. 3, the full-selection sample does not exhibit as significant a variation relative to the single-tag selected sample as that shown in Table 1. This is because, for events passing the single-tag selection but failing the signal-side PID selection, the kinematic variable $ M^2_{{\rm{inv}}}(\pi e) $ is often not properly assigned in the analysis program and therefore falls outside the displayed range of the horizontal axis.

In Fig. 3, we combine all event samples with $ c\tau=1 $ m, 5 m, and 10 m. The black, blue, and red curves correspond to the raw MC sample, the single-tag selected sample (including the $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $ selections in Table 1), and the final full-selection sample, respectively. A peak is observed at $ M^2_{{\rm{inv}}}(\pi e)=0.25 $ GeV$ ^2 $, which correctly corresponds to the HNL mass of 0.5 GeV. In the low invariant-mass region of the raw MC sample, around $ 0.1\; {\rm{GeV}}^{2} $, the invariant-mass distribution exhibits an anomalous enhancement. This is mainly because some low-energy hits are reconstructed as fake tracks, which are then incorrectly identified as e or π particles. Importantly, the distributions show that, with the selection cuts imposed on the event samples, the signal peak around 0.25 GeV² becomes more pronounced and exhibits a more Gaussian-like shape, with the tails on both sides decreasing. Therefore, the application of our selection criteria effectively suppresses such misreconstruction.

Similar to the approach employed in our phenomenological analysis, we perform a scan over the HNL mass in the range 0.2–$ 0.8\;{\rm{ GeV}} $ with a step size of $ 0.1\; {\rm{GeV}} $. However, we consider only three representative values of $ c\tau $ approximately matched to the detector scale, i.e. $ 1\; {\rm{m}} $, $ 5\; {\rm{m}} $, and $ 10\; {\rm{m}} $, to ensure sufficient statistics of HNL decays occurring within the fiducial volume. This choice is motivated by the fact that the reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ of the signal channel is intrinsically independent of $ c\tau $ and instead depends on the ALP mass, the HNL mass, and its radial decay position r within the tracker. Thus, instead of parameterizing the reconstruction efficiencies in terms of $ |V_{eN}|^2 $ or $ c\tau $, we use efficiencies as functions of r, which provide a more fundamental characterization and enable a reliable estimate without simulating prohibitively large event samples.

The reconstruction efficiencies $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ as functions of the transverse decay position $ r_{{\rm{decay}}} $ relative to the IP are shown in Fig. 4 for representative values of $ m_a $ and $ m_N $. To obtain Fig. 4, we divide the fiducial volume into radial bins of width $ 1.5 $ cm, covering the range $ 0\;{\rm{ cm}} \lt r \lt 85 $ cm, i.e., the region between the IP and the edge of the MDC. For each bin, we compute the ratio of the number of reconstructed signal events to the number of MC truth-level signal events. This ratio is defined as the signal-event reconstruction efficiency at the bin center and is shown as a point in Fig. 4.

Figure 4.  (color online) Reconstruction efficiencies $ \epsilon_{{\rm{rec}}} $ as functions of the transverse decay position from the IP, $ r_{{\rm{decay}}} $, for $ m_a=1.4 $ GeV (left) and $ m_a=1.6 $ GeV (right). For each benchmark ALP mass, the efficiencies are shown for representative HNL masses.

In this numerical simulation, we use benchmark $ c\tau $ values of $ {\cal{O}}(1) $ m, which are comparable to the detector length scale, to obtain large statistics for HNL decays within the fiducial volume. However, with this choice, the probability of both HNLs decaying inside the fiducial volume is non-negligible. Thus, when evaluating the reconstruction efficiency for one of the two HNLs produced in ALP decays, we require the other HNL to decay outside the entire fiducial volume to avoid contamination from decays occurring outside the tested bin region. This requirement is consistent with Eq. (7), which is used in our phenomenological analysis. Note that, in computing this efficiency, we select all HNL decay events within the unit volumes that pass the tag-side selection criteria, including those on $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $.

We observe that the reconstruction efficiencies are relatively low at small $ r_{{\rm{decay}}} $, reach a maximum at intermediate radii, and then gradually decrease at larger $ r_{{\rm{decay}}} $. The suppression at small $ r_{\rm decay} $ originates from the requirement that both the electron and the pion be identified as displaced tracks. For decays occurring close to the IP, the electron and pion can easily fail the requirements of $ d_0>0.4\; {\rm{cm}} $ and $ z_0>0.6\; {\rm{cm}} $ for both tracks. Moreover, lighter HNLs are typically more boosted, and their decay products tend to be more collimated along the HNL flight direction, further decreasing the probability of satisfying these requirements. At large distances from the IP, the reconstruction efficiencies decrease because of the degraded tracking capability of the tracker.

The detector-level reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a, m_N) $ is more realistic and reliable than the linear assumption $ L(r) $ in Eq. (5). Replacing $ L(r) $ in Eq. (8) with $ \epsilon_{{\rm{rec}}}(r, m_a,m_N) $ allows the detector-level signal acceptances to be extrapolated to any $ c\tau $ value. This approach has been adopted previously, for example, in a Belle search [122].

The results presented above are obtained from the full detector OSCAR simulation of the signal MC samples only. We have also analyzed inclusive MC samples of $ D^+D^- $ pairs, including all SM decay channels and corresponding to approximately $ {\cal{O}}(10^7) $ events, to obtain a qualitative understanding of the background. The results are consistent with the expectations in Sec. III: by applying the same selection criteria as those used in the signal MC simulation, the physical backgrounds can be effectively suppressed, and the remaining background is dominated by fake tracks arising from misreconstruction of detector responses, as well as by combinatorial backgrounds. These backgrounds can be further reduced by performing vertex fitting on the final-state tracks and by employing machine-learning techniques [123, 124], while maintaining the original reconstruction efficiency. It should be noted that, owing to computational resource limitations, only a limited inclusive MC sample is analyzed, which is insufficient to draw definitive conclusions about the background level corresponding to an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $. A more comprehensive background study and further optimization of the selection criteria are beyond the scope of the present work and are left for future experimental studies.

Based on the discussion in the previous and present sections, we assume vanishing background levels and derive the contours corresponding to 3 signal events, representing the sensitivity reach at the 95% confidence level (CL).

To simulate more accurately the dependence of the detector efficiency on the radial production position r presented in Sec. III, we employ the SNiPER-based OSCAR framework, which was specifically developed for STCF offline data analysis. OSCAR enables the generation of physical events, simulation of detector responses, event reconstruction, and analysis of the signal process.

During the event generation and simulation phase, we employ the KKMC generator [120, 121] to simulate $ e^+e^- $ collisions and produce the $ \psi(3770) $ resonance. We use the StcfEvtGen module to handle the production and decay processes of charm mesons in both the signal and tag channels.

In the phenomenological estimation of the total tagging efficiency discussed in the previous section, all known decay channels of the $ D^- $ mesons listed in Ref. [105] are taken into account. At the reconstructed level, we use the channel $ D^- \to K^+ \pi^- \pi^- $ as a representative benchmark for evaluating the reconstruction efficiencies of the signal channel, motivated by its relatively large branching ratio and clean charged final states. Assuming that the signal-event reconstruction efficiencies are largely independent of the specific tag channel chosen, we include all tag channels in the numerical computation of the signal-event yields and use the reconstruction efficiencies of the $ D^- \to K^+ \pi^- \pi^- $ channel for each of them.

Finally, we select candidate particles, extract relevant observables from kinematic properties, and apply selection criteria to identify signal events and compute the efficiencies.

In Table 1, we list the signal-event selection efficiency cut flows for a benchmark point with $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV, combining data samples with $ c\tau $ = 1 m, 5 m, and 10 m. On the tag side, in each event, we select a set of charged-particle candidates by using $ \Delta E_{{\rm{tag}}}= E_{K\pi\pi}-E_{{\rm{beam}}} $, the difference between the energy of the system of promptly produced charged particles and the beam energy $ E_{{\rm{beam}}}=1.8865\; {\rm{GeV}} $. Since the $ \Delta E_{{\rm{tag}}} $ distribution exhibits a sufficiently narrow peak, we impose a $ \pm 3\sigma $ requirement around its central value, where $ 1\sigma = 0.0045 $ GeV is obtained by fitting the $ \Delta E_{{\rm{tag}}} $ distribution with a Crystal Ball function. Subsequently, we compute the beam-constrained mass $ M_{{\rm{BC}}}=\sqrt{E^2_{{\rm{beam}}}-|{\boldsymbol{p}}_{D^-}|^2} $ for each candidate combination and apply a $ \pm 3\sigma $ window around the central value $ M_{{\rm{BC}}}=1.8697\; {\rm{GeV}} $ to the event yield, where $ 1\sigma = 0.0011 $ GeV is obtained in the same manner as for $ \Delta E_{{\rm{tag}}} $. Here, $ {\boldsymbol{p}}_{D^-} $ is the three-momentum of all candidate particles used to reconstruct the $ D^- $ meson in the COM system.

On the signal side, we first select charged particles originating from DVs by requiring $ d_0 \gt 0.4\; {\rm{cm}} $ and $ z_0 \gt 0.6\; {\rm{cm}} $, which are the transverse and longitudinal distances, respectively, of the point of closest approach to the IP. This requirement distinguishes displaced tracks from candidate particles on the tag side as well as from other prompt particles. We then require at least one displaced π candidate and at least one displaced e candidate, with the π and e carrying opposite electric charges. We also require at least one promptly produced π with a charge opposite to that of the tag-side $ D^- $ meson. Subsequently, a $ \Delta E_{{\rm{sig}}} $ variable, defined as the difference between the energy of the displaced charged particles and the beam energy, is used to further select these charged particles and obtain a set of signal-side candidates; however, no window requirement around the central value of $ \Delta E_{{\rm{sig}}} $ is imposed. Finally, we require the transverse momenta of both the displaced π and displaced e candidates to satisfy $ p_T \gt 0.1\; {\rm{GeV}} $ to improve the $ \pi/e $ separation capability and track-reconstruction precision. The final signal is observed through the distribution of the invariant mass $ M^2_{{\rm{inv}}}(\pi e) $ reconstructed from the $ \pi e $ pair, as illustrated in Fig. 3 for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV. We choose not to apply additional selections on $ M^2_{{\rm{inv}}}(\pi e) $, since doing so could potentially lead to a loss of signal peaks in a real experimental analysis.

Figure 3.  (color online) Distributions of the squared invariant mass of the $ \pi e $ system for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV at $ \sqrt{s}=3.773 $ GeV, obtained from OSCAR simulations. The black, blue, and red curves correspond to the raw MC samples, single-tagged samples, and full-selection samples, respectively. Samples with $ c\tau=1 $ m, 5 m, and 10 m are combined.

As shown in Table 1, after applying the PID selection criteria on the signal side, the corresponding selection efficiency decreases significantly. This reduction is mainly driven by the requirement that displaced pions and electrons be identified in front of the MDC. Specifically, this effect arises from both the decrease in detector reconstruction efficiency with increasing HNL radial decay position r and the fact that some long-lived HNLs in our event samples do not decay inside the fiducial volume. Meanwhile, in Fig. 3, the full-selection sample does not exhibit as significant a variation relative to the single-tag selected sample as that shown in Table 1. This is because, for events passing the single-tag selection but failing the signal-side PID selection, the kinematic variable $ M^2_{{\rm{inv}}}(\pi e) $ is often not properly assigned in the analysis program and therefore falls outside the displayed range of the horizontal axis.

In Fig. 3, we combine all event samples with $ c\tau=1 $ m, 5 m, and 10 m. The black, blue, and red curves correspond to the raw MC sample, the single-tag selected sample (including the $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $ selections in Table 1), and the final full-selection sample, respectively. A peak is observed at $ M^2_{{\rm{inv}}}(\pi e)=0.25 $ GeV$ ^2 $, which correctly corresponds to the HNL mass of 0.5 GeV. In the low invariant-mass region of the raw MC sample, around $ 0.1\; {\rm{GeV}}^{2} $, the invariant-mass distribution exhibits an anomalous enhancement. This is mainly because some low-energy hits are reconstructed as fake tracks, which are then incorrectly identified as e or π particles. Importantly, the distributions show that, with the selection cuts imposed on the event samples, the signal peak around 0.25 GeV² becomes more pronounced and exhibits a more Gaussian-like shape, with the tails on both sides decreasing. Therefore, the application of our selection criteria effectively suppresses such misreconstruction.

Similar to the approach employed in our phenomenological analysis, we perform a scan over the HNL mass in the range 0.2–$ 0.8\;{\rm{ GeV}} $ with a step size of $ 0.1\; {\rm{GeV}} $. However, we consider only three representative values of $ c\tau $ approximately matched to the detector scale, i.e. $ 1\; {\rm{m}} $, $ 5\; {\rm{m}} $, and $ 10\; {\rm{m}} $, to ensure sufficient statistics of HNL decays occurring within the fiducial volume. This choice is motivated by the fact that the reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ of the signal channel is intrinsically independent of $ c\tau $ and instead depends on the ALP mass, the HNL mass, and its radial decay position r within the tracker. Thus, instead of parameterizing the reconstruction efficiencies in terms of $ |V_{eN}|^2 $ or $ c\tau $, we use efficiencies as functions of r, which provide a more fundamental characterization and enable a reliable estimate without simulating prohibitively large event samples.

The reconstruction efficiencies $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ as functions of the transverse decay position $ r_{{\rm{decay}}} $ relative to the IP are shown in Fig. 4 for representative values of $ m_a $ and $ m_N $. To obtain Fig. 4, we divide the fiducial volume into radial bins of width $ 1.5 $ cm, covering the range $ 0\;{\rm{ cm}} \lt r \lt 85 $ cm, i.e., the region between the IP and the edge of the MDC. For each bin, we compute the ratio of the number of reconstructed signal events to the number of MC truth-level signal events. This ratio is defined as the signal-event reconstruction efficiency at the bin center and is shown as a point in Fig. 4.

Figure 4.  (color online) Reconstruction efficiencies $ \epsilon_{{\rm{rec}}} $ as functions of the transverse decay position from the IP, $ r_{{\rm{decay}}} $, for $ m_a=1.4 $ GeV (left) and $ m_a=1.6 $ GeV (right). For each benchmark ALP mass, the efficiencies are shown for representative HNL masses.

In this numerical simulation, we use benchmark $ c\tau $ values of $ {\cal{O}}(1) $ m, which are comparable to the detector length scale, to obtain large statistics for HNL decays within the fiducial volume. However, with this choice, the probability of both HNLs decaying inside the fiducial volume is non-negligible. Thus, when evaluating the reconstruction efficiency for one of the two HNLs produced in ALP decays, we require the other HNL to decay outside the entire fiducial volume to avoid contamination from decays occurring outside the tested bin region. This requirement is consistent with Eq. (7), which is used in our phenomenological analysis. Note that, in computing this efficiency, we select all HNL decay events within the unit volumes that pass the tag-side selection criteria, including those on $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $.

We observe that the reconstruction efficiencies are relatively low at small $ r_{{\rm{decay}}} $, reach a maximum at intermediate radii, and then gradually decrease at larger $ r_{{\rm{decay}}} $. The suppression at small $ r_{\rm decay} $ originates from the requirement that both the electron and the pion be identified as displaced tracks. For decays occurring close to the IP, the electron and pion can easily fail the requirements of $ d_0>0.4\; {\rm{cm}} $ and $ z_0>0.6\; {\rm{cm}} $ for both tracks. Moreover, lighter HNLs are typically more boosted, and their decay products tend to be more collimated along the HNL flight direction, further decreasing the probability of satisfying these requirements. At large distances from the IP, the reconstruction efficiencies decrease because of the degraded tracking capability of the tracker.

The detector-level reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a, m_N) $ is more realistic and reliable than the linear assumption $ L(r) $ in Eq. (5). Replacing $ L(r) $ in Eq. (8) with $ \epsilon_{{\rm{rec}}}(r, m_a,m_N) $ allows the detector-level signal acceptances to be extrapolated to any $ c\tau $ value. This approach has been adopted previously, for example, in a Belle search [122].

The results presented above are obtained from the full detector OSCAR simulation of the signal MC samples only. We have also analyzed inclusive MC samples of $ D^+D^- $ pairs, including all SM decay channels and corresponding to approximately $ {\cal{O}}(10^7) $ events, to obtain a qualitative understanding of the background. The results are consistent with the expectations in Sec. III: by applying the same selection criteria as those used in the signal MC simulation, the physical backgrounds can be effectively suppressed, and the remaining background is dominated by fake tracks arising from misreconstruction of detector responses, as well as by combinatorial backgrounds. These backgrounds can be further reduced by performing vertex fitting on the final-state tracks and by employing machine-learning techniques [123, 124], while maintaining the original reconstruction efficiency. It should be noted that, owing to computational resource limitations, only a limited inclusive MC sample is analyzed, which is insufficient to draw definitive conclusions about the background level corresponding to an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $. A more comprehensive background study and further optimization of the selection criteria are beyond the scope of the present work and are left for future experimental studies.

Based on the discussion in the previous and present sections, we assume vanishing background levels and derive the contours corresponding to 3 signal events, representing the sensitivity reach at the 95% confidence level (CL).

To simulate more accurately the dependence of the detector efficiency on the radial production position r presented in Sec. III, we employ the SNiPER-based OSCAR framework, which was specifically developed for STCF offline data analysis. OSCAR enables the generation of physical events, simulation of detector responses, event reconstruction, and analysis of the signal process.

During the event generation and simulation phase, we employ the KKMC generator [120, 121] to simulate $ e^+e^- $ collisions and produce the $ \psi(3770) $ resonance. We use the StcfEvtGen module to handle the production and decay processes of charm mesons in both the signal and tag channels.

In the phenomenological estimation of the total tagging efficiency discussed in the previous section, all known decay channels of the $ D^- $ mesons listed in Ref. [105] are taken into account. At the reconstructed level, we use the channel $ D^- \to K^+ \pi^- \pi^- $ as a representative benchmark for evaluating the reconstruction efficiencies of the signal channel, motivated by its relatively large branching ratio and clean charged final states. Assuming that the signal-event reconstruction efficiencies are largely independent of the specific tag channel chosen, we include all tag channels in the numerical computation of the signal-event yields and use the reconstruction efficiencies of the $ D^- \to K^+ \pi^- \pi^- $ channel for each of them.

Finally, we select candidate particles, extract relevant observables from kinematic properties, and apply selection criteria to identify signal events and compute the efficiencies.

In Table 1, we list the signal-event selection efficiency cut flows for a benchmark point with $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV, combining data samples with $ c\tau $ = 1 m, 5 m, and 10 m. On the tag side, in each event, we select a set of charged-particle candidates by using $ \Delta E_{{\rm{tag}}}= E_{K\pi\pi}-E_{{\rm{beam}}} $, the difference between the energy of the system of promptly produced charged particles and the beam energy $ E_{{\rm{beam}}}=1.8865\; {\rm{GeV}} $. Since the $ \Delta E_{{\rm{tag}}} $ distribution exhibits a sufficiently narrow peak, we impose a $ \pm 3\sigma $ requirement around its central value, where $ 1\sigma = 0.0045 $ GeV is obtained by fitting the $ \Delta E_{{\rm{tag}}} $ distribution with a Crystal Ball function. Subsequently, we compute the beam-constrained mass $ M_{{\rm{BC}}}=\sqrt{E^2_{{\rm{beam}}}-|{\boldsymbol{p}}_{D^-}|^2} $ for each candidate combination and apply a $ \pm 3\sigma $ window around the central value $ M_{{\rm{BC}}}=1.8697\; {\rm{GeV}} $ to the event yield, where $ 1\sigma = 0.0011 $ GeV is obtained in the same manner as for $ \Delta E_{{\rm{tag}}} $. Here, $ {\boldsymbol{p}}_{D^-} $ is the three-momentum of all candidate particles used to reconstruct the $ D^- $ meson in the COM system.

On the signal side, we first select charged particles originating from DVs by requiring $ d_0 \gt 0.4\; {\rm{cm}} $ and $ z_0 \gt 0.6\; {\rm{cm}} $, which are the transverse and longitudinal distances, respectively, of the point of closest approach to the IP. This requirement distinguishes displaced tracks from candidate particles on the tag side as well as from other prompt particles. We then require at least one displaced π candidate and at least one displaced e candidate, with the π and e carrying opposite electric charges. We also require at least one promptly produced π with a charge opposite to that of the tag-side $ D^- $ meson. Subsequently, a $ \Delta E_{{\rm{sig}}} $ variable, defined as the difference between the energy of the displaced charged particles and the beam energy, is used to further select these charged particles and obtain a set of signal-side candidates; however, no window requirement around the central value of $ \Delta E_{{\rm{sig}}} $ is imposed. Finally, we require the transverse momenta of both the displaced π and displaced e candidates to satisfy $ p_T \gt 0.1\; {\rm{GeV}} $ to improve the $ \pi/e $ separation capability and track-reconstruction precision. The final signal is observed through the distribution of the invariant mass $ M^2_{{\rm{inv}}}(\pi e) $ reconstructed from the $ \pi e $ pair, as illustrated in Fig. 3 for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV. We choose not to apply additional selections on $ M^2_{{\rm{inv}}}(\pi e) $, since doing so could potentially lead to a loss of signal peaks in a real experimental analysis.

Figure 3.  (color online) Distributions of the squared invariant mass of the $ \pi e $ system for $ m_a=1.4 $ GeV and $ m_N=0.5 $ GeV at $ \sqrt{s}=3.773 $ GeV, obtained from OSCAR simulations. The black, blue, and red curves correspond to the raw MC samples, single-tagged samples, and full-selection samples, respectively. Samples with $ c\tau=1 $ m, 5 m, and 10 m are combined.

As shown in Table 1, after applying the PID selection criteria on the signal side, the corresponding selection efficiency decreases significantly. This reduction is mainly driven by the requirement that displaced pions and electrons be identified in front of the MDC. Specifically, this effect arises from both the decrease in detector reconstruction efficiency with increasing HNL radial decay position r and the fact that some long-lived HNLs in our event samples do not decay inside the fiducial volume. Meanwhile, in Fig. 3, the full-selection sample does not exhibit as significant a variation relative to the single-tag selected sample as that shown in Table 1. This is because, for events passing the single-tag selection but failing the signal-side PID selection, the kinematic variable $ M^2_{{\rm{inv}}}(\pi e) $ is often not properly assigned in the analysis program and therefore falls outside the displayed range of the horizontal axis.

In Fig. 3, we combine all event samples with $ c\tau=1 $ m, 5 m, and 10 m. The black, blue, and red curves correspond to the raw MC sample, the single-tag selected sample (including the $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $ selections in Table 1), and the final full-selection sample, respectively. A peak is observed at $ M^2_{{\rm{inv}}}(\pi e)=0.25 $ GeV$ ^2 $, which correctly corresponds to the HNL mass of 0.5 GeV. In the low invariant-mass region of the raw MC sample, around $ 0.1\; {\rm{GeV}}^{2} $, the invariant-mass distribution exhibits an anomalous enhancement. This is mainly because some low-energy hits are reconstructed as fake tracks, which are then incorrectly identified as e or π particles. Importantly, the distributions show that, with the selection cuts imposed on the event samples, the signal peak around 0.25 GeV² becomes more pronounced and exhibits a more Gaussian-like shape, with the tails on both sides decreasing. Therefore, the application of our selection criteria effectively suppresses such misreconstruction.

Similar to the approach employed in our phenomenological analysis, we perform a scan over the HNL mass in the range 0.2–$ 0.8\;{\rm{ GeV}} $ with a step size of $ 0.1\; {\rm{GeV}} $. However, we consider only three representative values of $ c\tau $ approximately matched to the detector scale, i.e. $ 1\; {\rm{m}} $, $ 5\; {\rm{m}} $, and $ 10\; {\rm{m}} $, to ensure sufficient statistics of HNL decays occurring within the fiducial volume. This choice is motivated by the fact that the reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ of the signal channel is intrinsically independent of $ c\tau $ and instead depends on the ALP mass, the HNL mass, and its radial decay position r within the tracker. Thus, instead of parameterizing the reconstruction efficiencies in terms of $ |V_{eN}|^2 $ or $ c\tau $, we use efficiencies as functions of r, which provide a more fundamental characterization and enable a reliable estimate without simulating prohibitively large event samples.

The reconstruction efficiencies $ \epsilon_{{\rm{rec}}}(r,m_a,m_N) $ as functions of the transverse decay position $ r_{{\rm{decay}}} $ relative to the IP are shown in Fig. 4 for representative values of $ m_a $ and $ m_N $. To obtain Fig. 4, we divide the fiducial volume into radial bins of width $ 1.5 $ cm, covering the range $ 0\;{\rm{ cm}} \lt r \lt 85 $ cm, i.e., the region between the IP and the edge of the MDC. For each bin, we compute the ratio of the number of reconstructed signal events to the number of MC truth-level signal events. This ratio is defined as the signal-event reconstruction efficiency at the bin center and is shown as a point in Fig. 4.

Figure 4.  (color online) Reconstruction efficiencies $ \epsilon_{{\rm{rec}}} $ as functions of the transverse decay position from the IP, $ r_{{\rm{decay}}} $, for $ m_a=1.4 $ GeV (left) and $ m_a=1.6 $ GeV (right). For each benchmark ALP mass, the efficiencies are shown for representative HNL masses.

In this numerical simulation, we use benchmark $ c\tau $ values of $ {\cal{O}}(1) $ m, which are comparable to the detector length scale, to obtain large statistics for HNL decays within the fiducial volume. However, with this choice, the probability of both HNLs decaying inside the fiducial volume is non-negligible. Thus, when evaluating the reconstruction efficiency for one of the two HNLs produced in ALP decays, we require the other HNL to decay outside the entire fiducial volume to avoid contamination from decays occurring outside the tested bin region. This requirement is consistent with Eq. (7), which is used in our phenomenological analysis. Note that, in computing this efficiency, we select all HNL decay events within the unit volumes that pass the tag-side selection criteria, including those on $ \Delta E_{\rm{tag}} $ and $ M_{{\rm{BC}}} $.

We observe that the reconstruction efficiencies are relatively low at small $ r_{{\rm{decay}}} $, reach a maximum at intermediate radii, and then gradually decrease at larger $ r_{{\rm{decay}}} $. The suppression at small $ r_{\rm decay} $ originates from the requirement that both the electron and the pion be identified as displaced tracks. For decays occurring close to the IP, the electron and pion can easily fail the requirements of $ d_0>0.4\; {\rm{cm}} $ and $ z_0>0.6\; {\rm{cm}} $ for both tracks. Moreover, lighter HNLs are typically more boosted, and their decay products tend to be more collimated along the HNL flight direction, further decreasing the probability of satisfying these requirements. At large distances from the IP, the reconstruction efficiencies decrease because of the degraded tracking capability of the tracker.

The detector-level reconstruction efficiency $ \epsilon_{{\rm{rec}}}(r,m_a, m_N) $ is more realistic and reliable than the linear assumption $ L(r) $ in Eq. (5). Replacing $ L(r) $ in Eq. (8) with $ \epsilon_{{\rm{rec}}}(r, m_a,m_N) $ allows the detector-level signal acceptances to be extrapolated to any $ c\tau $ value. This approach has been adopted previously, for example, in a Belle search [122].

The results presented above are obtained from the full detector OSCAR simulation of the signal MC samples only. We have also analyzed inclusive MC samples of $ D^+D^- $ pairs, including all SM decay channels and corresponding to approximately $ {\cal{O}}(10^7) $ events, to obtain a qualitative understanding of the background. The results are consistent with the expectations in Sec. III: by applying the same selection criteria as those used in the signal MC simulation, the physical backgrounds can be effectively suppressed, and the remaining background is dominated by fake tracks arising from misreconstruction of detector responses, as well as by combinatorial backgrounds. These backgrounds can be further reduced by performing vertex fitting on the final-state tracks and by employing machine-learning techniques [123, 124], while maintaining the original reconstruction efficiency. It should be noted that, owing to computational resource limitations, only a limited inclusive MC sample is analyzed, which is insufficient to draw definitive conclusions about the background level corresponding to an integrated luminosity of $ 1\; {\rm{ab}}^{-1} $. A more comprehensive background study and further optimization of the selection criteria are beyond the scope of the present work and are left for future experimental studies.

Based on the discussion in the previous and present sections, we assume vanishing background levels and derive the contours corresponding to 3 signal events, representing the sensitivity reach at the 95% confidence level (CL).

We present the sensitivity reach of STCF with 1 ab⁻¹ to long-lived HNLs as a function of $ m_N $ in Fig. 5. We fix the couplings $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹, ensuring prompt ALP decays. For the two benchmark values of $ m_a $ considered, 1.4 GeV (dashed curves) and 1.6 GeV (dot-dashed curves), we show the numerical results obtained from the phenomenological analysis (blue) and the realistic OSCAR-simulation analysis (red). The realistic sensitivities are weaker than those from the simplified phenomenological analysis. This is mainly because the position-dependent efficiencies estimated in the OSCAR-based analysis are somewhat lower than those associated with the linear function, Eq. (5), used in the simplified analysis.

Figure 5.  (color online) Sensitivity reach of the STCF with an integrated luminosity of 1 ab⁻¹ for long-lived HNLs from prompt ALP decays produced in $ D^+ $-meson decays. We fix $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹. The results shown correspond to the benchmark values $ m_a=1.4 $ GeV (dashed) and 1.6 GeV (dot-dashed). The red and blue curves are obtained with the full OSCAR simulation and a fast truth-level simulation, respectively. The existing bounds from PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown in dark gray, and the approximate CHARM search reinterpretation results, obtained under simplified assumptions about the production and kinematic distributions, are displayed as gray dot-dashed lines.

The existing bounds obtained at PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown as the dark-gray region, while the approximate reinterpretation of the CHARM search [110], obtained under simplified assumptions about the production mechanism and kinematic distributions of the HNLs, is shown by the dot-dashed gray lines. Furthermore, we plot a light-gray band corresponding to the parameter region targeted by the type-I seesaw mechanism, for active-neutrino masses between 0.05 eV and 0.12 eV. The lower and upper bounds on the active-neutrino masses originate from neutrino-oscillation experiments [125] and cosmological observations [126], respectively.

These results show that STCF is expected to probe the mixing parameter $ |V_{eN}|^2 $ about one to two orders of magnitude beyond the existing bounds for $ m_N $ between 0.5 GeV and 0.8 GeV, although it does not yet reach the theoretical parameter space predicted by the type-I seesaw mechanism (the light-gray band in Fig. 5). We note, however, that the approximate reinterpretation of the CHARM search [110] in the present ALP–HNL framework is subject to model dependence, as it relies on assumptions about the production mechanism and kinematic distributions of the HNLs. In particular, the presence of an intermediate ALP modifies both the production and decay kinematics of the HNLs, which can affect the acceptance of beam-dump experiments. Therefore, the reinterpretation shown here should be regarded as indicative rather than a definitive exclusion. A dedicated recast including a full detector simulation would be required for a robust comparison, which is beyond the scope of the present work. Refs. [66, 67] investigated the expected signals of the ALP-HNL model at other future experiments, such as LHC far detectors and SHiP, as well as at the ongoing Belle II experiment, providing exclusion limits on $ |V_{eN}|^2 $. Ref. [67] studied a benchmark scenario with an ALP produced in B-meson decays rather than D-meson decays, while Ref. [66] investigated benchmark scenarios including ours. However, even in the latter case, the ALP mass values considered differ from ours. Therefore, our study should be regarded as highly complementary in terms of parameter-space coverage. For these reasons, we do not explicitly overlay such existing results in our sensitivity plot. Nevertheless, it can be inferred from the results in Ref. [66] that, for our benchmark scenario, future LHC far-detector and SHiP experiments are expected to probe values of $ |V_{eN}|^2 $ orders of magnitude below those to which STCF is sensitive. In this context, STCF provides a complementary probe of this scenario, with a clean experimental environment and well-controlled kinematics at threshold.

For smaller values of $ g^u_{2,1}/\Lambda $, the sensitivities are expected to weaken. Quantitatively, if $ g^u_{2,1}/\Lambda $ is reduced by a factor of 10, the sensitivity reach for $ |V_{eN}|^2 $ is expected to weaken by a factor of 100.

We present the sensitivity reach of STCF with 1 ab⁻¹ to long-lived HNLs as a function of $ m_N $ in Fig. 5. We fix the couplings $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹, ensuring prompt ALP decays. For the two benchmark values of $ m_a $ considered, 1.4 GeV (dashed curves) and 1.6 GeV (dot-dashed curves), we show the numerical results obtained from the phenomenological analysis (blue) and the realistic OSCAR-simulation analysis (red). The realistic sensitivities are weaker than those from the simplified phenomenological analysis. This is mainly because the position-dependent efficiencies estimated in the OSCAR-based analysis are somewhat lower than those associated with the linear function, Eq. (5), used in the simplified analysis.

Figure 5.  (color online) Sensitivity reach of the STCF with an integrated luminosity of 1 ab⁻¹ for long-lived HNLs from prompt ALP decays produced in $ D^+ $-meson decays. We fix $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹. The results shown correspond to the benchmark values $ m_a=1.4 $ GeV (dashed) and 1.6 GeV (dot-dashed). The red and blue curves are obtained with the full OSCAR simulation and a fast truth-level simulation, respectively. The existing bounds from PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown in dark gray, and the approximate CHARM search reinterpretation results, obtained under simplified assumptions about the production and kinematic distributions, are displayed as gray dot-dashed lines.

The existing bounds obtained at PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown as the dark-gray region, while the approximate reinterpretation of the CHARM search [110], obtained under simplified assumptions about the production mechanism and kinematic distributions of the HNLs, is shown by the dot-dashed gray lines. Furthermore, we plot a light-gray band corresponding to the parameter region targeted by the type-I seesaw mechanism, for active-neutrino masses between 0.05 eV and 0.12 eV. The lower and upper bounds on the active-neutrino masses originate from neutrino-oscillation experiments [125] and cosmological observations [126], respectively.

These results show that STCF is expected to probe the mixing parameter $ |V_{eN}|^2 $ about one to two orders of magnitude beyond the existing bounds for $ m_N $ between 0.5 GeV and 0.8 GeV, although it does not yet reach the theoretical parameter space predicted by the type-I seesaw mechanism (the light-gray band in Fig. 5). We note, however, that the approximate reinterpretation of the CHARM search [110] in the present ALP–HNL framework is subject to model dependence, as it relies on assumptions about the production mechanism and kinematic distributions of the HNLs. In particular, the presence of an intermediate ALP modifies both the production and decay kinematics of the HNLs, which can affect the acceptance of beam-dump experiments. Therefore, the reinterpretation shown here should be regarded as indicative rather than a definitive exclusion. A dedicated recast including a full detector simulation would be required for a robust comparison, which is beyond the scope of the present work. Refs. [66, 67] investigated the expected signals of the ALP-HNL model at other future experiments, such as LHC far detectors and SHiP, as well as at the ongoing Belle II experiment, providing exclusion limits on $ |V_{eN}|^2 $. Ref. [67] studied a benchmark scenario with an ALP produced in B-meson decays rather than D-meson decays, while Ref. [66] investigated benchmark scenarios including ours. However, even in the latter case, the ALP mass values considered differ from ours. Therefore, our study should be regarded as highly complementary in terms of parameter-space coverage. For these reasons, we do not explicitly overlay such existing results in our sensitivity plot. Nevertheless, it can be inferred from the results in Ref. [66] that, for our benchmark scenario, future LHC far-detector and SHiP experiments are expected to probe values of $ |V_{eN}|^2 $ orders of magnitude below those to which STCF is sensitive. In this context, STCF provides a complementary probe of this scenario, with a clean experimental environment and well-controlled kinematics at threshold.

For smaller values of $ g^u_{2,1}/\Lambda $, the sensitivities are expected to weaken. Quantitatively, if $ g^u_{2,1}/\Lambda $ is reduced by a factor of 10, the sensitivity reach for $ |V_{eN}|^2 $ is expected to weaken by a factor of 100.

We present the sensitivity reach of STCF with 1 ab⁻¹ to long-lived HNLs as a function of $ m_N $ in Fig. 5. We fix the couplings $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹, ensuring prompt ALP decays. For the two benchmark values of $ m_a $ considered, 1.4 GeV (dashed curves) and 1.6 GeV (dot-dashed curves), we show the numerical results obtained from the phenomenological analysis (blue) and the realistic OSCAR-simulation analysis (red). The realistic sensitivities are weaker than those from the simplified phenomenological analysis. This is mainly because the position-dependent efficiencies estimated in the OSCAR-based analysis are somewhat lower than those associated with the linear function, Eq. (5), used in the simplified analysis.

Figure 5.  (color online) Sensitivity reach of the STCF with an integrated luminosity of 1 ab⁻¹ for long-lived HNLs from prompt ALP decays produced in $ D^+ $-meson decays. We fix $ g^u_{2,1}/\Lambda $ at $ 2\times 10^{-4} $ TeV⁻¹ and $ g_N/\Lambda $ at $ 10^{-3} $ GeV⁻¹. The results shown correspond to the benchmark values $ m_a=1.4 $ GeV (dashed) and 1.6 GeV (dot-dashed). The red and blue curves are obtained with the full OSCAR simulation and a fast truth-level simulation, respectively. The existing bounds from PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown in dark gray, and the approximate CHARM search reinterpretation results, obtained under simplified assumptions about the production and kinematic distributions, are displayed as gray dot-dashed lines.

The existing bounds obtained at PIENU [108], KENU [109], CHARM [110], NA62 [111, 112], T2K [113], and BEBC [114] are shown as the dark-gray region, while the approximate reinterpretation of the CHARM search [110], obtained under simplified assumptions about the production mechanism and kinematic distributions of the HNLs, is shown by the dot-dashed gray lines. Furthermore, we plot a light-gray band corresponding to the parameter region targeted by the type-I seesaw mechanism, for active-neutrino masses between 0.05 eV and 0.12 eV. The lower and upper bounds on the active-neutrino masses originate from neutrino-oscillation experiments [125] and cosmological observations [126], respectively.

These results show that STCF is expected to probe the mixing parameter $ |V_{eN}|^2 $ about one to two orders of magnitude beyond the existing bounds for $ m_N $ between 0.5 GeV and 0.8 GeV, although it does not yet reach the theoretical parameter space predicted by the type-I seesaw mechanism (the light-gray band in Fig. 5). We note, however, that the approximate reinterpretation of the CHARM search [110] in the present ALP–HNL framework is subject to model dependence, as it relies on assumptions about the production mechanism and kinematic distributions of the HNLs. In particular, the presence of an intermediate ALP modifies both the production and decay kinematics of the HNLs, which can affect the acceptance of beam-dump experiments. Therefore, the reinterpretation shown here should be regarded as indicative rather than a definitive exclusion. A dedicated recast including a full detector simulation would be required for a robust comparison, which is beyond the scope of the present work. Refs. [66, 67] investigated the expected signals of the ALP-HNL model at other future experiments, such as LHC far detectors and SHiP, as well as at the ongoing Belle II experiment, providing exclusion limits on $ |V_{eN}|^2 $. Ref. [67] studied a benchmark scenario with an ALP produced in B-meson decays rather than D-meson decays, while Ref. [66] investigated benchmark scenarios including ours. However, even in the latter case, the ALP mass values considered differ from ours. Therefore, our study should be regarded as highly complementary in terms of parameter-space coverage. For these reasons, we do not explicitly overlay such existing results in our sensitivity plot. Nevertheless, it can be inferred from the results in Ref. [66] that, for our benchmark scenario, future LHC far-detector and SHiP experiments are expected to probe values of $ |V_{eN}|^2 $ orders of magnitude below those to which STCF is sensitive. In this context, STCF provides a complementary probe of this scenario, with a clean experimental environment and well-controlled kinematics at threshold.

For smaller values of $ g^u_{2,1}/\Lambda $, the sensitivities are expected to weaken. Quantitatively, if $ g^u_{2,1}/\Lambda $ is reduced by a factor of 10, the sensitivity reach for $ |V_{eN}|^2 $ is expected to weaken by a factor of 100.

In this work, we considered an ALP with off-diagonal couplings to the charm and up quarks and to a pair of long-lived Majorana HNLs. We performed both fast and dedicated MC simulations to assess the search prospects for long-lived HNLs at the proposed Super Tau–Charm Facility.

We focused on the center-of-mass energy mode of $ \sqrt{s}=3.773 $ GeV at STCF and on the decays $ D^+\to \pi^+ a $ and $ a\to NN $. The ALP a is on shell and decays promptly and exclusively into a pair of N's. Using both a simple phenomenological analysis and a dedicated reconstruction-level analysis, we determined the sensitivity reach for the active-neutrino–HNL mixing parameter $ |V_{eN}|^2 $ as a function of the HNL mass, with the branching ratio of $ D^+\to \pi^+ a $ fixed at the experimentally allowed upper bound by assuming a sufficiently small value of $ g^u_{2,1}/\Lambda $.

Our numerical results show that STCF can probe values of $ |V_{eN}|^2 $ down to approximately one to two orders of magnitude below the existing bounds, while an approximate reinterpretation of a CHARM search is found to surpass the STCF sensitivities. A complete and dedicated recasting and reinterpretation procedure would be required to establish such exclusion bounds; however, this is beyond the scope of the present work. Moreover, for benchmark scenarios in which the HNL mixes purely with $ \nu_\mu $ ($ \nu_\tau $) instead of $ \nu_e $, the sensitivities are expected to be comparable (weaker), given the similar (distinct) HNL decay products.

Finally, we note that our analysis is readily generalizable to additional theoretical scenarios, such as a scenario with both a new scalar and a dark photon, and the dark axion portal scenario.

In this work, we considered an ALP with off-diagonal couplings to the charm and up quarks and to a pair of long-lived Majorana HNLs. We performed both fast and dedicated MC simulations to assess the search prospects for long-lived HNLs at the proposed Super Tau–Charm Facility.

We focused on the center-of-mass energy mode of $ \sqrt{s}=3.773 $ GeV at STCF and on the decays $ D^+\to \pi^+ a $ and $ a\to NN $. The ALP a is on shell and decays promptly and exclusively into a pair of N's. Using both a simple phenomenological analysis and a dedicated reconstruction-level analysis, we determined the sensitivity reach for the active-neutrino–HNL mixing parameter $ |V_{eN}|^2 $ as a function of the HNL mass, with the branching ratio of $ D^+\to \pi^+ a $ fixed at the experimentally allowed upper bound by assuming a sufficiently small value of $ g^u_{2,1}/\Lambda $.

Our numerical results show that STCF can probe values of $ |V_{eN}|^2 $ down to approximately one to two orders of magnitude below the existing bounds, while an approximate reinterpretation of a CHARM search is found to surpass the STCF sensitivities. A complete and dedicated recasting and reinterpretation procedure would be required to establish such exclusion bounds; however, this is beyond the scope of the present work. Moreover, for benchmark scenarios in which the HNL mixes purely with $ \nu_\mu $ ($ \nu_\tau $) instead of $ \nu_e $, the sensitivities are expected to be comparable (weaker), given the similar (distinct) HNL decay products.

Finally, we note that our analysis is readily generalizable to additional theoretical scenarios, such as a scenario with both a new scalar and a dark photon, and the dark axion portal scenario.

In this work, we considered an ALP with off-diagonal couplings to the charm and up quarks and to a pair of long-lived Majorana HNLs. We performed both fast and dedicated MC simulations to assess the search prospects for long-lived HNLs at the proposed Super Tau–Charm Facility.

We focused on the center-of-mass energy mode of $ \sqrt{s}=3.773 $ GeV at STCF and on the decays $ D^+\to \pi^+ a $ and $ a\to NN $. The ALP a is on shell and decays promptly and exclusively into a pair of N's. Using both a simple phenomenological analysis and a dedicated reconstruction-level analysis, we determined the sensitivity reach for the active-neutrino–HNL mixing parameter $ |V_{eN}|^2 $ as a function of the HNL mass, with the branching ratio of $ D^+\to \pi^+ a $ fixed at the experimentally allowed upper bound by assuming a sufficiently small value of $ g^u_{2,1}/\Lambda $.

Our numerical results show that STCF can probe values of $ |V_{eN}|^2 $ down to approximately one to two orders of magnitude below the existing bounds, while an approximate reinterpretation of a CHARM search is found to surpass the STCF sensitivities. A complete and dedicated recasting and reinterpretation procedure would be required to establish such exclusion bounds; however, this is beyond the scope of the present work. Moreover, for benchmark scenarios in which the HNL mixes purely with $ \nu_\mu $ ($ \nu_\tau $) instead of $ \nu_e $, the sensitivities are expected to be comparable (weaker), given the similar (distinct) HNL decay products.

Finally, we note that our analysis is readily generalizable to additional theoretical scenarios, such as a scenario with both a new scalar and a dark photon, and the dark axion portal scenario.