Sensitivity Study of the Tau Lepton Electric Dipole Moment at the Super Tau-Charm Facility

The electric dipole moment (EDM) is a physical quantity that characterizes the distribution of electric charge within a system. Theoretically, the EDM of a point-like particle is expected to be zero. Within the Standard Model (SM), due to charge-parity-conjugation (CP) violation, the intrinsic EDM of the τ lepton is predicted to be approximately $ {10}^{-37}e\,\mathrm{cm} $ [1−3], far below the experimental sensitivity. From the perspective of symmetry, a non-zero EDM signifies CP violation, which is a crucial prerequisite for the observed matter-antimatter asymmetry in the current observable universe. The CP violation predicted by the Standard Model is minute and insufficient to account for the observed asymmetry [4]. However, some Beyond the Standard Model (BSM) theories predict a significantly larger value, around $ {10}^{-19}e\,\mathrm{cm} $ [5, 6], making it feasible to measure the τ EDM experimentally and potentially reveal new physics.

For stable particles like neutrons and electrons, EDM measurements typically employ precise experimental methods based on the spin precession phenomenon in strong electric or magnetic fields [7, 8]. When a particle with a non-zero EDM interacts with an electric field, its spin direction will undergo a slight deflection, which can be detected by highly sensitive experimental apparatuses. To date, numerous experiments have set extremely stringent upper limits on the EDMs of neutrons and electrons, further corroborating the predictions of the SM while also providing crucial experimental constraints for new physics searches beyond the SM.

For short-lived particles such as heavy quarks and leptons, direct EDM measurements face significant challenges due to their extremely short lifetimes, which generally preclude traditional spin precession experiments [9]. Consequently, indirect approaches are adopted, wherein deviations in scattering cross-sections or decay rates from SM predictions are used to search for potential EDM signals of these particles. This approach not only broadens the applicability of EDM measurements but also opens new experimental avenues for exploring BSM physics [10].

The most precise measurement of the τ EDM to date comes from the Belle experiment in Japan, which studied the process $ e^+e^-\rightarrow \tau^+\tau^- $ at the KEKB collider, setting an upper limit of $ {10}^{-17}e\,\mathrm{cm} $ on the τ EDM [11, 12]. To further explore the potential of measuring the τ EDM in the tau-charm energy region, we are considering the next-generation electron-position collider, the Super Tau-Charm Facility (STCF), which is under design and construction in China. The STCF is designed to operate in the energy range of $ 2.0-7.0\,\mathrm{GeV} $ with a peak luminosity of $ 0.5\times{10}^{35}\,\mathrm{cm^{-2}s^{-1}} $, allowing an integrated luminosity of $ 1\,\mathrm{ab^{-1}} $ per year. At center-of-mass energies (CMEs) of $ 4.2\,\mathrm{GeV} $ and $ 7.0\,\mathrm{GeV} $, the STCF is expected to produce $ 3.5\times{10}^9 $ and $ 1.7\times{10}^9 $ τ pairs per year, respectively, significantly enhancing the sensitivity of τ EDM measurements. Additionally, the STCF will utilize a vertex detector to improve vertex resolution [13]. Compared to the high-energy LEP experiments and the $ 10\,\mathrm{GeV} $ Belle experiment, the tau-charm energy region offers reduced radiative τ pair events, higher τ pair production cross-sections, and more accurate charged particle identification and photon reconstruction, leading to a higher reconstruction efficiency for τ leptons.

The Lagrangian related to the EDM of τ lepton is given by [14, 15]:

Including the EDM term, the scattering amplitude for the tau-pair production at electron-positron collider $ e^{+}(p_1) e^{-}(p_2) \rightarrow \tau^{+}\left(p_3,s_+\right) \tau^{-}\left(p_4,s_-\right) $ is given by:

