Prospect for measurement of the CP-violating phase $\phi_{s}$ in the $B_{s}\rightarrow J/\psi\phi$ channel at a future $Z$ factory

The CP-violating phase $\phi_{s}$, $B_{s}$ decay width $\Gamma_{s}$, and the width difference $\Delta\Gamma_{s}$ between the heavier and lighter $B_{s}$ meson eigenstates have been thoroughly investigated in experiments conducted by ATLAS  [5, 6], CDF [7], CMS [8, 9], D0 [10], and LHCb [11, 12, 13, 14, 15, 16]. The decay channel $B_{s}\rightarrow J/\psi(\rightarrow\mu^{+}\mu^{-})\phi(\rightarrow K^{+}K^{-})$ is particularly notable due to its sizeable branching fraction and the final state consisting entirely of charged tracks. This decay channel benefits from the narrow decay widths of the $J/\psi$ and $\phi$ particles, effectively suppressing the combinatorial background. It stands as the most prominent channel for measuring $\phi_{s}$, and it also allows for the concurrent extraction of $\Delta\Gamma_{s}$ and $\Gamma_{s}$.

The time and angular distribution of $B_{s}\rightarrow J/\psi\phi$ is a sum of ten terms corresponding to the three polarization amplitudes and the non-resonant S-wave, together with their interference terms:

$$\frac{d^{4}\Gamma(B_{s}\rightarrow J/\psi\phi)}{dtd\Omega}\propto\sum_{k=1}^{10}h_{k}(t)f_{k}(\Omega),$$

(1)

where

$$\begin{array}[]{rcl}h_{k}(t|B_{s})&=&\displaystyle N_{k}e^{-\Gamma_{s}t}\Biggl{[}a_{k}\cosh(\frac{1}{2}\Delta\Gamma_{s}t)+b_{k}\sinh(\frac{1}{2}\Delta\Gamma_{s}t)\\ &+&\displaystyle c_{k}\cos(\Delta m_{s}t)+d_{k}\sin(\Delta m_{s}t)\Biggl{]},\\ h_{k}(t|\bar{B}_{s})&=&\displaystyle N_{k}e^{-\Gamma_{s}t}\Biggl{[}a_{k}\cosh(\frac{1}{2}\Delta\Gamma_{s}t)+b_{k}\sinh(\frac{1}{2}\Delta\Gamma_{s}t)\\ &-&\displaystyle c_{k}\cos(\Delta m_{s}t)-d_{k}\sin(\Delta m_{s}t)\Biggl{]},\end{array}$$

and $f_{k}(\Omega)$ is the amplitude function.

In the formulation of $h_{k}(t)$, the term $\Delta m_{s}$ represents the mass difference between the $B_{s}$ mass eigenstates, while $N_{k}$ denotes the amplitude of the component at $t=0$. The phase $\phi_{s}$ is encapsulated within the parameters $a_{k},b_{k},c_{k}$, and $d_{k}$. For an in-depth explanation of these parameters, one can refer to the LHCb publication [12]. The values of $\phi_{s}$, $\Delta\Gamma_{s}$, and $\Gamma_{s}$ could be obtained by fitting the time and angular distributions of $B_{s}\rightarrow J/\psi\phi$ decay events.

Additionally, when determining $\phi_{s}$, $\Delta\Gamma_{s}$ and $\Gamma_{s}$, parameters such as $\Delta m_{s}$ can also be simultaneously derived from this fitting. However, the precision of these parameters is beyond the scope of this work and will not be discussed here.

The CEPC and the baseline detector (CEPC-v4) [18] are taken as an example to study the precision of $\phi_{s}$, $\Delta\Gamma_{s}$ and $\Gamma_{s}$. As a baseline, the CEPC is assumed to run in the Tera-$Z$ mode, i.e., produces $10^{12}$ $Z$ bosons during its lifetime. The CEPC baseline detector consists of a vertex system, a silicon inner tracker, a TPC, a silicon external tracker, an electromagnetic calorimeter, a hadron calorimeter, a solenoid of 3 Tesla, and a Return Yoke.

A Monte Carlo signal sample is generated to analyze the geometric acceptance and the reconstruction efficiency of the $B_{s}\rightarrow J/\psi\phi$ decay. Additionally, this sample is also used for the examination of the proper decay time resolution for the $B_{s}$, which has a direct correlation with the spatial resolution at the $B_{s}$ decay vertex.

