$Z$-boson decays into an invisible dark photon at the LHC, HL-LHC and future lepton colliders

We study the decay of the $Z$ vector boson into a photon and a massless (invisible) dark photon in high-energy collisions. The photon can be used as trigger for the event, while the dark photon is detected indirectly as missing momentum in the event final state. We investigate the possibility of searching for such a dark photon at the LHC, HL-LHC and future lepton colliders, and compare the respective sensitivities. As expected, the best result is found for the lepton colliders running at the $Z$ mass, FCC-ee and CEPC, with a final sensitivity to branching ratios of order $O(10^{-11})$. We also discuss how to use the photon angular distribution of the events in lepton collisions to discriminate between the dark photon and a pseudo-scalar state like the axion.


I. MOTIVATIONS AND SYNOPSIS
The two-body Z boson decay into an isolated photon γ and a stable (or meta-stable) neutral particle beyond the Standard Model (SM), effectively coupled to the Z and γ gauge bosons, can give rise to a peculiar experimental signature at high-energy colliders.
The kinematical properties of the detected photon and of the neutral particle, which are both monochromatic in the rest frame of the decaying Z, makes this process quite attractive in the search for new physics effectsmost notably in the case of lepton colliders, for which the monochromatic photon energy is smeared only by Bremsstrahlung radiation and detector effects. This feature is lessened at hadron colliders because of the additional challenges in reconstructing the rest frame of the Z-boson, due to the characteristic uncertainties in the transverse missing momentum measurement.
What are the possible candidates for such a neutral particle?
Within supersymmetric theories, neutral states can be either fermions (neutralinos and gluinos), which, like the SM neutrinos, would only contribute to Z three-body decays, or scalars (sneutrinos), which are too heavy to be produced in the Z-boson decay.
Other theoretical frameworks remain open. Here we consider the case of the dark photon, which is a particle belonging to a dark sector-a sector comprising particles interacting only feebly with SM states, and fashioned as a generalization of dark matter (see e.g. [1] for a recent review).
A dark photonγ is the gauge boson of a U (1) D group under which dark matter as well as all other dark particles are charged. In particular, a massless dark photon (corresponding to an unbroken U (1) D gauge group) does not couple directly to the SM currents, but only through a dipole-like operator of dimension-five [2,3], a distinguishing characteristic that makes the massless case very different from the massive one. As far as its phenomenology is concerned, the massive dark photon interacts (via mixing effects) like a SM photon. Since the decay of the Z boson into two photons is forbidden by the Landau-Yang theorem [4,5], a Z boson cannot decay into a photon and a massive dark photon via mixing. On the other hand, the decay of the Z boson into a photon and a massless dark photon is mediated by a one-loop diagram of SM fermions containing theγ dark dipole operator and the usual electromagnetic vertex for the emission of the SM photon. Then the Landau-Yang theorem does not apply to this case because the interaction vertices of the dark and the ordinary photon are distinguishable. A branching ratio (BR) for Z → γγ between O(10 −9 ) and O(10 −11 ) can be expected [6].
In the present Z decay channel, the missing energy could be carried away by neutral bosons other than the dark photon. For example, an axion-like particle (ALP) has been considered and found to have a BR as large as O(10 −4 ) [7,8]. Also more exotic cases have been suggested: a Kaluza-Klein graviton in models of large extradimensions (with a BR around O(10 −11 )) [9] and an unparticle (for which a BR as large as O(10 −6 ) is possible) [10]. The signature in the latter two cases will be different since the photon is not resonant due to the continuum spectrum of the missing mass. The same signature is shared by the irreducible background Z → γνν [11].
At the experimental level, the process where X 0 stands for undetected neutral particles, was explored at the Large Electron-Positron Collider (LEP) and a limit of 10 −6 at the 95% CL was found [12] for the corresponding BR when considering a massless X 0 in the arXiv:2006.15945v1 [hep-ph] 29 Jun 2020 final state. In this paper, we investigate the possibility of searching for such a massless dark photon in Z decays at the Large Hadron collider (LHC), the High Luminosity LHC (HL-LHC), and future lepton colliders, and compare the respective sensitivities.
The paper is organized as follows. In section II, we provide the theoretical framework for the effective Zγγ couplings, and the expression for the corresponding effective Lagrangian. In section III, we study, through Monte Carlo (MC) simulations, the upper limit on the BR of the Z decay into a photon and a dark photon, that can be reached using data collected at the LHC. Performing a simple extrapolation from the result obtained for the LHC, we also estimate the limit on the BR which can be obtained at the HL-LHC. In section IV, we extend the study to Future Circular Colliders, namely the FCC-ee and the CEPC. At the electron-positron colliders, as expected, the lower level of background compared to hadron colliders and the production of large samples of Z bosons provide the most stringent limit on the BR. In section V, assuming that the decay has been observed, we look at the angular distribution of the events to establish the spin of the particle carrying away the missing energy (see [13] for a discussion of the spin dependence of the signature), and determine a lower bound on the number of events needed to distinguish between the case of the spin-1 dark photon from a spin-0 pseudo-scalar. In section VI, we give our conclusions.