where $ \mathcal{M}_{\rm{SM}} $ is the contribution from the SM without the EDM and $ \mathcal{M}_{\rm{Re}} $ is the contribution from the EDM operator. Note that we introduce extra v in the second term in order to balance the dimensionality of $ \mathcal{M}_{\rm{SM}} $ and $ \mathcal{M}_{\rm{Re}} $. For simplicity, we only keep the real part of $ {\rm{Re}}(d_\tau) $. Then the squared matrix element contains three contributions:

where $ |\mathcal{M}_{\rm{inter}}|^2=2{\rm{Re}}(\mathcal{M}_{\rm{SM}}^*\mathcal{M}_{\rm{Re}}) $ is the interference between the SM contribution and the EDM operator and will be the key component for retrieving the information about the EDM of τ lepton. Each term in the matrix element depends on the spin of $ \tau^+ $ and $ \tau^- $:

The optimal observable [16] is used to maximize the sensitivity to $ d_\tau $ which is given by:

The mean value of the observable $ \mathcal{O}_{\mathrm{Re}} $ is given by:

where the integration is performed over the available phase space Π and we ignore the higher order terms ($ |\mathcal{M}_{\rm{Re}}|^2 $). The mean value of the optimal observable is thus a linear function of $ {\rm{Re}}(d_\tau)v $ [11, 14]:

where

Therefore, the electric dipole moment $ d_\tau $ can be obtained from $ |\mathcal{M}_{\mathrm{SM}}|^2 $ and $ |\mathcal{M}_{\operatorname{inter}}|^2 $ of which the computation depends on the τ lepton momenta and spin vectors, which will be further constructed from the momenta of the tau decay products.

The $ \rho\rho $ mode of τ pair production, $ e^+e^-\rightarrow \tau^+\tau^- \left(\tau^+\rightarrow \pi^+\pi^0\bar{\nu}_\tau,\tau^-\rightarrow \pi^-\pi^0\nu_\tau,\pi^0\rightarrow \gamma\gamma\right) $, is the dominant channel and provides the most promising avenue for probing the intrinsic EDM of the τ lepton [11]. Consequently, this study focuses on this process as the signal process. By simulating the electron-positron annihilation process, we employ multivariate analysis to optimize selection criteria, filter signal events, and ensure a low background rate. We analyze the kinematic properties of the final state particles, compute the momentum and spin of the τ lepton, and obtain the optimal observables and their relationship with the EDM, facilitating the measurement of the τ EDM at the STCF.

The STCF detector system is a sophisticated assembly designed to maximize physics potential in the τ-charm energy region. From the interaction point outward, it features a tracking system comprising an inner tracker (ITK) using radiation-hard technologies like μRWELL-based gaseous detectors or MAPS-based silicon pixels, followed by a large main drift chamber (MDC) with helium-based gas for precise momentum measurement. Particle identification is achieved via a barrel RICH detector and an endcap time-of-flight system (DTOF), providing kaon-pion separation up to $ 2\,\mathrm{GeV/c} $. The electromagnetic calorimeter (EMC) uses pure CsI crystals for high-resolution photon detection and energy measurement. A superconducting solenoid generates a $ 1\,\mathrm{T} $ magnetic field for tracking, surrounded by an iron yoke for structural support and flux return. The outermost layer is a muon detector (MUD) combining resistive plate chambers and plastic scintillators for efficient muon identification. The system is optimized for high luminosity, with advanced data acquisition handling event rates up to $ 400\,\mathrm{kHz} $.

