Unveiling Time-Varying Signals of Ultralight Bosonic Dark Matter at Collider and Beam Dump Experiments

The ultralight boson represents a promising dark matter candidate exhibiting unique wave-like behaviors. These properties could transfer to the dark mediator, such as the kinetic mixing dark photon, which can be a link between the dark and Standard Model sectors, resulting in periodic oscillations of its mass. We propose a method to detect ultralight dark matter using dark mediators in collider and beam dump experiments, distinguishing it from conventional atomic, molecular, and optical methods. The time-varying nature of dark mediator mass exhibits a double-peak spectrum, reducing traditional constraints by 1 to 2 orders of magnitude, due to decreased luminosity exposure in each resonant mass bin. To enhance sensitivity, we utilize event time-stamps in the CMS Open Data and demonstrate that this technique boosts sensitivity by approximately one order of magnitude compared to the time-blind method. Moreover, it proves effective in detecting the invisible decay of the dark mediator.


INTRODUCTION
Ultralight bosons are a promising candidate for dark matter (DM) due to their wave-like behavior, which could alter DM activity at small scales.As a result, they have become a popular area of research [1].Recent proposals have suggested that ultralight DM could produce time-varying fundamental constants [2][3][4].These changes can be detected through atomic, molecular, and optical physics, the Oklo phenomenon, and astrophysical experiments [5][6][7].Ultralight bosonic DM may also cause oscillations in Standard Model (SM) gauge couplings and fermion masses, where the oscillation period is related to the DM mass [8][9][10], as well as temporal changes due to topological defects [11][12][13].
The kinetic mixing dark photon model [14][15][16] is an illustrative example that mediates the SM and the dark sector with a mixing strength ϵ [2,[17][18][19][20].Particle experiments have placed strict limits on both ϵ and mass m A ′ [21,22].In particular, the parameter space that explains the recent muon (g − 2) µ excess [23,24] has been ruled out [22,25].However, if an ultralight scalar DM ϕ is charged under the dark U (1) ′ , it can induce periodic oscillations of the A ′ mass, leading to a significant reduction in the existing collider and beam-dump bounds.This is especially true for the dilepton resonance searches [26][27][28][29][30][31][32][33], where even the previously highly excluded dark photon solution to muon (g − 2) µ can become viable again.Furthermore, we propose the use of event-by-event time stamps of recorded events and demonstrate with CMS Open Data, which could improve the sensitivity on the signal.
In this work, we investigate the implications of the time oscillation of the mediator mass in the dark sector on high energy colliders and beam dump experiments, instead of its direct coupling with the SM sector.This phenomenon has the potential to significantly alter the experimental constraints, as it generally reduces the luminosity exposure in each resonant mass bin.Additionally, the invariant mass spectrum exhibits a multi-peak (typically double-peak) feature, which is distinct from the traditional single-peak resonance.

Time-varying mass of the particle
We assume the resonant particle has a time-dependent mass m res (t) due to environmental effects and further take a time oscillating form with a period of τ , For the resonant searches, the invariant mass of the event changes with time.Thus, the usual strategy of looking for resonance in a fixed bin suffers from reduced time exposure and leakage into other bins.If the data takes time t exp ≫ τ and the experiment analyzes the full data in a time-blind way, then the relevant physical quantity is the time exposure ∆t i in the ith mass bin [m i , m i+1 ], with the expression Instead of a narrow resonance, the signal has a spread template fully determined by dt/dm res .Then, the event arXiv:2206.14221v3[hep-ph] 28 Aug 2023 number in ith bin is where σ (i) res and ϵ i are the resonant production crosssection and cut-efficiency for the ith bin, respectively, while L is the integrated luminosity.Since particle production and decay happen very quickly, a single event's resonant mass is unchanged.Therefore, the difference between our analysis and the previous resonant analysis is fully described by the time exposure fraction ∆t i /t exp .
There is a double-peak feature in the time-varying mass scenario that the peaks must show up at the minimum and maximum of resonant mass.Assuming the function m res (t) is continuous and differentiable, the physical mass contains global minimum and maximum regardless of its periodic feature.Thus, it must have dm res /dt = 0 at two extreme points; therefore, the time exposure blows up accordingly.Additional local extrema can also contribute to peaks, leading to the multi-peak scenario.Other observable, if it is a function of timevarying mass, will inherit this property.

Model setup
We consider a kinetic mixing dark photon A ′ with U (1) ′ interaction, where ϵ is the kinetic mixing strength which controls the strength of A ′ coupling to electromagnetic current J em .The mass m 0 is a constant from U (1) ′ spontaneously breaking.In addition, we consider a complex scalar DM ϕ with small charge Q ϕ under U (1) ′ .The ultralight scalar DM obtains its relic abundance through misalignment mechanism [34][35][36][37], satisfying the equation of motion and at the late time, it is locally described by the classical wave function ϕ(t), where ϕ 1,2 are the complex field strengths, satisfying ϕ in the non-relativistic limit.In scalar quantum electrodynamics (QED), one can have the following four-point vertex It effectively leads to time oscillating A ′ mass today as where κ is the amplitude of the oscillation, Thus, the oscillation mass is fully determined by three parameters, m0 , κ, and m ϕ , with the phase removed by the definition of t = 0.
The mass ratio y(t) ≡ m A ′ (t)/ m0 has minimum y min = 1 and maximum y max = √ 1 + κ respectively, and oscillates with time period of τ ≡ π/m ϕ .From Eq. ( 2), the invariant mass bin has a time exposure proportional to where a factor of 2 is multiplied since each mass appears twice in one period.The probability density function (PDF), is normalized between y ∈ [y min , y max ].Indeed, the time exposure diverges at the minimum and maximum of the resonant mass bin in Fig. 1.After including the detector resolution, it becomes finite and shows the double-peak feature.The right peak contains a larger probability than the left because f (y) ∝ y.One can evaluate the probability difference contained in the two peaks, which is a factor of 2 difference for κ ∼ O (15).
If the data taking duration lasts much longer than the oscillation period, t exp ≫ τ , the events will run between m0 and √ 1 + κ m0 many times.In this case, the normalized mass spectrum f (y) fully describes the data distribution without explicit dependence on t, initial oscillation phase, or m ϕ .Since Lyman-alpha constraints require m ϕ ≳ 2 × 10 −20 eV [38] suggesting the longest oscillation period of about one day, most of the experiments satisfy the t exp ≫ τ condition.
Recasting via the double-peak method (DPM) The dark photon has been searched for in the dilepton channels A ′ → ℓ + ℓ − (ℓ = e, µ) with various production mechanisms [26][27][28][29][30][31].The general strategy is to fit the dilepton invariant mass spectrum m ℓℓ with a given signal hypothesis in a mass window broader than a few times of the energy resolution σ re .The upper limits on the production cross section can be translated into bounds on ϵ 2 as a function of m A ′ .To recast the BaBar and LHCb results, for a given mass window, we fit the m ℓℓ spectrum with a quadratic or cubic function and compare it to the observed data with and without the signal events to get the likelihoods L and L 0 , respectively.Then we require the log-likelihood ratio LLR ≡ −2 ln(L/L 0 ) = 3.84 to obtain the upper limit for signal event number corresponding to 95% confidence level at rejecting a signal hypothesis [40].
We first recast the traditional single-peak resonance signal method for dilepton experiments following the experimental setups summarized in Table I.Our recast results agree with the experimental results quite well, and the detailed process is in Methods and Supplemental Notes 2 and 4 [22,26,[28][29][30][41][42][43][44][45][46].We next apply our time-varying resonance signal model to fit the background data.According to the double-peak feature of the signal, for a fixed mass window centering around m A ′ , there are two signal peaks to fit, the minimum m A ′ = m 0 and the maximum m A ′ = √ 1 + κm 0 .Thus, one can obtain two sets of ϵ 2 constraints as a function of m 0 , and the stronger one is adopted as the ϵ 2 limit for a given m 0 .For a given m 0 , the best limit usually comes from the maximum.Note that one can even constrain m 0 below the dilepton mass threshold via the maximum peak; e.g., the dimuon limit of LHCb extends to m 0 much smaller than 2M µ , as shown in Fig. 2.
In Fig. 2, we apply LLR calculation to experiments in Table I.Other experiments such as APEX [47], HADES [48], KLOE [49][50][51][52], PHENIX [53] and WASA [54], generally provide relatively weaker constraints on ϵ 2 in interested parameter space.Hence, we recast their results by rescaling the limits on ϵ 2 according to the time exposure in a bin.We checked the robustness of such a simplified method by applying it to the LHCb experiment and obtaining a good agreement with the full-fitting LLR method, whose details are discussed in Supplemental Note 1.
The constraints on our time-varying signal model from the above experiments are plotted in Fig. 2. For m 0 ≳ 10 −2 GeV and κ ∼ O (10), current experiments constrain ϵ 2 ≳ 10 −7 − 10 −5 , around one order weaker than the traditional single-peak bounds, whose envelope is shown as gray dashed line labeled as κ = 0. Especially, the excluded (g − 2) µ solution becomes viable at O(10) MeV when the (g − 2) µ red band crosses with the BaBar NA48/2 bounds, as depicted in Fig. 3.
There are beam dump experiments E774 [41], E141 [42] and NA64 [43].They set limits on A ′ → ℓ + ℓ − based on the signal event number N (ϵ, m A ′ ), for given ϵ and m A ′ .The A ′ is produced at beam dump and propagates a distance according to its lifetime and decay to ℓ + ℓ − .We can translate the upper limit on the event number of A ′ decay, e.g., 17 events for E774, to our scenario by simply time averaging the signal events as then compare it with the upper limit as shown in Fig. 2.
The detailed estimation of N (ϵ, m A ′ ) is given in Methods.In addition, the LHCb limits for the displaced A ′ are shown in Fig. 2 with details given in Supplemental Note 2. For A ′ invisible decaying [55][56][57][58], the doublepeak method can apply in the same way.Improving by the time-dependent method (TDM) Previous calculations only consider two parameters, κ and m 0 , while do not exploit the recorded time-stamp of events, t.The experiments can reanalyze the data using all the above three parameters, because the signal events only happen at certain time t and mass m A ′ (t), as shown in Fig. 4. In principle, the experiments can figure out κ, m 0 , m ϕ and the initial phase via a two-dimension fit on the t-m ℓℓ plane; for example, dark matter mass m ϕ can be extracted by the the period of mass modulation of dark photon via τ = π/m ϕ , while the dark photon bare mass m 0 and κ can be extracted via the minimum and maximum of signal m ℓℓ distribution, etc.Note that a direct probe of the ultralight dark matter mass m ϕ is possible only in TDM through the time-varying feature of the dark mediator A ′ .Without the time information, we can assume the observed data has a flat probability in time and estimate the signal sensitivity.We adopt the same mass grid as the experiment, and it automatically generates the time grid according to the signal mass function m A ′ (t).Specifically, if the total number in the ith mass bin is N i , we have the number of data in ith mass bin and jth time bin as N i,j = N i ∆t j /τ .For a fixed mass bin (horizontal gray shaded), we pick up the two red bins in Fig. 4, which contain the signal.Adding the data in red bins together, we have forming an updated set of data N red i .Then, our previous calculation using the DPM can apply.
In this method, the signal event number is unaffected while the background event is suppressed by a factor of 1 τ mi+1 mi dt dm A ′ dm A ′ .Fig. 5 shows the projected sensitivities using TDM for NA48/2, BaBar, and LHCb as examples, assuming no excess is detected.Indeed it improves the reach by 1-2 orders of magnitude compared with the DPM, and the viable regions of (g−2) µ at O (10) MeV in Fig. 2 can be probed now.We then propose the experimentalists to reanalyze the data using TDM to probe the time-varying signal model.Reducing the invariant mass resolution and hence the size of m ℓℓ bin can reduce the number of background events in the corresponding mass bin, which enhances the sensitivity to signal, as the number of signal events is unaffected when picking up the red bins along the oscillation trajectory in Fig. 4.However, if the mass bin is too small, the analysis will suffer from statistical error due to small N red i .

The invisible dark photon
Due to small g ′ Q ϕ , the dominant decay channels of A ′ are SM fermions in our minimal setup.However, in an extended dark sector, A ′ may decay to invisible particles dominantly.For example, if some dark fermion χ that carries dark U (1) ′ charge, then A ′ → χ χ channel can dominate the A ′ decay with the subsequent decay of χ to DM.In that case we have an invisible decaying dark photon scenario, and it has been searched by several experiments including BaBar [55], BES-III [56], NA64 [57], and NA62 [58].We will briefly discuss how the time-varying scenario can affect the results.BaBar and BES-III have studied the monophoton channel, e + e − → A ′ γ with invisible A ′ .The photon energy is s being total collision energy.With time-varying m A ′ (t), E γ extends to a spectrum determined by the following differential, where Then, the analysis is similar to visible A ′ by substituting the invariant mass bin for the photon energy bin, and it is expected to weaken the limit.The exception happens for very small m 0 satisfying m 2 0 < 2 √ sσ γ /κ, with σ γ being the photon energy resolution.In this case, the limit will be unchanged because both E min and E max fall into the same photon energy bin.
The electron beam dump experiment NA64 [57] studies the process e − Z → e − ZA ′ for invisible A ′ .The signal events are selected with E miss > 50 GeV and other cuts on electromagnetic and hadronic calorimeter energies.Since E miss is much larger than interested m A ′ , the cut efficiency should not significantly depend on m A ′ .Therefore, the dominant effect of time-varying m A ′ will show up in the production rate, scaling as m 2 e /m 2 A ′ [45].Thus one can take the time average for this factor to estimate the weakening of the limits.Another beam dump experiment NA62 [58] focuses on invisible A ′ from π 0 → γA ′ , where A ′ mass is reconstructed using the process K ± → π ± π 0 as m 2 res = (p K ± − p π ± − p γ ) 2 .The double-peak method can analyze the data and set updated limits.
The results are shown in Fig. 6, in which the (g − 2) µ parameter space has been entirely excluded, even for the time-varying signal.Therefore, we shall assume the invisible decay channel is subdominant to have a dark photon (g − 2) µ solution.

Analyzing CMS Open Data
The 2012 dimuon events from CMS Open Data [59,60] can be used to justify the two analyzing methods, DPM and TDM.We recast the CMS analysis with the traditional single-peak model shown as κ = 0, in Fig. 7.
For the real data, there is a non-trivial feature: the instant luminosity L(t) is non-uniform.It introduces time dependence from human operation aside from the time oscillation from the signal.Therefore, the invariant mass distribution becomes where ϵ S and σ 0 are the cut efficiency and the production cross-section at m A ′ = m ℓℓ , and t ± i are the two solutions in the i-th time period satisfying m A ′ (t) = m ℓℓ .When the oscillation period τ is smaller than the typical L(t) variation time (about a few hours), L(t) becomes a constant in each period.Therefore, τ 2 i L(t + i ) + L(t − i ) = L and we can apply DPM to CMS Open Data.In addition, we can also apply TDM to CMS Open Data after filling the data into t-m ℓℓ grids.More details for the analysis of CMS Open Data are given in Methods [61][62][63][64][65].
The results from DPM and TDM are plotted in Fig. 7 for κ = 15, 24, ϕ 0 = 0, and m ϕ = 10 −19 eV.As expected, the limits from DPM are weaker than our recast for the traditional limit (κ = 0).Moreover, the limits for TDM are indeed better than DPM by 1-2 orders.We show that TDM can work for the real collider data.

Other constraints
Besides collider and beam dump experiments, there are other constraints to clarify.Firstly, the coupling g ′ Q ϕ should be small to avoid ϕ thermalization via selfscattering and scattering with normal matter via f ϕ → f ϕ or f f ↔ ϕϕ * .Secondly, the A ′ → ϕϕ * decay could freeze-in ϕ as a hot relic, which should be very small.For interested parameter space m 0 ∼ O(0.1) GeV, the coupling g ′ Q ϕ around 10 −6 -10 −10 , is small enough to satisfy the above requirements for m ϕ ∼ 10 −20 -10 −17 eV, while keeping κ ∼ O (10).Moreover, in the early Universe, the field value of ϕ is much larger than today.Thus, a heavy A ′ mass helps to evade the thermalization and freezein constraints.For the ultralight scalar, the black hole super-radiance can exclude some mass regions, but not all the interested regions [66,67].
At 1-loop level, the SM fermion mass can receive a QED-like correction from A ′ interaction, with m A ′ ≫ m f .This leads to a logarithmic coupling between (ϕ * ϕ) and the fermion mass operator.Its Taylor expansion can not be naively truncated due to large κ, hence is different from the linear and quadratic couplings.Especially for experiments relying on certain ϕ field distribution around massive objects, such as Cassini stochastic, binary pulsars tests [68,69], atomic clocks [70,71], torsion balances [72][73][74], and MICROSCOPE space experiment [75], the constraints do not simply apply.As for the fifth force experiments [76,77], even for the quadratic coupling, it only provides a loose bound [10,78], which can be easily satisfied for small ϵ and g ′ Q ϕ .
The constraints from Big Bang nucleonsynthesis due to the enhanced ϕ value [10,79,80], can be easily evaded in our scenario thanks to the logarithmic coupling.

CONCLUSION
We have introduced an innovative strategy to search for ultralight bosonic DM, namely the time-varying resonance from the dark sector at the collider and beam dump experiments.We found it can lead to double-peak feature in the invariant mass spectrum, which can help to evade the dilepton and missing mass resonant searches at collider and beam dump experiments.Moreover, the mass spectrum is independent of time as long as the oscillation period is short compared with the variation time scale of instant luminosity.A concrete model is discussed with ultralight complex scalar DM inducing an oscillating mass for kinetic mixing dark photon.For mass around tens of MeV, the already excluded muon (g − 2) µ solution from A ′ becomes viable again; this parameter region can be further tested by reanalyzing the existing data with event-by-event time stamps, as shown by the time-dependent method.We use CMS Open Data for a time-dependent resonance search and justify that our method works as expected, even with the complexity of a non-uniform instant luminosity.We also demonstrate its application in the case of invisible decay dark mediator at the electron-positron collider.In general, the analysis strategy can search the time-varying signal from the ultralight DM oscillation at collider and beam dump experiments.

METHODS
We show the detailed log-likelihood ratio (LLR) calculations for the collider and beam dump experiments.Next, the recasts for the existing experiments are in good agreement with the official results.Therefore our calculations using double-peak (DPM) and time-dependent methods (TDM) are robust.Aside from providing projected limits for existing experiments, we use CMS Open Data to provide real collider limits using DPM and TDM.Lastly, we constrain time-varying A ′ , which decays to the invisible final state using the DPM.

The LLR calculation for collider experiments
In this section, we will show the details of LLR calculation which sets limits on dark photon signal cross section or mixing parameter ϵ.Taking the BaBar experiment as an example, we recast the limits of the traditional single-peak resonance model (κ = 0) and compare them with official results.Good agreements are obtained.We then adopt the time-varying signal model and derive the updated bounds for both the double-peak and timedependent analysis strategies.The details of the recast and updated bounds of LHCb, and A1 experiments are given in Supplemental Notes 2, 3 and 4, including the long-lived dark photon at the LHCb experiment for the time-varying scenario.
The BaBar collaboration collected 514 fb −1 data at the vicinity of the Υ(4S), Υ(3S) and Υ(2S) resonances to search for e + e − → γA ′ process, with A ′ → e + e − and µ + µ − decay channels within a mass range 0.02 GeV < m A ′ < 10.2 GeV [26].A total of N e = 5704 (N µ = 5370) mass hypotheses are searched in the e + e − (µ + µ − ) channel respectively.For a given mass m A ′ , an interval of m A ′ ± 10 σ re is used to perform the fits, where σ re is the energy resolution at m A ′ , varying from 1.5 to 8 MeV in the whole m A ′ range.Even though the full data is not given in Ref. [26], there are available data of m ee and m R = m 2 µµ − 4M 2 µ spectrum up to ∼ 10 GeV with a uniform bin size of 100 MeV in the main text, and the zoomed-in spectrum for m ee ∈ [17.8, 62.2] MeV (m R ∈ [0.522, 77.4] MeV) with a bin size of 0.5 MeV (1.0 MeV) in the appendix.We refer to the former and latter as the low-and high-granularity datasets, respectively.
We first generate the artificial high-granularity data for the whole mass spectrum to recast the BaBar results.For m ee < 62.2 MeV (m R < 77.4 MeV), we adopt the high-granularity data themselves, while for higher mass, we use the interpolation of the low-granularity data, including appropriate statistical smearing.We assume the bin size increases linearly with resonant mass and keep the total number of the data to be N e (N µ ), respectively.Moreover, we assume the resolution σ re increases linearly with the resonant mass.Finally, for a given mass point m A ′ , we fit the artificial data in the mass window m A ′ ± 10 σ re .The LLR is defined as where B i is the background event number in the ith mass bin, and is the background fitting function with m i being the center value of ith bin for m ee or m R , while is the normalized Gaussian distribution and S is the total signal event number.f G (m i ) is the signal template without the time-varying effect, and after detector smearing, it is defined as In the LLR calculation, only statistical error is considered, and there is an extra Jacobi factor for the dimuon channel from the definition of m R .
After requiring LLR = 3.84 (the 95% confidence level at rejecting a signal hypothesis), we obtain the limits on the allowed signal total event S. They can be used to unfold the limit on σ(e + e − → γA ′ ) via the acceptance factor 0.15 (0.35) in the dielectron (dimuon) channel respectively [26].We found our recast results are consistent with BaBar's official results, as shown in Fig. 8(a) for both e + e − and µ + µ − channels.To test our recast result, we define a ratio R = σ Recast S /σ Official S , which is shown as the low panels of Fig. 8(a).In Fig. 8(b), we show the distribution of log 10 R for our simulated points.As expected, the log 10 R is distributed around 0. We fit the distribution with a standard Gaussian function , where N tot is the total number of points, µ and σ is the mean value and standard deviation respectively.For e + e − and µ + µ − channels, we found similar Gaussian shapes with σ = 0.28 and µ ≈ 0. Therefore, our recasts are quite close to BaBar's official results.Our deviation from the official result is within a factor of 10 0.28 ≈ 1.9 at 1σ level.As a result, our recasts are quite robust.Therefore, our LLR calculation and the projected limits can be trusted.
Regarding the time-varying scenario, the signal invariant mass spectrum follows the probability density function as defined in Eq. ( 12) of the main text, with y = m ℓℓ /m 0 and y takes value from 1, √ 1 + κ .With Gaussian smearing, the signal template becomes with m min = m 0 and m max = √ 1 + κm 0 .In the µ + µ − channel, the additional Jacobi factor should be considered as m i refers to m R , not m µµ .The time-varying scenario is then fitted using f S in Eq. (20).As f S peaks

BaBar
Our Recast around m min and m max , for a given m A ′ we perform two independent fits for m 0 = m A ′ and m 0 = m A ′ / √ 1 + κ respectively.Therefore, for a given m A ′ , we can obtain two sets of limits on the allowed S corresponding to the left and right peaks, respectively.
To reduce the systematic uncertainties, we take the ratio of S from the time-varying resonant mass scenario f S and the traditional resonant scenario f G to rescale the ϵ 2 constraints from the BaBar measurement [26].Since the right peak of f S is usually higher than the left peak, for a given m 0 , the most stringent constraint usually comes from the mass window around m A ′ = √ 1 + κm 0 .Therefore, this is adopted as the BaBar constraints on our time-varying resonance scenario.

The detailed limits calculation for beam dump experiments
The beam dump experiments E774 [41], E141 [42] and NA64 [43], are all electron fixed-target experiments.The dark photon A ′ are dominantly produced by electron Bremsstrahlung process, e − Z → e − ZA ′ [44,45].The A ′ will travel some distance and decay before reaching the detector.Since there is shielding behind the collision target, A ′ must have a lifetime with macroscopic scale.To simplify our analysis, we employ the estimate of the A ′ signal event number following Refs.[22,[44][45][46], A ′ e −a1L sh Γ A ′ (1 − e −a2L dec Γ A ′ ), (25) where N e is the total electron number in experiments, C ′ is a parameter defined in Ref. [45] with typical value of 10, a 1 and a 2 are fitting parameters, Γ A ′ is the decay width of A ′ , L sh is the distance of the end of the shield and L dec is the distance of detector from the collision point.The constraints are obtained by requiring N (ϵ, m A ′ ) equal to the allowed signal events, for example, 17 events for E774.We adjust the fitting parameter a 1 and a 2 to reproduce the original results from the experiments.The results are shown in Fig. 9, and one can see the fits are quite well.
Then, with the obtained a 1,2 we constrain the timevarying scenario by replacing N (ϵ, m A ′ ) with its time average, namely Finally, setting N (ϵ, m 0 , κ) to the allowed signal events, we obtain the limits for the time-varying scenario.

The Analysis of the CMS Open Data
We test the time-varying scenario using real collider data.Unfortunately, the BaBar data are unavailable, while the LHCb and ATLAS Open Data are inadequate for scientific study.Only the CMS Open Data provides the full complexity of the collision data and is adopted in our analysis.We analyze the CMS dimuon sample in the AOD format, collected in 2012 with a collision energy of √ s = 8 TeV [59,60].The CMS has recorded an integrated luminosity of 20.6 fb −1 for 8 TeV in 2012 [64], but only provide 50% to the public in the form of CMS Open Data [65] according to the CMS data policy [62].
For a time-varying signal, the greatest challenge is that the actual experimental data are not taken continuously and uniformly.The data delivering follows the scheduled program, which leads to a non-uniform instant luminosity L(t), as shown in the left panel of Fig. 10.In the right panel of Fig. 10, we plot the number of dimuon events in a small time interval (10 seconds) and divide it by the integrated luminosity in this interval.As expected, the event rate per unit luminosity is roughly constant since SM cross-section is constant with time.The instant luminosity induces time dependence from human operation aside from the intrinsic time-oscillation of the signal.Therefore, we must improve our previous analysis, which assumed a uniform and continuous instant luminosity, and adapt improved strategies to the experiment data.
For the double-peak method, we perform a onedimensional analysis on the dimuon invariant mass distribution m ℓℓ with ℓ = µ.In general, the distribution should be where ϵ S (m ℓℓ ) and σ 0 (m ℓℓ ) are the cut efficiency and the production cross section of pp → A ′ (→ ℓ + ℓ − )X with m A ′ = m ℓℓ , δ(...) is the Dirac δ function, m A ′ (t) = m 0 1 + κ cos 2 (m ϕ t + ϕ 0 ) is the time-dependent reso-nant mass varying from m 0 to √ 1 + κm 0 , and the oscillation period is τ = π/m ϕ .The δ function can be rewritten as where y = m ℓℓ /m 0 , f (y) is the probability density function for the invariant mass spectrum, t ± i are the two solutions for m A ′ (t) = m ℓℓ for the ith oscillation period, with i ∈ Z. Their explicit expressions are The expected event number is the time integral of the product of cross section dσ S /dm ℓℓ and instant luminosity L(t).Substituting Eq. ( 28) into Eq.( 27), one gets the signal event number in an invariant mass bin during a time interval [t 1 , t 2 ] where i is the ith time period within [t 1 , t 2 ].If the integrated luminosity is denoted as L, we can define the probability density function for the real data f (y) as When the oscillation period is smaller than the typical L(t) variation time (about a few hours), the instant luminosity in each period is roughly constant, e.g., L(t) = L 0 i in the ith oscillation period.Therefore, for τ ≲ a few hours, we have In Fig. 11, we use the instant luminosity L(t) from CMS Open Data to calculate the distribution f (y) (blue)   and compare it with the ideal signal distribution f (y) (yellow).We choose κ = 15, the initial phase ϕ 0 = 0 and three different m ϕ = 10 −20 , 10 −19 , 10 −18 eV as benchmark points.For the blue lines of f (y), we see significant fluctuations and deviations from the yellow line f (y), for m ϕ = 10 −20 eV, which originates from the non-uniform instant luminosity of CMS Open Data.For m ϕ = 10 −19 eV (τ = 5.75 hours), the oscillation time scale is short compared with the luminosity variation time scale.Hence, the real signal spectrum f (y) is very close to the ideal spectrum f (y).For even shorter oscillation period, m ϕ = 10 −18 eV, f (y) is almost the same as f (y).Therefore, we demonstrated that for oscillation period τ smaller than a few hours, the signal invariant mass distribution is not affected by the non-uniform instant luminosity.This conclusion remains unchanged when varying the initial phase ϕ 0 .
As a result, for a small time oscillation period, e.g., τ ≲ a few hours, the double-peak method described in the main text can apply to CMS Open Data without any further modification.We select the 2012 dimuon sample in the CMS Open Data set and apply minimal cuts as p µ T > 15 GeV, |η| < 2.4.We aim to compare the limits for the traditional scenario κ = 0 and the time-varying scenario using double-peak and time-dependent methods.In the recast for the traditional scenario , our procedure and limit using CMS Open Data are similar to the recast of Ref. [63] using CMS data in Ref. [64].
For the DPM, we smear the signal distribution with an energy resolution σ re = 0.026 GeV + 0.013 m ℓℓ [61], and take bins of m ℓℓ with a width of ∆m = σ re .For each m A ′ , we fit the signal model for m 0 = m A ′ and m 0 = m A ′ / √ For the time-dependent method, we read the invariant mass m ℓℓ and time-stamp t of each event and then fulfill the events into a two dimensional t-m ℓℓ bins.We use the letter i to denote the ith time bin and j to denote the jth invariant mass bin.In the 2D grid, we only take the bins along the m A ′ (t) trajectory, as illustrated in Fig. 3 of the main text, which roughly satisfies m j ≈ m A ′ (t i ).We denote these bins along the signal trajectory as (i, j) pairing, where j is determined by i.Therefore, the expected signal event within such (i, j) bin is where the first factor is the Gaussian smearing of the m ℓℓ resolution, the second term is the instant luminosity, and the last term is signal cut efficiency times the production cross-section.For the two-dimensional binning, we use ∆m = σ re and ∆t = τ /8.We have tried more fined grids for t, but the result has not improved.Because the current grid is already small enough, the observed event number in each grid is quite small, or even zero.Therefore, further refining the grid will suffer more from statistical error.Given the signal and background event numbers in each bin, we use Eq. ( 20) to derive the LLR by summing up all the (i, j) bins.Finally, we plot the limit from the TDM in Fig. 5 in the main text for κ = 15, 24, ϕ 0 = 0, and m ϕ = 10 −19 eV.The limit from TDM is indeed better than DPM by 1 − 2 orders of magnitude.
In summary, we have successfully placed the limits using the actual collider data for a time-oscillating signal, which opens up the time-oscillating analysis to the exotic category searches at the collider.It is a pathfinder search for this broad class of time-varying signals, not limited to changing mass but also couplings, production cross-sections, and decay rates.Furthermore, we demonstrated the feasibility and the potential of a renewed class of exotic signals, namely the time oscillation signal, which is useful for experimental analysis.

DATA AVAILABILITY
The data extraction code of CMS Open Data is included in https://github.com/gitguojh/cms.open.data.code.git,which is based on the example code of the dimuon spectrum from a CMS 2011 primary dataset [82].By analyzing the AOD file of CMS Open Data, aside from the basic information of dimuons, such as the invariant mass m µµ , transverse momentum p T , and pseudorapidity η, we also extract the the time-stamp and instant luminosity information of each event using the structure timespec in Level1TriggerScalers.h and the function instantLumi in LumiScalers.h.With this event-by-event information in hand, both the DPM and TDM can be applied.Additional details can be accessed in README.md,demoanalyzer cfg.py, and scr/DimuonSpectrum2011.cc at the aforementioned Github link.And all relevant data are available from the corresponding author upon request.limits on the time-varying scenario can be similarly obtained by setting nA ′ ob [m 0 , κ, ϵ 2 ] ≤ nA ′ ex [m 0 , κ, ϵ 2 ].In Fig. S3, we plot the excluded regions for κ = 15, 24 of the time-varying scenario in the cyan-shaded region and compare it with the traditional resonance model (κ = 0) in gray dashed contours.First, the exclusion regions are much smaller than κ = 0. Second, the exclusion region shift to the left by a factor around √ 1 + κ.These two features can be understood from the double-peak feature of the signal.The figure shows the sensitivity region of LHCb in the gray dashed contours for traditional resonance scenario, around m ′ A = m µ + µ − ∼ 250 MeV and ϵ 2 ∼ 10 −9 .A time-varying signal has a double-peak feature in the invariant mass spectrum, one is m 0 = m ′ A , and the other is √ 1 + κm 0 = m ′ A .The right peak of m A ′ has a probability larger than the left peak by ∼ (1 + κ) 1/4 , as shown in Eq. ( 13) in the main text.Therefore, the exclusion region of the time-varying scenario shrinks and shifts to the left around m 0 ∼ 250/ √ 1 + κ MeV.
SUPPLEMENTARY NOTE 4: A ′ → e + e − DECAY AT A1 An analogous LLR calculation is performed to A1 experiment in the fixed-target electron scattering process e − Z → e − ZA ′ , with prompt decay A ′ → e + e − for 0.040 GeV ≲ m A ′ ≲ 0.300 GeV in 2014 [30].With the data taken from Ref. [30], we have recasted the exclusion limits with Gaussian smearing.In the left panel of Fig. S2, one can see the recast is very close to the A1 official result for the fixed resonant scenario.The A1 experiment [30] provides m e + e − data down to 30 MeV.However, the number of events there is quite small and thus will suffer from large statistical errors.Therefore, their official analysis stops at 40 MeV.In addition, we plot the ratio R between our recast and the A1 official result and the distribution of log 10 R in Fig. S2.The results show that our recast agrees with the official results well.The deviation is tiny, within a factor of 10 0.16 = 1.45 at 1σ level.Like BaBar and LHCb experiment, the double-peak method is again applied to obtain new limits on ϵ 2 as a function of m 0 .

y 15 fFigure 1 .
Figure1.The normalized probability density function (PDF).The blue solid and red dashed lines represent the PDF before and after smearing with a detector resolution of 3%, respectively.

Figure 2 .
Figure2.Limits on mixing strength ϵ 2 as a function of mass parameter m0.The colored shaded regions in (a) and (b) represent the excluded parameter space using the doublepeak method for κ = 15 and 24, respectively, and the gray dashed line represents the traditional limit (κ = 0).The redshaded region can provide a solution to (g − 2) µ .

Figure 3 .
Figure3.Limits on mixing strength ϵ 2 in the mass range m0 ∈ [0.01, 0.07] GeV.The colored shaded regions in (a) and (b) represent the excluded parameter space using the double-peak method for κ = 15 and 24, respectively, and the gray dashed line represents the traditional limit (κ = 0).The red-shaded region can provide a solution to (g − 2) µ .

2 κ=15Figure 4 .
Figure 4. Time grids used in the time-dependent method.The resonant mass curve is represented by the black solid line.For a fixed row (gray), two particular bins (red) are chosen in the time-dependent method.
DPM: double-peak method TDM: time-dependent method

2 CMSFigure 7 .
Figure 7. Limits on mixing strength ϵ 2 from CMS Open Data.The initial phase is ϕ0 = 0 and the DM mass is m ϕ = 10 −19 eV, while κ = 15 (a) and 24 (b).We show the recast of traditional limit (κ = 0) in gray dashed, and DPM and TDM limits in purple and blue.

Figure 8 .
Figure 8. Recast of the BaBar data.(a): The upper limit on signal cross-section without the time-varying effect, as a function of m A ′ in dielectron (e + e − ) and dimuon (µ + µ − ) channels for the BaBar experiment [26].The existing and our recast limits are plotted as solid and dashed.(b): the function R shows the ratio between our recast and official results.(c)(d): The distribution of log 10 R for e + e − and µ + µ − channels.

Figure 9 .
Figure 9. Recast of the beam dump experiments.The constraints obtained from our estimation Eq. (25) are shown in dashed lines, comparing with corresponding results obtained from experiments and early analysis shown in solid lines.

Figure 10 .
Figure 10.Luminosity records from the CMS Open Data.(a): the instant luminosity L(t) recorded by the CMS dimuon data 2012 [59, 60], with time t = 0 starting from 2012.07.07, 23:14:58.(b): the ratio between the event rate and instant luminosity is approximately constant.

Figure 11 .
Figure 11.Realistic signal normalized spectrum based on the CMS Open Data.The realistic spectrum f (y) is shown in blue for y = m ℓℓ /m0 based on L(t) from CMS Open Data, compared with the ideal spectrum f (y) (shown in orange) assuming uniform instant luminosity.Three different oscillation periods are shown, with m ϕ = 10 −20 eV (a), 10 −19 eV (b) and 10 −18 eV (c).With small τ , f (y) becomes close to f (y).

Table I .
The summary table for experiments using the dilepton resonance to search for A ′ .