II. THEORETICAL FRAMEWORK
We consider the effective coupling of the photon γ and the Z gauge boson to a massless dark photonγ. The SM fermions couple to a massless dark photonγ only through radiative corrections induced by loops of dark-sector particles. The starting point is thus given by considering the leading contribution provided by the magnetic-and electric-dipole operators, described by the effective Lagrangian where the sum runs over all SM fields ψ f , B µν and e D are the U (1) D dark photon field strength and elementary charge (we assume universal couplings), and Λ the effective scale of the dark sector. The magnetic-(d f M ) and electric-dipole (d f E ) factors can arise for instance at one-loop order, in the framework of a UV completion of the dark sector, in which there are messenger fields providing an interaction between SM and dark fields, as discussed, for instance, in [14]. In this case, the scale Λ will be associated to the characteristic mass scale of the new physics running in the loop. The scale Λ and the couplings d f M,E can be considered as free parameters.
As shown in [6], a non-vanishing effective coupling Zγγ can be generated at one loop, inducing the manifestly gauge invariant effective Lagrangian where e is the unit of electric charge, Λ is the scale of the new physics, the dimension-six operators O i are given by , darkphoton (B µ ) and photon (A µ ) fields, respectively, and F µν ≡ ε µναβ F αβ is the dual field strength. The coefficients C i are dimensionless couplings that can be computed from the UV completion of the theory. For the case of the Lagrangian in equation (2) they are a function of the couplings d f M , the U (1) D unit of charge e D , the SM fermion masses, and the Z mass [6].
Analogously, the Lagrangian induced by the electricdipole moment is where the dimension-six operator is and the expression for the coefficient C E for the Lagrangian in equation (2) can be found in [6]. The operators in equation (3) and equation (7) are CP even and odd, respectively. The amplitudes in momentum space for the decay Z → γγ can be found by taking into account the effective Lagrangians in equation (3) and equation (7); the total width is given by [6] where C M = i C i . In the following we study BR(Z → γγ) = Γ(Z → γγ)/Γ Z tot as the observable providing a direct probe to the Z → γγ process, investigating its discovery potential both at present and future hadron and lepton colliders.

III. HADRON COLLIDERS
The LHC [15] is a circular superconducting protonsynchrotron situated at the CERN laboratory, which accelerates and collides protons at a center-of-mass (CM) energy of 13 TeV. LHC hosts two general purpose experiments which study the collision products: ATLAS [16] and CMS [17]. In this paper we assume an ATLASlike detector, simulating events at 13 TeV with a total integrated luminosity of 140 pb −1 , a choice that reproduces conditions similar to those obtained during Run-2 at the LHC. The HL-LHC [18] is the foreseen upgrade of the LHC collider to achieve instantaneous luminosities a factor of five larger than the present LHC nominal value. The upper limit on BR(Z → γγ) at the HL-LHC has been derived from the one obtained for the LHC by taking into account the increase in luminosity [18].