In this study, we utilize the Monte Carlo (MC) method with the KKMC and Tauola generators [17, 18] to generate $ 5,567,300 $ $ e^+e^-\rightarrow \tau^+\tau^- $ events at a CME of $ 4.68\,\mathrm{GeV} $ based on the theoretical reaction cross-sections and branching ratios, and simulate the detector's electronic signals (including timing, amplitude, etc.) under a fast simulation package [19]. The process $ e^+e^-\rightarrow \tau^+\tau^- \left(\tau^+\rightarrow \pi^+\pi^0\bar{\nu}_\tau,\tau^-\rightarrow \pi^-\pi^0\nu_\tau\right) $ accounts for approximately (25%)²=6.2% (i.e., the signal fraction). The main background decays include $ \tau^{\pm} \xrightarrow{17.8\%} \nu_{\tau} \mathrm{e}^{\pm} {\nu}_{\mathrm{e}}, \tau^{\pm} \xrightarrow{17.4\%} \nu_{\tau} \mu^{\pm} {\nu}_{\mu}, \tau^\pm \xrightarrow{10.8\%} \pi^\pm \nu_\tau, \tau^\pm\xrightarrow{9.3\%}\pi^\pm \pi^0\pi^0 \nu_\tau $ [20]. The reconstruction software is then used to derive physical quantities such as the momentum and energy of the final state particles from these electronic signals.

The MC simulated data includes truth information (such as particle species) that is not available in actual experiments. This additional information helps optimize the parameter selection in the analysis algorithms, which in turn enhances the effectiveness of these algorithms when applied to real experimental data.

In the final state of the studied process, the detectable particles are $ \pi^+ $, $ \pi^- $, and γ. The analysis algorithm filters events by reconstructing tracks in the main drift chamber (MDC), selecting those with exactly two charged tracks (corresponding to $ \pi^+ $ and $ \pi^- $) and a total charge sum of zero. It is important to note that only tracks within a small distance from the electron-positron collision point are considered valid to eliminate cosmic rays and beam-related backgrounds. Further, the analysis algorithm uses Particle Identification (PID) based on ionization energy loss ($ \mathrm{d}E/\mathrm{d}x $) and Time of Flight (TOF) to identify $ \pi^+ $ and $ \pi^- $ particles. Events are then selected if they contain exactly one $ \pi^+ $ and one $ \pi^- $. This step can efficiently suppress background events where the final state contains e or μ.

For each photon, the algorithm scans all charged tracks to find the minimum angle between the photon and the tracks. A photon is considered valid only if this minimum angle exceeds 20° and its energy is greater than $ 0.04\,\mathrm{GeV} $. This criterion helps to reject noise photons produced by hadronic showers. The expected number of final state photons in signal events is $ 4 $ from two $ \pi^0 $s. Considering the presence of noise photons in the electromagnetic calorimeter (EMC), the algorithm selects events with at least $ 4 $ photons. This will suppress events without $ \pi^0 $ in the τ decay.

To further eliminate noise photons, a machine learning approach is employed to distinguish signal photons from noise photons. The variables used in the machine learning analysis include:

● gam_energy: Energy of the EMC photon cluster (GeV).

● gam_secmom: Second moment of the EMC cluster, describing the shower shape (energy spread).

● gam_a42mom: $ a_{42} $ moment of the EMC cluster, a higher-order shower shape variable.

● gam_hits: Number of crystals hit in the cluster.

● gam_var2: $ (E_{3\times3} - E_{\text{seed}}) / E_{3\times3} $, energy concentration in $ 3\times3 $ array.

● gam_var3: $ (E_{\text{total}} - E_{\text{seed}}) / (N_{\text{hits}} - 1) / E_{\text{seed}} $, average non-seed energy normalized by seed energy.

Their distributions for signal and noise photons are shown in Fig. 1. The distinct differences between signal and noise photons allow for effective photon selection using machine learning. We train models using Boosted Decision Trees (BDT) and Boosted Decision Trees with Gradient Boosting (BDTG) with the TMVA Toolkit [21], and the resulting Receiver Operating Characteristic (ROC) curves are shown in Fig. 2. BDTG performs slightly better than BDT, so BDTG is used in subsequent analysis.

Figure 1.  Variable distributions for signal and noise photons.

Figure 2.  ROC curves for BDT and BDTG for photon selection.

The BDTG response distributions for training and testing datasets are shown in the left panel of Fig. 3, indicating no significant overfitting as the distributions agree quite well. The BDTG cut efficiency is presented in the right panel of Fig. 3. Considering both signal efficiency and signal-to-noise ratio, a BDTG cut of $ 0.6 $ is chosen, leaving approximately 84% of signal photons and 37% of noise photons.