Using the WHIZARD [19] generator, roughly 6000 events of $Z\rightarrow b\bar{b}\rightarrow B_{s}(\bar{B}_{s})+X$ are simulated. The $B_{s}(\bar{B}_{s})$ particles are then forced to decay through the $B_{s}(\bar{B}_{s})\rightarrow J/\psi(\rightarrow\mu^{+}\mu^{-})\phi(\rightarrow K^{+}K^{-})$ process using PYTHIA 8 [20], with a uniform distribution in phase space.

The simulation of particle transport within the detector utilizes MokkaC, the simulation software for the CEPC study, based on the GEANT 4 [21]. Based on Monte Carlo truth data, the reconstructed particles are categorized into hadrons, muons, and electrons.

The $J/\psi$ candidates are reconstructed from every pairing of a positively charged muon with a negatively charged muon. They are then selected based on the invariant mass window, ranging from $3.07$ to $3.14$ $\mathrm{GeV}/c^{2}$. The $\phi$ candidates are reconstructed from every possible combination of a positively charged kaon and a negatively charged kaon. The $\phi$ candidates are selected within the mass window from $1.017$  to $1.023$ $\mathrm{GeV}/c^{2}$. The $B_{s}(\bar{B}_{s})$ meson is reconstructed by combining all pairs of $J/\psi$ and $\phi$ candidates identified in the preceding steps. The four-momentum of the $B_{s}(\bar{B}_{s})$ meson is determined using the four-momentum of the $J/\psi$ and $\phi$, ensuring conservation of energy and momentum. They are selected within a mass window ranging from $5.28$ to $5.46$ $\mathrm{GeV}/c^{2}$. Following the reconstruction of the $B_{s}(\bar{B}_{s})$ meson, a decay vertex is constructed using the tracks associated with the $B_{s}(\bar{B}_{s})$.

As the CEPC was initially designed as a Higgs factory, the secondary vertex reconstruction algorithm and the flavor tagging algorithm are not in the standard CEPC software chain. A vertex reconstruction procedure and a simple flavor tagging algorithm were specially implemented for this study, which will be described in Sect. 2.5.

An additional sample of $Z\rightarrow b\bar{b}\rightarrow X$ is generated to show that a low background level is achievable through appropriate event selection criteria for the $\phi_{s}$ measurement. The detector simulation and event reconstruction processes are consistent with those applied to the signal sample.

Assuming that all $b\bar{b}$ events can be selected with high purity, the background events in this work are the $b\bar{b}$ events that do not contain $B_{s}\rightarrow J/\psi(\rightarrow\mu^{+}\mu^{-})\phi(\rightarrow K^{+}K^{-})$ signal. The branching fraction of $b\bar{b}$ hadronized to $B_{s}$ is 10% [22]. The branching ratio of $B_{s}\rightarrow J/\psi\phi$ is $1.08\times 10^{-3}$. And the branching ratio of $J/\psi\rightarrow\mu^{+}\mu^{-}$, $\phi\rightarrow K^{+}K^{-}$ are $6\%$ and $50\%$ respectively [22]. If the background is not suppressed by any event selection criteria, the number of background events is $1/(10\%\times 1.08\times 10^{-3}\times 6\%\times 50\%)=3.1\times 10^{5}$ times larger than the number of signal events. Applying the invariant mass selection criteria described in Sect. 2.2 to the background sample, the probability of reconstructing a fake $B_{s}$ candidate from a $Z\rightarrow b\bar{b}\rightarrow X$ event is $6.7\times 10^{-6}$. Therefore, after the event selection, background statistics are of the same magnitude as the signal statistics.

The combinatorial background events that pass the invariant mass selection criteria are further suppressed by using vertex information. In the background events, the fake $B_{s}$ candidates come from four arbitrarily combined tracks, two of which are lepton tracks and two of which are hadron tracks. The lepton usually has a large impact parameter, and the hadron has a small one. It is difficult to reconstruct a high-quality vertex with arbitrarily combined tracks. The $D_{xy}^{2}$ is used to measure the quality of the vertex reconstruction, where

$$D_{xy}^{2}=\sum_{\mathrm{tracks}}d_{xy}^{2}.$$

The $d_{xy}$ in the formula represents the distance from the reconstructed vertex to the track in the plane perpendicular to the beam direction. The $D_{xy}^{2}$ distributions of the signal and background are shown in Fig. 1. The vertex $D_{xy}^{2}$ of signal is usually very small. And the $D_{xy}^{2}$ of background is distributed over an extensive range. With a very loose cut at $D_{xy}^{2}<0.1~{}\mathrm{mm^{2}}$, $95\%$ of the signals are selected and $99.2\%$ of the backgrounds are discarded.