The ATLAS detector
The ATLAS detector, used in the following as a benchmark experimental framework for simulations at hadronic colliders, is a 46 m long cylinder, with a diameter of 25 m. It consists of six different cylindrical subsystems wrapped concentrically in layers around the collision point to record the trajectories, momenta, and energies of the particles produced in the collision final states. A magnet system bends the paths of the charged particles so that their momenta can be measured. The four major components of the ATLAS detector are the Inner Detector (ID), the Calorimeter, the Muon Spectrometer and the Magnet System. The ID components, embedded in a solenoidal 2 T magnetic field cover a pseudorapidity range of |η| < 2.5. The calorimeters cover the range |η| < 4.9, using different techniques suited to the widely varying requirements of the physics processes of interest and of the radiation environment. The total thickness, including 1.3 λ from the outer support, is 11 λ at η = 0. In the muon spectrometer, over the range |η| < 1.4, magnetic bending is provided by the large barrel toroid. For 1.6 < |η| < 2.7, muon tracks are bent by two smaller end-cap magnets. Over the so-called transition region 1.4 < |η| < 1.6, magnetic deflection is provided by a combination of barrel and end-cap fields [16].
A. Methods

Monte Carlo simulation
At the LHC, we study the discovery potential for the decay Z → γγ through the signal process pp → Z → γγ. The Z-boson production fiducial cross section σ f id was calculated at the Leading Order (LO) in QCD using MadGraph5 aMC@NLO [19] with parton distributions from NNPDF23 [20], summing over all fermionic pp → Z/γ * → ff contributions and imposing on the Z decays the same fiducial cuts required on the signal single detected photon, that is, a minimum p T of 10 GeV and |η| < 2.5. The obtained Z-boson production fiducial cross-section turns out to be equal to with an expected number N Z of Z-bosons produced at the LHC at its design luminosity of 300 fb −1 of N Z = 7.5×10 9 . At the HL-LHC, at design luminosity of 3 ab −1 , the predicted N Z turns out to be 7.5 × 10 10 , that is, ten times more. The crucial ingredient in searching in proton collisions for the decay Z → γγ is the separation of the signal from the background processes, which is in general very challenging.
The three processes taken into account as main background sources are Regarding process (11a), the unreconstructed energy from jet clustering can mimic missing energy coming from aγ. In process (11b), each neutrino ν has the same signature of a massless dark photon, appearing as missing momentum in the Z → γ + X final state. The same holds true for the neutrinos in the process (11c), where electrons are wrongly reconstructed as photons and pass the photon identification requirements [21]. Both signal and background events were generated at the leading order (LO) using Mad-Graph5 aMC@NLO, which allows to compute the full hard process matrix element, including spin correlations. The Zγγ interaction was modeled adapting the UFO [22] from [23], which in our notation corresponds to assuming C i = C in equation (3) by fixing three out of the nine d f M couplings, thus effectively simulating a simplified version of the full Lagrangian in equation (3) with one single free parameter C, whose value was fixed during simulation to get conventionally BR(Z →γγ) = 0.5.
Parton shower and hadronization effects were simulated using Pythia8 [24,25], while the ATLAS detector response simulation was performed using Delphes [26]. The simulated processes and the corresponding number of generated events and cross-sections are shown in Table I.
In order to optimize the generation step, all samples were produced by requiring a single isolated photon with p T > 10 GeV and a minimum parton transverse momentum p T > 20 GeV in the hard process, where relevant.
In the same fashion, three distinct samples were simulated for the process p p → γ+jets, referred to in Table  Process Slice Nsim σ (pb) Selection Eγ > 300 GeV p p → γ+jets III 3860000 (2.468 ± 0.005) × 10 4 p γ T > 30 GeV Table I. Simulated processes with the corresponding number of generated events (Nsim) and cross-sections. The signal crosssection is conventionally set following BR(Z → γγ) = 0.5. Selection cuts for the single hard photon required during generation are also shown, when relevant. The uncertainties quoted on cross-sections are purely statistical.
I and in the following as slices: slice I was generated without a cut in the photon energy, slice II required a minimum photon energy E min γ before detector smearing effects ("particle level") of 300 GeV, while slice III required a minimum photon transverse momentum p min T,γ of 30 GeV. The three slices were eventually merged into a single final sample, properly weighting by the corresponding cross sections events in the overlapping regions.