Figure 3.  BDTG response (left panel) and cut efficiency (right panel) for photon selection.

After the machine learning selection, the effective number of photons for $ \rho\rho $ mode events is typically $ 4 $. Consequently, the analysis algorithm selects events with $ 4 $ photons.

While we can reconstruct the tracks of final-state $ \pi^+ $, $ \pi^- $, and photons, it is unclear which photons are signal photons and how these photons associate with the corresponding $ \pi^+ $ or $ \pi^- $ in $ \rho\rho $ mode. The objective of this step is to determine the best pairing of these particles. Various methods were designed and compared for selecting the optimal pairing, with the best approach determined by the $ \chi^2 $ value from a joint kinematic fitting, as detailed in Section III. The joint kinematic fitting not only determines the best pairing but also improves the signal-to-noise ratio by selecting events with $ \chi^2<10 $.

To further enhance the signal-to-noise ratio, event-level machine learning was applied to select signal events. The variables used in this analysis include the momentum of the $ \pi^+ $, $ \pi^- $, and the four photons (in sequence). Similar to the photon selection using machine learning, the BDTG response and cut efficiency are shown in Fig. 4. A BDTG cut of $ 0.2 $ was chosen, resulting in about 73% of the signal events and 32% of the background events remaining.

Figure 4.  BDTG response (left panel) and cut efficiency (right panel) for event selection.

To obtain the optimal observables, it is necessary to reconstruct for the momentum of the τ lepton, as detailed in Section IV A. During this process, scenarios may arise with two solutions, one solution, or no solution. Events with no solution are considered background and are excluded.

The selection algorithms applied reduced the initial $ 5,567,300 $ events to $ 27,125 $ events. The selection efficiency at each step is summarized in Table 1. The overall selection efficiency is 0.49%, and the signal ($ \tau^+\rightarrow \pi^+\pi^0\bar{\nu}_\tau,\tau^-\rightarrow \pi^-\pi^0\nu_\tau $) efficiency is 6.3%, with the signal purity increased from 6.2% to 80.0%.

No.	Step	Percentage of previous step		Signal purity
		Inclusive events	Signal events
0	Total events	-	-	6.2%
1	Number of charged tracks = 2, total charge = 0	58.3%	76.5%	8.1%
2	Number of photons = 4	7.2%	23.7%	26.7%
3	Number of $ \pi^+ $ = 1, Number of $ \pi^- $ = 1	81.8%	92.1%	30.0%
4	Passed the particle pairing	25.2%	52.3%	62.5%
5	Passed event-level machine learning selection	57.9%	73.4%	79.3%
6	Passed the τ momentum reconstruction	97.0%	97.7%	80.0%

Table 1.  Event selection results for $ 5,567,300 $ simulated events.

Detailed event type analysis on the $ 27,125 $ selected events with a generic topology analysis package, TopoAna [22], shows that background decays $ \tau^{\pm} \rightarrow \nu_{\tau} \mathrm{e}^{\pm} {\nu}_{\mathrm{e}}, \tau^\pm\rightarrow \pi^\pm \nu_\tau $ have been almost entirely filtered out. The dominant background processes after selection involves more $ \pi^0 $, such as $ \tau^\pm\rightarrow \pi^\pm \pi^0\pi^0 \nu_\tau $ (about 14%), which may require further selection optimization.