By employing a combination of invariant mass and vertex cut, the acceptance $\times$ efficiency of the signal is 75%, while the background is maintained at 1% of the signal level.

Due to potential particle misidentification, a small peaking background may be present in the signal region. Implementing a strict threshold on the hadron ID can reduce this peaking background; however, it would also result in diminished efficiency. This loss of efficiency is not considered in the present analysis due to the excellent PID performances of the CEPC.

At CEPC, the electron tracking performance is as good as that of muon tracking. The $J/\Psi$ could be reconstructed via the $J/\Psi\rightarrow e^{+}e^{-}$ channel as well. Consequently, the total effective sample size is considered to be roughly twice as much as when only the $J/\Psi\rightarrow\mu^{+}\mu^{-}$ decay channel is considered.

The precision of $\phi_{s}$ is affected by the inaccurate determination of the decay time. The proper decay time for the $B_{s}$ meson is determined using the decay vertex position and the transverse momentum of the $B_{s}$ as follows:

$$t=\frac{m_{s}r}{p_{\mathrm{T}}},$$

(5)

where $r=\sqrt{x^{2}+y^{2}}$ represents the $B_{s}$ decay vertex position in the transverse plane, $p_{\mathrm{T}}$ denotes the transverse momentum of the $B_{s}$ meson, and $m_{s}$ represents the mass of the $B_{s}$. The $B_{s}$ decay vertex is constructed from the four tracks from the $B_{s}$ decay. It is assumed that the primary vertex, the production point of the $B_{s}$, is located at the origin. The resolution for determination of the primary vertex is considered negligible, given that the abundance of tracks available to reconstruct the primary vertex far exceeds those available for determining the $B_{s}$ decay location.

The decay point of the $B_{s}$ meson is determined by minimizing the $\chi^{2}$ value, which is calculated as the sum of the squares of the shortest distances from the vertex to each of the four helical tracks. This technique is an adaptation of the method described in Ref. [25].

To simplify the minimization process, the helical track is approximated as a straight line that is tangent to the helix at point $P^{(0)}$. This point $P^{(0)}$ is the closest to a reference point on the helix. In this algorithm, the true decay vertex of $B_{s}$ is used as the reference point, denoted by $v_{\text{ref}}$.

Subsequently, each of the tracks is parameterized by a point $r_{i}$ and a direction $a_{i}$. Consequently, the minimization process is replaced by the solving of a matrix equation

$$v=v_{\mathrm{ref}}+(\sum H_{i})^{-1}\cdot\sum r_{i},$$

(6)

where $v$ is the $B_{s}$ decay point position and

$$H=\begin{pmatrix}a_{y}^{2}+a_{z}^{2}&&-a_{x}a_{y}&&-a_{x}a_{z}\\ -a_{x}a_{y}&&a_{x}^{2}+a_{y}^{2}&&-a_{y}a_{z}\\ -a_{x}a_{z}&&-a_{y}a_{z}&&a_{x}^{2}+a_{y}^{2}\\ \end{pmatrix}.$$

The difference between the reconstructed $B_{s}$ decay vertex position ($r_{\mathrm{reco}}$) and the truth vertex position ($r_{\mathrm{sim}}$) is shown in Fig. 3.

The transverse momentum distribution of $B_{s}$ mesons is shown in Fig. 4. The majority of these mesons have large transverse momentum because they hadronize from a high-energy $b$-quark, which carries nearly half of the beam energy. Consequently, a substantial transverse momentum corresponds to a large Lorentz boost factor, which substantially enhances the resolution of the proper decay time.

Figure 5 shows the distribution of the differences between $t_{\text{reco}}$ and $t_{\text{sim}}$ for the same event. Both $t_{\text{reco}}$ and $t_{\text{sim}}$ represent the proper decay time as calculated from Eq. (5). To derive $t_{\text{sim}}$, detector effects are not considered, and the vertex and $p_{\mathrm{T}}$ are obtained directly from the Monte Carlo truth record. In contrast, $t_{\text{reco}}$ is determined after the particles undergo detector simulation, as well as track and vertex reconstruction, using the reconstructed vertex position and transverse momentum to calculate the proper decay time. The distribution is fitted using the sum of three Gaussian functions with the same mean value. The effective time resolution is combined as

$$\sigma_{t}=\sqrt{-\frac{2}{\Delta m_{s}^{2}}\ln(\sum_{i}f_{i}e^{-\frac{1}{2}\sigma_{i}^{2}\Delta m_{s}^{2}})},$$

where $f_{i}$ and $\sigma_{i}$ are the fraction and width of the i-th Gaussian function. The effective resolution of the decay time, achieved through the combined effect of the three Gaussian resolution models, is ${4.7~{}\mathrm{fs}}$.