Event reconstruction
In order to model a detector as close as possible to ATLAS, the following conditions were specifically implemented. Charged particles were assumed to propagate in a magnetic field of 3.8 T enclosed in a cylinder of radius 1.29 m and length 6 m. Photons were reconstructed from clusters of simulated energy deposits in the electromagnetic calorimeter measured in projective towers with no matching tracks. Photons were identified and isolated by requiring the energy deposits in the calorimeters to be within a cone of size ∆R = ∆η 2 + ∆φ 2 = 0.1 (12) around the cluster barycentre. Candidate photons were required to have transverse energy E T > 10 GeV, in order to simulate ATLAS photon efficiency calibrations selections [21], and to be within |η| ≤ 2.5. Table II. Cuts applied to maximize s/ √ b. When deriving limits, the cut on MT is applied only when using Eγ as discriminating variable (see subsection III A 4 ).
Electrons and muons were reconstructed from clusters of energy deposits in the electromagnetic calorimeter matched to a track within ∆R = 0.5, without any  attempt to simulate independenlty the muon detector response. Jets were reconstructed with the anti-k T algorithm [27] with a radius parameter R = 0.5 from clusters of energy deposits at the electromagnetic scale in the calorimeters. This scale reconstructs the energy deposited by electrons and photons correctly but does not include any corrections for the loss of signal for hadrons due to the noncompensating character of the ATLAS calorimeters. A correction to calibrate the jet energy was then applied, such that a sample of hadrons of a given energy is reconstructed with that same energy. Candidate jets were required to have p T > 20 GeV.
The missing transverse energy E miss T was computed as the transverse component of the negative vector sum of the momenta of the candidate physics objects.

Event selection
Only events with at least one reconstructed photon and no jets in the final state were selected, in order to improve discrimination against the main backgrounds. This requirement, together with the loose cuts applied at generation level, will be referred to in the following as "preselection".
In order to maximize the sensitivity of the search we investigated the possible application of several kinematic cuts targeting an increase of the signal over square root of background yields ratio s/ √ b (see Figures 1 and 2). The resulting cuts are summarized in Table II, where the invariant transverse mass M T built with the photon and the missing transverse energy is defined through the formula where p γ T in the transverse component of the photon momentum, and ∆φ(γ, E miss T ) is the angle between the photon and the missing transverse momenta in the transverse plane.
In the following, the cuts reported in Table II will be referred to as "selection".  The cut efficiencies due to preselection and selection are reported in Table III.

Statistical Methods
We use a simple binned likelihood function L(µ, θ) constructed as a product of Poisson probability terms over all bins considered [28]. Either M T or E γ are used as discriminating variable, whichever gives the better limits.
The likelihood function is implemented in the RooStats package [29]. It depends on the signal-strength parameter µ, a multiplicative factor that scales the number of expected signal events, and θ, a set of nuisance parameters (NPs) that encode the effect of systematic uncertainties on the background expectations, which are implemented in the likelihood function as Gaussian constraints. One should notice that, in our conventions, we can identify the parameter µ with the BR(Z → γγ) in the limit µ 1, which always holds throughout the paper when deriving limits. Individual sources of systematic uncertainty are considered to be uncorrelated. The statistical uncertainty of the MC events is not taken into account in L(µ, θ), while increasing statistics of MC samples in specific regions of phase space when needed (see e.g. subsection III A 1).
The test statistic q µ is defined as the profile likelihood ratio whereμ andθ are the values of the parameters that maximize the likelihood function (with the constraint 0 ≤ µ ≤ µ), andθ µ are the values of the NPs that maximize the likelihood function for a given value of µ. The test statistic q µ is implemented in the RooFit package [30].
Assuming the absence of any significant excess above the background expectation, upper limits on BR(Z → γγ) are derived by using q µ and the CL s method [31,32].

B. Results
We can now derive the upper limit on the BR(Z → γγ) attainable at the LHC.
Let us first look at the best case scenario in which systematic uncertainties are negligible. We find Though the BR in equation (15) might compete with the LEP result BR(Z → γγ) < 10 −6 [12] after taking into account the combination of the ATLAS and CMS experiments, or assuming HL-LHC luminosities, this result is weakened by the effect of systematic uncertainties which are unavoidable in the real case.