Considering that at the CME of $ 4.2\,\mathrm{GeV} $, the STCF is expected to produce approximately $ 3.5 \times 10^9 $ τ lepton pairs per year [13], of which around $ 2.2 \times 10^8 $ are signal events of the $ \rho\rho $ mode, we estimate a signal yield of $ 1.4 \times 10^7 $ per year after selection, with a signal efficiency of 6.3% and a signal purity of 80.0%. For comparison, up to 2022, Belle experiments have $ 5.2 \times 10^7 $ events of the $ \rho\rho $ mode. After selection, the signal yield is $ 3.3 \times 10^6 $, with a signal efficiency of 6.3% and a signal purity of 82.4% [11]. The current feasibility study shows that the tau pair selection efficiency is comparable to that of the B factory. After 10 years of operation, the STCF will collect $ 1.4 \times 10^8 $ tau pairs of $ \rho\rho $ mode after reconstruction, which is $ 2 $ orders of magnitude higher than Belle. In addition, the event selection at STCF can be further optimized with the vertex detector. This increase in statistics will improve the sensitivity of the τ lepton EDM measurement.

This section introduces and compares four methods for determining the pairings of $ \pi^+, \gamma_{(1)}, \gamma_{(2)} $ and $ \pi^-, \gamma_{(3)}, \gamma_{(4)} $. The main method we use is joint kinematic fitting.

Let the number of photons detected in the electromagnetic calorimeter (EMC) for the current event be n ($ n \geqslant 4 $, and $ n=4 $ in this study). The total number of possible pairings is $ C_n^2 C_{n-2}^2 $. The analysis algorithm performs kinematic fitting on all possible pairings. The kinematic constraints are imposed based on the conservation laws, including the total energy-momentum conservation and the mass constraints of the intermediate states τ and $ \pi^0 $. The specific equations are as follows:

where:

Eqs. (11) and (12) are automatically satisfied by the reconstructed data, while Eq. (13) imposes the $ \pi^0 $ mass constraint on the photon energy-momentum. Eqs. (9), (10), and (14) constrain the unknown neutrino (ν) energy-momentum. The kinematic fitting adjusts the known quantities, solves for the unknowns, and yields a fitting error $ \chi^2 $ (the deviation of fitted and original values over the error). The correct pairing will have a smaller $ \chi^2 $ and a smaller $ p_{\nu}=p_{\nu_{(1)}} + p_{\nu_{(2)}} $, whereas incorrect pairings lead to situations where the constraint equations have no solutions or yield non-physical solutions, resulting in a larger $ \chi^2 $ and a larger $ p_{\nu} $.

From the truth information in the simulated data, it is observed that the sum of the magnitudes of the momenta of the two neutrinos is generally less than $ E^\text{Total}/2 $. Therefore, it can be required that the kinematic fitting solution satisfies:

Under this constraint, the algorithm selects the pairing with the smallest $ \chi^2 $ as the chosen pairing. The kinematic fitting can also be used to reject events that fail the fit or have $ \chi^2 \geqslant 10 $ or $ p_{\nu}\geqslant E^\text{Total}/2 $, thus improving the signal-to-noise ratio.

The truth information in the simulated data contains the actual pairing information of photons with $ \pi^+, \pi^- $. By comparing selected pairings with this information (pairings with relative and absolute errors of each photon momentum component in order within an acceptable range are considered correct), the pairing correct rate can be obtained. Additionally, the correctness of the photon pairing forming the $ \pi^0 $ can be examined, that is, only checking whether $ \gamma_{(1)}, \gamma_{(2)} $ are paired to form a $ \pi^0_{(1)} $ and whether $ \gamma_{(3)}, \gamma_{(4)} $ are paired to form a $ \pi^0_{(2)} $, regardless of the pairing selection of $ \pi^0_{(1)}, \pi^0_{(2)} $ with $ \pi^+, \pi^- $. For signal events, the pass rate and pairing correct rate of the kinematic fitting are shown in Table 2, with a fully correct particle pairing rate of up to 82.5%.

Signal events	Percentage
Passed the joint kinematic fitting	52.3%
Correct $ \pi^0 $ pairing	95.6%
Fully correct particle pairing	82.5%

Table 2.  Kinematic fitting pass rate and pairing correct rate.