The possible impact of non-uniform decay time-dependent efficiency, known as decay time acceptance, on the precision of $\Gamma_{s}(\Delta\Gamma_{s})$ measurement has been evaluated. While it was considered that different time acceptance profiles at hadron and lepton colliders might markedly influence the precision of $\Gamma_{s}$ and $\Delta\Gamma_{s}$ measurements, our findings indicate that the effect is not substantial.

Figure 6 shows the reconstruction efficiency of the $B_{s}$ meson as a function of its proper decay time at the CEPC. The efficiency decreases at larger decay times due to the tracks of the $B_{s}$ with larger flight distances deviating increasingly from the interaction point, which complicates the reconstruction process.

The efficiency is parameterized with a 2nd-order polynomial function $f_{\mathrm{acc}}^{\mathrm{CEPC}}=0.83-0.0061\times t-5.25\times 10^{-5}\times t^{2}$. The time acceptance profile from LHCb, as referenced in Ref. [26] is $f_{\mathrm{acc}}^{\mathrm{LHCb}}=f_{1}\times f_{2}$, where

$$f_{1}=1-0.037\times t+0.001\times t^{2},$$

and

$$f_{2}=\frac{[1.589\times(t-0.097)]^{1.150}}{1+[1.589\times(t-0.097)]^{1.150}}.$$

A toy Monte Carlo simulation is employed to explore the impact of the time acceptance, as elaborated in the Appendix. It is found that, with the two distinct time acceptance profiles, the fitted parameters differ by only $17.7\%$, which will be neglected in the final results.

The above simulations show that in future $Z$-factories, the proper decay time resolution can reach $4.7~{}\mathrm{fs}$, the detector
acceptance$\times$efficiency can be as good as $75\%$, and the flavor tagging power can be $17.4\%$ under a conservative assumption on the PID performance. In addition, the acceptance and efficiency is almost flat in decay time. Assuming the future $Z$-factory operating in Tera-$Z$ mode (i.e., $10^{12}$ $Z$), the scaling factor $\xi_{FE}$ is $0.0021$. The expected $\phi_{s}$ resolution is $\sigma(\phi_{s},\text{FE})=\xi_{\text{FE}}\times\sigma(\phi_{s},\text{LHCb})/\xi_{\text{LHCb}}=4.6~{}\mathrm{mrad}$.

The $\Gamma_{s}$ and $\Delta\Gamma_{s}$ depend weakly on tagging power and decay time resolution. The 3.7 times better flavor tagging power and 1.92 times better time resolution factor of CEPC, in contrast to $\phi_{s}$, have negligible effects on these observables. The estimated resolution is $2.4~{}\mathrm{ns^{-1}}$ for $\Delta\Gamma_{s}$ and $0.72~{}\mathrm{ns^{-1}}$ for $\Gamma_{s}$. The measured resolution of $\Gamma_{s}-\Gamma_{d}=0.0024~{}\mathrm{ps^{-1}}$[12] is taken as the resolution of $\Gamma_{s}$.

Figure 7 shows the expected confidential range ($68\%$ confidential level) of $\Delta\Gamma_{s}-\phi_{s}$. The black dot is the prediction of the standard model from CKMFitter group [1] and HQE theory calculation [4]. The green curves and red curves represent the expected precision of Tera-$Z$ CEPC and 10-Tera-$Z$ CEPC, respectively. The different line styles represent different PID performance assumptions, where the solid line represents the conservative PID performance assumption, which is a degradation of $30$% of the intrinsic PID performance. The blue dashed curve represents the LHCb at the HL-LHC, projected in this study. The magenta and yellow curves are the projections from ATLAS and CMS at the HL-LHC [27, 28], respectively.

The $\phi_{s}$ resolution at the 10-Tera-$Z$ CEPC can reach the current precision of SM prediction. All the future experiment measurements of $\Delta\Gamma_{s}$ can provide stringent constraints on the HEQ theory. The CEPC could do a better job on the measurement of the $\phi_{s}$ than the measurement of the $\Delta\Gamma_{s}$ because the flavor tagging and decay time resolution are excellent.

If the penguin diagram is considered in the $B_{s}$ decay, the relation between $\phi_{s}$ and $\beta_{s}$ should be corrected as

$$\phi_{s}=-2\beta_{s}+\Delta\phi_{s}(a,\theta).$$

(7)

The shift $\Delta\phi_{s}$ could be expressed as