Observable
Systematic uncertainty ci E miss Table IV. Main systematic uncertainty sources on the background yields and corresponding overall relative impact. The sources of systematic uncertainties were taken from [34].
The level of uncertainty on the background yields will depend on several factors, possibly decreasing with the life of the accelerator due to welcome efforts of the collaborations in understanding and constraining systematic effects. Here an estimate of the overall relative impact c i = ∆b i /b of the systematic i on the total background yields is attempted, by analysing a subset of possible systematic uncertainty sources chosen, among the set affecting a search with analogue signature [34], as the three with highest impact: the uncertainty on E miss T , the jet energy scale uncertainty and the theoretical uncertainty on the modeling of σ b . In order to do this, we varied up and down the source of the uncertainty by one standard deviation, taken from [34]: we eventually computed the resulting value of c i using the highest variation ∆b i = max ∆b up i , ∆b down i , as summarized in Table IV. A dedicated uncertainty on F e→γ , the e → γ fake-rate, was assigned to the background from electrons wrongly reconstructed and misidentified as photon candidates [21]. This systematic uncertainty conservatively covers both converted and unconverted photon fake rates in the central region.
After properly assigning a NP to each source of systematic uncertainty from Table IV in the likelihood function of subsection III A 4, we compute the 95% CL upper limit on BR(Z → γγ) considering L = 140 fb −1 to be BR(Z → γγ) < 8.0 × 10 −6 .
The M T distribution after selection cuts, with signal yields normalized according to equation (16), is shown in Figure 3. This result, compared with equation (15), highlights how the study and improvement in the control of systematic uncertainties will be a key feature at the LHC. The corresponding HL-LHC upper limit on BR(Z → γγ) can be estimated under the assumption that the systematic uncertainties will decrease by a factor 1/ √ L. Therefore the upper limit on BR(Z → γγ) in equation (16) can simply be multiplied by a factor L LHC /L HL-LHC = 140 fb −1 /3000 fb −1 , obtaining which represents the estimate for the best upper limit reachable by a single experiment at the HL-LHC. From this discussion, it is clear that a search program for the process pp → Z → γγ will hardly compete with the LEP result, mainly because of the large background contamination. This negative result does not come unexpected but the exercise of this section is still useful in showing quantitatively the kind of problems one encounters in trying to study this particular process at a hadron collider.
More promising results can be obtained at future lepton colliders, to which we now turn.

IV. FUTURE LEPTON COLLIDERS
Two circular lepton colliders have currently been proposed. The first one is part of the Future Circular Collider project (FCC), whose integrated program foresees operations in two stages: initially an electron-positron collider (FCC-ee) serving as a Higgs and electroweak factory running at different CM energies, followed by a proton-proton collider at a collision energy of 100 TeV.
The FCC-ee [35] is a high-luminosity, high-precision, 100 km circumference storage ring collider, designed to provide e + e − collisions with centre-of-mass energies from 88 to 365 GeV. The CM operating points with most physics interest are around 91 GeV (Z-boson pole), 160 GeV (W ± pair-production threshold), 240 GeV (ZH production) and 340-365 GeV (tt threshold and above). The machine should deliver peak luminosities above 10 34 cm −2 s −1 per experiment at the tt threshold and the highest ever luminosities at lower energies, with an expected total integrated luminosity L = 150 ab −1 at the Z pole.
The other planned electron-positron machine is the Circular Electron Positron Collider (CEPC [36]), which is expected to collide electrons and positrons at different CM energies: 91.2 GeV, 160 GeV and 240 GeV [37]. This machine is expected to provide a total integrated luminosity L = 16 ab −1 at the Z pole.
These new lepton colliders will produce huge statistics samples of events (of the order of Tera Z bosons), allowing many measurements with unprecedented accuracy, and the discovery and study of rare Z, W , Higgs boson and top decays. Besides offering the ultimate investigations of electroweak symmetry breaking, these precision measurements will be highly sensitive to the possible existence of yet unknown particles, with masses up to about 100 TeV. Sensitive searches for particles with couplings much smaller than weak, such as sterile neutrinos, can be envisioned as well.
We look into the feasibility of a search for a massless dark photon at these new machines, focusing on the center of mass energy of 91.2 GeV, which represents the most promising setup for the process here considered.