We examine the invariant mass distributions of $ (\pi^+, \pi^0_{(1)}) $ and $ (\pi^-, \pi^0_{(2)}) $, as shown in Fig. 5. It can be seen that the invariant masses are concentrated around the mass of the $ \rho(770) $, as expected. Besides, we compare the distribution of $ \chi^2 $ and $ p_{\nu}=p_{\nu_{(1)}} + p_{\nu_{(2)}} $ for selected pairings and all pairings in Fig. 6, which also meet expectations.

Figure 5.  Invariant mass distribution of $ (\pi^+,\pi^0_{(1)}) $ (left panel) and $ (\pi^-,\pi^0_{(2)}) $ (right panel).

Figure 6.  Distribution of $ \chi^2 $ (left panel) and $ p_{\nu}=p_{\nu_{(1)}} + p_{\nu_{(2)}} $ (right panel) for selected pairings and all pairings.

We tested three additional methods for particle pairing (detailed in Appendix A): stepwise kinematic fitting, kinematic fitting with ρ resonance mass constraint, and kinematic fitting with momentum direction. Among the methods tested, joint kinematic fitting achieves the highest accuracy. Therefore, all subsequent steps are based on the data selected using this method.

Due to the presence of neutrinos among the τ lepton decay products, which cannot be detected by the instruments, it is impossible to fully reconstruct the final-state particles to determine the momentum of the τ lepton. Therefore, a special method is required to calculate the τ lepton momentum.

After event selection and particle pairing, the momenta and energies of $ \pi^+, \pi^-, \pi^0_{(1)} $, and $ \pi^0_{(2)} $ have been obtained from the track information of the main drift chamber and the photon signals in the electromagnetic calorimeter. To accurately determine the τ lepton momentum, analytical computation is used to derive the solutions.

In electron-positron collider experiments, the beams are designed to collide at a small angle at the interaction point to optimize collision performance and data collection efficiency. In the laboratory frame, the total momentum of the system is given by $ \boldsymbol{p}^\text{Total} = (p_x, 0, 0) $, where $ p_x $ is a small value that can be determined from the total energy $ E^\text{Total} $ and the design parameters of the experiment. By applying a Lorentz transformation to all particles to move into the center-of-mass frame, we have:

and the energies of the $ \tau^+ $ and $ \tau^- $ leptons become:

The magnitude of the $ \boldsymbol{p}_{\tau^\pm} $ can then be derived from $ E_{\tau} $ as follows:

Next, we solve for the direction of $ \boldsymbol{p}_{\tau^\pm} $. By combining Eqs. (9), (10), (14), and (17), the cosine of the angle $ \theta_{\pm} $ between the vectors $ \boldsymbol{p}_{\tau^\pm} $ and $ \boldsymbol{p}_{h^\pm} = \boldsymbol{p}_{\pi^\pm} + \boldsymbol{p}_{\pi^0} $ is given by [23]:

where:

In the center-of-mass frame, we have:

Given that the angle between the vectors $ \boldsymbol{p}_{\tau^+} $ and $ \boldsymbol{p}_{h^+} $ is $ \theta_+ $, and the angle between the vectors $ -\boldsymbol{p}_{\tau^-} $ and $ -\boldsymbol{p}_{h^-} $ is $ \theta_- $, $ \boldsymbol{p}_{\tau^+} $ must lie on the intersection of two cones: one with axis $ \boldsymbol{p}_{h^+} $ and half-angle $ \theta_+ $, and the other with axis $ -\boldsymbol{p}_{h^-} $ and half-angle $ \theta_- $. Using Mathematica, two analytic solutions for $ \boldsymbol{p}_{\tau^+} $ were obtained. There are square root operations in the expressions for the solutions. If the term under the square root is negative (small negative values due to experimental uncertainties are treated as zero in this study), the corresponding case has no physical solution and is therefore discarded. If the term is positive, it yields two distinct solutions. Among these, only one corresponds to the physical reality. Experimentally, this ambiguity cannot be resolved, so following the same approach as the Belle experiment [11], this study adopts the method of taking the average of the two solutions.