$$\tan(\Delta\phi_{s})=\frac{2\epsilon a\cos\theta\sin\gamma+\epsilon^{2}a^{2}\sin(2\gamma)}{1+2\epsilon a\cos\theta\cos\gamma+\epsilon^{2}a^{2}\cos(2\gamma)},$$

(8)

where $a$ and $\theta$ are penguin parameters, $\epsilon=\lambda^{2}/(1-\lambda^{2})$ is defined through a Wolfenstein parameter $\lambda$, and $\gamma$ is the angle $\gamma$ of the Unitarity Triangle.

Control channels, such as $B\rightarrow J/\Psi\rho$ and $B_{s}\rightarrow J/\Psi K^{*}$, were employed to determine the penguin parameters $a$ and $\theta$, proposed by Ref. [29]. In this study, the LHCb measurements involving $B_{s}\rightarrow J/\Psi K^{*}$ are utilized to estimate the expected precision of $\Delta\phi_{s}$[30, 31]. This is based on the assumption that the findings from the $B_{s}\rightarrow J/\Psi\phi$ measurements can be directly applied to the $B_{s}\rightarrow J/\Psi K^{*}$ analysis, despite the topological differences between the two decay channels.

The observables in $B_{s}\rightarrow J/\Psi K^{*}$ measurements are

$$A^{\mathrm{CP}}=-\frac{2a\sin{\theta}\sin{\gamma}}{1-2a\cos{\theta}\cos{\gamma}+a^{2}},$$

(9)

and

$$H=\frac{1-2a\cos{\theta}\cos{\gamma}+a^{2}}{1+2\epsilon a\cos{\theta}\cos{\gamma}+\epsilon^{2}a^{2}},$$

(10)

where $A^{\mathrm{CP}}$ is the $CP$ asymmetry and $H$ is an observable constructed containing the branching fraction information, assuming the $\mathrm{SU}(3)$ symmetry.

The parameters $A^{\mathrm{CP}}$ and $H$ are polarization-dependent. The transverse components are measured at LHCb as

$$A^{\mathrm{CP}}_{\perp}=-0.049\pm 0.096,H_{\perp}=1.46\pm 0.14.$$

The constraints on the penguin parameters $a$ and $\theta$, as defined by Eq. 9 and Eq. 10, and within the ranges of $A^{\mathrm{CP}}_{\perp}=-0.049\pm 0.096$ and $H_{\perp}=1.46\pm 0.14$, are shown in Fig. 8 as blue and yellow bands, respectively. At future Tera-$Z$ $Z$-factory, the $H$ is expected to improve according to the Eq. (4), while the $A^{\mathrm{CP}}$ improves according to the Eq. (2). Therefore, the $A^{\mathrm{CP}}$ and $H$ are expected to be measured at CEPC with the precision

$$\sigma(A^{\mathrm{CP}})=0.0090,\sigma(H)=0.035,$$

taking into account the expected improvement of LHCb Run 2 compared with LHCb Run 1. The anticipated constraints on the penguin parameters $a$ and $\theta$ at CEPC are represented by green and red bands in Fig. 8, with the central values of $A^{\mathrm{CP}}$ and $H$ from LHCb measurements. The expected uncertainty of $a$ and $\theta$ is obtained by a $\chi^{2}$ fit to Eq. (9) and (10), resulting in

$$a=0.436\pm 0.023,\theta=3.057\pm 0.016^{\circ}.$$

The black contour in Fig. 8 outlines the region that falls within one standard deviation of the fitted value.

With an error propagation neglecting the correlation between $a$ and $\theta$, the precision of the penguin shift is estimated as $\sigma(\Delta\phi_{s})=2.4~{}\mathrm{mrad}$.

However, the $\mathrm{SU}(3)$ symmetry does not always hold, and thus controlling $\sigma(H)$ requires additional theoretical efforts. The degradation of $\sigma(\Delta\phi_{s})$ alongside the degradation of $\sigma(H)$ is shown in Fig. 9. To obtain $\sigma(\Delta\phi_{s})$, the expected resolution of $\sigma(A^{\mathrm{CP}})=0.0090$ at CEPC is used. The procedure for determining $\sigma(\Delta\phi_{s})$ with different $\sigma(H)$ follows the same methodology as the aforementioned statements. The rightmost point on the Fig. 9 corresponds to $\sigma(H)=0.28$, which reflects the theoretical uncertainty from the calculations in Ref. [32]. It is demonstrated that $\sigma(\Delta\phi_{s})$ is roughly linearly dependent on $\sigma(H)$, and clearly, without improved theoretical input, the control of penguin contamination will be far from satisfactory.