An FCC-ee detector: IDEA
The IDEA detector concept [35], developed specifically for FCC-ee, is based on established technologies resulting from years of R&D. Additional work is, however, needed to finalise and optimise the design.
The detector comprises a silicon pixel vertex detector, a large-volume extremely-light short-drift wire chamber surrounded by a layer of silicon micro-strip detectors, a thin, low-mass superconducting solenoid coil, a preshower detector, a dual-readout calorimeter, and muon chambers within the magnet return yoke. Electrons and muons with momenta above 2 GeV and unconverted photons with energies above 2 GeV can be identified with efficiencies of nearly 100% and with negligible backgrounds. The photon energy should be measured to a precision of 11%/ √ E ⊕ 1%.

The CEPC detector
Two primary detector concepts have been developed for the construction of the CEPC detector: a baseline, with two approaches for the tracking systems, and an alternative one, with a different strategy for meeting the jet resolution requirements. In the following we briefly describe the baseline approach, being the one used for simulations.
The baseline detector concept [37] incorporates the particle flow principle with a precision vertex detector, a Time Projection Chamber (TPC) and a silicon tracker, a high granularity calorimetry system, a 3 Tesla superconducting solenoid followed by a muon detector embedded in a flux return yoke. In addition, five pairs of silicon tracking disks are placed in the forward regions at either side of the Interaction Point (IP) to enlarge the tracking acceptance from | cos θ| < 0.99 to | cos θ| < 0.996.
The performance of the CEPC baseline detector concept have been investigated with full simulation, as summarized in the following with focus solely on features which have some impact in the study we are presenting here. Electrons and muons with momenta above 2 GeV and unconverted photons with energies above 5 GeV can be identified with efficiencies of nearly 100% and with negligible backgrounds. The photon energy should be measured to a precision better than 20%/ √ E ⊕ 1%. The ionization energy loss (dE/dx) will be measured in the TPC, allowing the identification of low momentum charged particles. Combining the measurements from the silicon tracking system and the TPC, the track momentum resolution will reach ∆(1/p T ) ∼ 2 × 10 −5 GeV −1 .
A. Methods

Monte Carlo simulation samples
The fiducial cross section for the SM Z-boson production at lepton colliders with a CM energy of 91.2 GeV and, following the same procedure described in subsection III A 1, with decays within |η| < 3 is σ(e + e − → Z/γ * → ff ) = (6.19 ± 0.01) × 10 4 pb . (18) The signal process e + e − → Z → γγ has a distinctive experimental signature. Both the photon and the massless dark photon are monochromatic with an energy of M Z /2 at the Z-factory. The massless dark photon has a neutrino-like signature, appearing as missing momentum in the Z → γ + X 0 final state, in association with a peak of photon events around the mentioned energy values.
Two main processes can contribute to background events: e + e − → γ νν, and (19a) In process (19a) the photon is the result of initial state radiation by either the electron or the positron, and the νν pair is produced either by the decay of a Zboson produced in the s-channel or by W -exchange in the t-channel. The radiative Bhabha reaction of process (19b) contribute to background events when both the final state electron and positron escape detection. However the number of background events from this process strongly depends on the geometry of the detector and on the presence of un-instrumented regions, as observed at LEP [12] and shown in the following sections for the CEPC and FCC-ee detectors.
The hard process, parton-shower and hadronization steps were simulated closely following the lines of subsection III A 1. We checked explicitly that relevant kinematic distributions from reconstruction of events with the two different detector configurations of FCC-ee and CEPC closely match in event yields, up to Monte Carlo statistical fluctuations, the minor differences coming only from the slightly different energy resolutions. If not stated explicitly, in the following the configuration for the baseline CEPC detector concept must be understood, with results from the two detectors differing only by the corresponding integrated luminosity.
The simulated processes and the corresponding number of generated events and cross-sections assuming collisions between positron and electron beams with E beam =45.6 GeV are reported in Table V. For each background process, two slices have been simulated with different E min γ during generation, respectively of 18 and 30 GeV for the e + e − → γνν and e + e − → γe + e − process (see subsection III A 1).