After determining $ \boldsymbol{p}_{\tau^\pm} $ in the center-of-mass frame, a Lorentz transformation is applied to convert it back to the laboratory frame, resulting in the reconstructed momentum of the τ lepton. For signal events, the relative deviation of the reconstructed transverse momentum $ p_T = \sqrt{p_x^2 + p_y^2} $ and the angular distribution between the reconstructed and truth τ lepton momentum are shown in Fig. 7. The full width at half maximum (FWHM) of the relative deviation in transverse momentum is $ 0.10 $, and the peak opening angle is $ 7^\circ $ with a root mean square (RMS) of $ 10^\circ $.

Figure 7.  Relative deviation of the transverse momentum (left panel) and the direction (right panel) of the τ lepton.

The decay products keep the information of the spin of the tau lepton, which can be used to reconstruct the polarimeter vector. For the hadronic decay channel of the tau lepton, we have the decay matrix element as [18]:

with:

where $ p_\pm $ are the momenta of $ \tau^\pm $, $ N^\pm $ are the momenta of the corresponding neutrinos. J will be different for different decay channels:

Combining the decay and the production, we obtain for the each component of the produciton matrix element:

Then the optimal observable can be constructed by all the final state momentum. However, as discussed above, we can only determine the momenta of the neutrino/tau lepton up to a two-fold ambiguity. Then the matrix element calculated from these two solutions will be averaged when calculating the optimal observable.

MadGraph [24] is used along with custom UFO model files to simulate the production of τ lepton pairs for different values of $ d_\tau $. Following the reconstruction of the final states in previous sections, a fast simulation of the detector response is performed with Delphes [25]. The distributions of the optimal observables for several different choices of $ d_\tau $ from the $ \rho\rho $ channel are shown in the left panel of Fig. 8. The mean values of the optimal observables are calculated and a linear fit is performed based on Eq. (7), which is shown in the right panel of Fig. 8. The corresponding result of the linear fit is given in Table 3 together with the estimated error on the measurement of the mean value of the optimal observables. The estimated sensitivity after 10 years of operation at the STCF is $ |d_\tau| < 3.89\times 10^{-18}\,e\cdot\mathrm{cm} $ at 68% confidence level. This provides a basis for estimating the EDM of the τ lepton in future experiments at STCF.

Figure 8.  Left panel: Distribution of $ \mathcal{O}_{\operatorname{Re}} $ for different $ d_{\tau} $ (Red: $ d_{\tau}<0 $; Green: $ d_{\tau}>0 $). The ratio r is defined as the value of $ \dfrac{1}{\sigma}\dfrac{d\sigma}{d\mathcal{O}_{\operatorname{Re}}} $ for nonzero $ d_{\tau} $ divided by that for $ d_{\tau}=0 $. Right panel: Relationship between $ \left<\mathcal{O}_{\operatorname{Re}}\right> $ and $ d_{\tau} $.

Mode	$ a_{\rm{Re}}\,[(10^{-18}{e\cdot \rm cm})^{-1}] $	$ b_{\rm{Re}} $	$ \delta\langle \mathcal{O}_{\rm{Re}}\rangle $
$ \rho\rho $	$ 3.06\times10^{-5}\pm1.43\times 10^{-7} $	$ -1.85\times10^{-5}\pm 4.87\times 10^{-5} $	$ 1.088\times 10^{-4} $

Table 3.  The result of the linear fit on the mean values of the optimal observables.

In this study, the process of τ lepton pair production in $ e^+e^- $ collisions was simulated using Monte Carlo methods. Signal photons and signal events were effectively selected through machine learning techniques, significantly improving the signal-to-noise ratio. Based on the simulation data, the event selection algorithm was optimized and its parameters adjusted, resulting in a final signal event selection efficiency of approximately 6.3%, with the signal purity increased from 6.2% to 80.0%. Regarding particle pairing selection, four methods were compared in terms of their principles and results. A novel approach using joint kinematic fitting was introduced, which greatly enhanced the accuracy of particle pairing, achieving a correct pairing rate of 82.5% for signal events. Additionally, analytical computation was employed to derive the solutions for the τ lepton momentum, leading to a more precise determination of the number of solutions and their numerical accuracy. Subsequently, the polarimeter vector of the τ lepton and the squared spin density matrix were calculated. The optimal observables and their relationship with the EDM were obtained, with the estimated sensitivity of $ |d_\tau| < 3.89\times 10^{-18}\,e\cdot\mathrm{cm} $ at 68% confidence level, laying the groundwork for determining the τ EDM from experimental data in future STCF experiments.