Event reconstruction
Charged particles were assumed to propagate in a magnetic field of 3.5 T homogeneously filling a cylindrical region of radius 1.81 m and length of 4.70 m.
Photons were reconstructed from clusters of energy deposits in the electromagnetic calorimeter measured in projective towers with no matching tracks. Photons were identified and isolated by requiring the energy deposits in the calorimeters to be within a cone of size ∆R = ∆η 2 + ∆φ 2 = 0.5 (20) around the cluster barycenter. Candidate photons were required to have E > 2 GeV and to be within | cos θ| ≤ 0.987. Electrons and muons were reconstructed from clusters in the electromagnetic calorimeter with a matching a track. The criteria for their identification were similar to those used for photons.
Jets were reconstructed with the anti-k T algorithm with a radius parameter R = 0.5 from clusters of energy deposits in the calorimeters (hadronic and electromagnetic). Candidate jets were required to have p T >20 GeV.
The missing energy vector p miss was computed as the negative vector sum of the momenta of the candidate physics objects.  Photon energy Eγ distribution after the cut in η relative to the signal (e + e − → γγ) and to the background processes (e + e − → γνν and e + e − → γe + e − ).The signal distribution is normalized to the total estimated background yield.
Events were initially selected to have at least one photon and no charged particles in the final state. This requirement, together with the loose cuts applied at generation level, will be referred to in the following as "preselection". Selection cuts were then applied to increase the ratio s/ √ b and improve the upper limit on the BR(e + e − → γγ), as summarized in Table VI.
At leptonic colliders the initial state of the process is fully determined (apart from initial-state-radiation ef-  Table V. Simulated processes with the corresponding number of generated events (Nsim) and cross-sections. As in Table I, the signal cross-section is conventionally set following BR(Z → γγ) = 0.5. See subsection III A 1 for the definition of the slices I and II.
fects) by the beam parameters, contrary to the hadron collider case, and the center of mass frame coincides with the laboratory frame. Yet we find instructive, in order to better understand the results of section III B, to keep track of the same kinematical observables defined therein.
The transverse invariant mass M T in events with one single photon and missing energy simplifies at particle level to the expression The M T and photon energy E γ distributions after selection are shown in Figures 4 and 5. The preselection and selection efficiencies are reported in Table VI: no event from the process e + e − → γe + e − passed selection requirements, as expected from kinematic arguments, thanks to the absence of uninstrumented regions in both CEPC and FCC-ee detector designs.

B. Results
A lepton collider working at a a CM energy of 91.2 GeV is an actual Z-factory and the ideal place where to look for the small values of the BR(Z → γγ) we are after.
Let us consider first the case of no systematic uncertainties. The lepton colliders perform very well. The CEPC, with a total integrated luminosity L = 16 ab −1 yields The FCC-ee, with an expected total integrated luminosity of L = 150 ab −1 , gives The results in equation (22) and equation (23) are not modified by taking into account the uncertainties on both the σ f id Z (0.01 × 10 4 pb) and the luminosity ∆L/L = 10 −4 [37]. Systematic effects on the upper limit on BR(e + e − → γγ) are found to be negligible.
The results summarized in Table VII show that, at both the CEPC and FCC-ee lepton colliders, running at the Z pole mass, the limit on the BR(e + e − → γγ) will significantly improve the present LEP bound.

V. SPIN ANALYSIS
Having discussed the discovery potential of a dark photon produced in association with a photon in Z-boson decays, in this last section we investigate how to establish the spin of such a new neutral state. Since nothing prevents the Z boson to decay into a photon and a pseudoscalar massless neutral particle a, the latter can be used as test hypothesis against the J P = 1 − nature of the dark photon.
No attempt is made here to optimize the search strategy for a pseudoscalar signal Z → aγ, meaning that results presented in previous sections do not necessarily apply to this test hypothesis.