We thank Zhipeng Xie, Yupeng Pei, Zekun Jia, Bo Wang, and Mingyi Liu for valuable discussions.

This method has two steps. First, $ \gamma_{(1)},\gamma_{(2)} $ are paired into $ \pi^0_{(1)} $, and $ \gamma_{(3)},\gamma_{(4)} $ into $ \pi^0_{(2)} $, followed by a kinematic fitting (using the constraints in Eqs. (11) and (13)). Considering the equivalence of $ \pi^0_{(1)} $ and $ \pi^0_{(2)} $, there are $ C_n^2 C_{n-2}^2/2 $ possible pairings. The algorithm traverses all pairings and selects the one with the smallest sum of $ \chi^2 $ from the kinematic fitting. Second, $ \pi^0_{(1)} $ is paired with $ \pi^+ $ and $ \pi^0_{(2)} $ with $ \pi^- $ for another kinematic fitting (with the constraints in Eqs. (9), (10), (12), and (14)); the reverse pairings are also attempted, with the fit yielding the smaller $ \chi^2 $ being selected (while still requiring that the neutrino momentum satisfies Eq. (16)). The particle pairing accuracy in this method is lower than Method I (joint kinematic fitting) due to the two-step kinematic fitting process, where separate $ \chi^2 $ determinations lead to a sequence of locally optimal choices that may not be globally optimal.

In τ decays ($ \tau^+ \rightarrow \pi^+ \pi^0 \nu, \tau^- \rightarrow \pi^- \pi^0 \nu, \pi^0 \rightarrow 2\gamma $), $ \pi^\pm $ and $ \pi^0 $ usually form a ρ resonance [20], allowing the inclusion of the ρ mass constraint:

Given the broad mass range of the ρ resonance (indicating a large uncertainty), directly adding this constraint to the kinematic fitting would force the fit to adjust physical quantities to strictly satisfy the constraint, which is not desirable. Instead, using the $ \pi^0_{(1,2)} $ energies and momenta obtained from the kinematic fitting in the second step of Method II, one can calculate the corresponding ρ masses using Eq. (A1), compare them to the theoretical mass, and select the pairing with the smallest deviation. Using this method, the particle pairing accuracy is still lower than Method I.

In this method, after performing the kinematic fitting for all possible pairings, the neutrino momenta $ \boldsymbol{p}_{\nu_{(1,2)}} $ are obtained, and the τ lepton momenta are calculated as $ \boldsymbol{p}_{\tau^+} = \boldsymbol{p}_{\pi^+} + \boldsymbol{p}_{\gamma_{(1)}} + \boldsymbol{p}_{\gamma_{(2)}} + \boldsymbol{p}_{\nu_{(1)}} $ and $ \boldsymbol{p}_{\tau^-} = \boldsymbol{p}_{\pi^-} + \boldsymbol{p}_{\gamma_{(3)}} + \boldsymbol{p}_{\gamma_{(4)}} + \boldsymbol{p}_{\nu_{(2)}} $. For the correct pairing, the angle between the momenta of the final-state particles and the τ lepton momenta (considering only the acute angle) should be smaller. Let $ \boldsymbol{n}_T $ denote the unit vector along $ \boldsymbol{p}_{\tau^+} $, and define:

Larger T values are more likely to correspond to the correct physical scenario [26], so this method selects the pairing with the largest T value as the chosen pairing. Results indicate that this method yields a higher pairing accuracy than Methods II and III, but it is still inferior to Method I.