A. Methods
We now assume that the discovery ofγ in Z-boson decays has occured at a future e + e − collider with √ s = M Z . In this scenario, a good observable discriminating between the two J P = 1 − , 0 − hypotheses is the cosine of the polar angle θ of the detected photon [13]. The corresponding distributions ( Figure 6) can be produced using the linear realization of the model [38], where the two relevant dimension-five operators regarding γa final states at the Z peak are with the latter contributing through interference with the first. We have checked that the independent contributions from O aZγ , O aγγ and their interference in the cos θ distribution are indistinguishable in shape, as expected. One can then apply the following analysis to any combination of the cB and cW couplings, and in particular for the common and phenomenologically appealing cB = − tan θ W cB choice [38], which corresponds to set g aγγ = 0. Following the statistical treatment described in subsection III A 4, a likelihood function L(J P , µ, θ) that depends on the spin-parity assumption of the signal is constructed as a product of conditional probabilities over the binned distribution of the discriminating observable cos θ: where, given the clean lepton collider environment, a good starting approximation is to assume the best case scenario in which all nuisance parameters θ are sufficiently constrained by auxiliary measurements through the functions A(θ), such to allow to safely neglect the contribution of systematic uncertainties on the background, which we assume subtracted from total yields with dedicated methods (see e.g. [39]).
Closely following [40], a proper test statistic q is then chosen to be the logarithm of the likelihood ratio In this case, given the rather simple form of both equation (26) and (27), they have been implemented directly in a ROOT [41] script. The distributions of the test statistic q for both signals shown in Figure 7 were obtained, as an example, using n toys = 160000 pseudoexperiments and with the specific choice of N = 10 signal events.
The distributions of q are used to determine the corresponding p 0 values p 0 (J P = 1 − ) and p 0 (J P = 0 − ). For instance, the tested hypothesis p 0 (J P = 0 − ), the expected and the observed p 0 value is obtained by integrating the corresponding test-statistic distribution, respectively, above the J P = 1 − q distribution median and above the observed value of q. The exclusion of the J P = 0 − hypothesis in favor of the dark photon J P = 1 − hypothesis is evaluated in terms of the corresponding CL s (0 − ), defined as: In the following, we always assume the observed value of the test statistics to be q = 0.

B. Results
The example described in Figure 7 gives an expected p 0 (J P = 0 − ) value of 3.9×10 −3 and an observed value of  gives an expected exclusion of the J P = 0 − hypothesis at the 99% CL, whereas an observed value of p 0 (J P = 1 − ) = 1.5 × 10 −1 gives an observed exclusion at 89% CL.
Repeating the procedure for all values in the range N ∈ [2, 25] we estimate the lower bound for the expected and observed number of signal events needed to exclude the p 0 (J P = 0 − ) test hypothesis under the p 0 (J P = 1 − ) assumption to be, respectively, N = 6 and N = 17 at the 95% CL.

VI. SUMMARY AND OUTLOOK
The Z-boson decay into a SM photon and a dark photon would be a most striking signature for the existence of dark photons.
The Z → γγ experimental signature is quite simple and distinctive. In the Z-boson CM frame, both the photon and the dark photon are monochromatic with an energy of M Z /2. A massless dark photon has a neutrinolike signature in a typical experiment, and appears as missing momentum in the Z → γ + X final state.
In this paper we present the estimate of the best exclusion limits on BR(Z → γγ) for the Z decay into a photon and a dark photon, comparing several present and future collider scenarios: the LHC (at a CM energy of 13 TeV with an integrated luminosity of 140 fb −1 ), the HL-LHC with 3000 fb −1 , and the two future circular leptonic machines FCC-ee and CEPC, at the specific design CM energy of 91.2 GeV (Z-factory). A summary of the 95% confidence level upper limits on the BR(Z → γγ) obtained under the expected conditions of these colliders are collected in Table VII. It is then straightforward to compare these results with the present LEP bound BR(Z → γγ) < 10 −6 at the 95% CL.
The impact of systematic uncertainties at the LHC and the challenge of large QCD backgrounds, intrinsic to hadron colliders, make all but impossible to match the LEP performance. Only a search at the HL-LHC could yield a result competing with the LEP limit.
At the Z-factories possibly realized at FCC-ee and/or CEPC, instead, limits better than the LEP one could be obtained, thanks to higher luminosities and efficiencies. The much cleaner environment could allow a sensitivity to BR(Z → γγ) of order O(10 −11 ). Such a value comes close to those predicted in dark-sector models where the effective coupling of the dark photon to the Z-boson can be computed. BR(Z → γγ) √ s L (ab −1 ) MT Eγ LHC 13 TeV 0.14 8 × 10 −6 5 × 10 −5