Search for vector mediator of Dark Matter production in invisible decay mode

A search is performed for a new sub-GeV vector boson ($A'$) mediated production of Dark Matter ($\chi$) in the fixed-target experiment, NA64, at the CERN SPS. The $A'$, called dark photon, could be generated in the reaction $ e^- Z \to e^- Z A'$ of 100 GeV electrons dumped against an active target which is followed by the prompt invisible decay $A' \to \chi \overline{\chi}$. The experimental signature of this process would be an event with an isolated electron and large missing energy in the detector. From the analysis of the data sample collected in 2016 corresponding to $4.3\times10^{10}$ electrons on target no evidence of such a process has been found. New stringent constraints on the $A'$ mixing strength with photons, $10^{-5}\lesssim \epsilon \lesssim 10^{-2}$, for the $A'$ mass range $m_{A'} \lesssim 1$ GeV are derived. For models considering scalar and fermionic thermal Dark Matter interacting with the visible sector through the vector portal the 90% C.L. limits $10^{-11}\lesssim y \lesssim 10^{-6}$ on the dark-matter parameter $y = \epsilon^2 \alpha_D (\frac{m_\chi}{m_{A'}})^4 $ are obtained for the dark coupling constant $\alpha_D = 0.5$ and dark-matter masses $0.001 \lesssim m_\chi \lesssim 0.5 $ GeV. The lower limits $\alpha_D \gtrsim 10^{-3} $ for pseudo-Dirac Dark Matter in the mass region $m_\chi \lesssim 0.05 $ GeV are more stringent than the corresponding bounds from beam dump experiments. The results are obtained by using tree level, exact calculations of the $A'$ production cross-sections, which turn out to be significantly smaller compared to the one obtained in the Weizs\"{a}cker-Williams approximation for the mass region $m_{A'} \gtrsim 0.1$ GeV.

A search is performed for a new sub-GeV vector boson (A ) mediated production of Dark Matter (χ) in the fixed-target experiment, NA64, at the CERN SPS. The A , called dark photon, can be generated in the reaction e − Z → e − ZA of 100 GeV electrons dumped against an active target followed by its prompt invisible decay A → χχ. The experimental signature of this process would be an event with an isolated electron and large missing energy in the detector. From the analysis of the data sample collected in 2016 corresponding to 4.3 × 10 10 electrons on target no evidence of such a process has been found. New stringent constraints on the A mixing strength with photons, 10 −5 10 −2 , for the A mass range m A 1 GeV are derived. For models considering scalar and fermionic thermal Dark Matter interacting with the visible sector through the vector portal the 90% C.L. limits 10 −11 y 10 −6 on the dark-matter parameter y = 2 αD( mχ m A ) 4 are obtained for the dark coupling constant αD = 0.5 and dark-matter masses 0.001 mχ 0.5 GeV. The lower limits αD 10 −3 for pseudo-Dirac Dark Matter in the mass region mχ 0.05 GeV are more stringent than the corresponding bounds from beam dump experiments. The results are obtained by using exact tree level calculations of the A production cross-sections, which turn out to be significantly smaller compared to the one obtained in the Weizsäcker-Williams approximation for the mass region m A 0.1 GeV.

I. INTRODUCTION
Despite the intensive experimental searches dark matter (DM) still is a great puzzle. The difficulty so far is that DM can be probed only through its gravitational interaction with visible matter. An exciting possibilities is that in addition to gravity, a new force between the dark and visible matter transmitted by a new vector boson, A , called dark photon, might exist [1][2][3][4]. The A can * Corresponding author, Sergei.Gninenko@cern.ch have a mass m A 1 GeV, and couple to the standard model (SM) via kinetic mixing with the ordinary photon, described by the term 2 F µν F µν and parameterized by the mixing strength . The Lagrangian of the SM is extended by the dark sector in the following way: where the massive vector field A µ is associated with the spontaneously broken U D (1) gauge group, F µν = ∂ µ A ν − ∂ ν A µ , e D is the coupling constant of the U (1) D arXiv:1710.00971v2 [hep-ex] 22 Mar 2018 gauge interactions, and m A , m χ are the masses of the dark photon and DM particles, respectively. Here, we consider as an example the Dirac spinor field, χ, which is treated as Dark Matter fermions coupled to A µ by the dark portal coupling constant e D . The mixing term of (1) results in the interaction: of dark photons with the electromagnetic current J µ em with a strength e, where e is the electromagnetic coupling and 1 [5][6][7]. Such small values of can be obtained naturally in GUT from loop effects of particles charged under both the dark and SM U (1) interactions with a typical 1-loop value = ee D /16π 2 10 −2 − 10 −4 [7], while 2-loop contributions result in the range 10 −3 −10 −5 . An additional hint for the existence of A is suggested by the 3.6 σ deviation from the SM prediction of the muon anomalous magnetic moment g µ − 2 [8], which can be explained by a sub-GeV A with coupling 10 −3 [9][10][11], as well as by hints of astrophysical signals of DM [3]. This has motivated a worldwide experimental and theoretical effort towards dark forces and other portals between the visible and dark sectors, see Refs. [4,[12][13][14][15] for a review.
Since there are no firm predictions for the A , its experimental searches have been performed over a wide range of A masses and decay modes. If the A is the lightest state in the dark sector, then it would decay mainly visibly, i.e., typically to SM leptons l (or hadrons) with the rate given by which can be used to detect it. Here, α = e 2 /4π and m l is the lepton mass. Such dark photons in the mass region below a few GeV has been mainly searched for in beam dump, fixed target, collider and rare meson decay experiments, which already put stringent limits on the mixing 2 10 −7 of such dark photons excluding, in particular, the parameter region favored by the g µ − 2 anomaly [15]- [23].
However, in the presence of light dark states χ, in particular, DM with the masses m χ < m A /2, the A would predominantly decay invisibly into those particles provided that coupling g D > e. The decay rate of A →χχ in this case is given by Various dark sector models motivate sub-GeV scalar and Majorana or pseudo-Dirac fermion DM coupled to dark photons [14,15,[24][25][26][27][28][29][30]. To interpret the observed abundance of thermal relic density, the requirement of the thermal freeze-out of DM annihilation into visible matter through γ − A kinetic mixing allows one to derive a relation among the parameters α D 0.02f 10 −3 2 m A 100 M eV 4 10 M eV m χ 2 (5) where α D = e 2 D /4π, f 10 for a scalar [24], and f 1 for a fermion [25]. This prediction combined with the fact that the intrinsic scale of the dark sector could be smaller than, or comparable to, that of the visible sector, provide an important target for the ( , m A ) parameter space which can be probed at energies attainable at the CERN SPS. Models introducing such invisible A also offered possibilities to explain the g µ − 2 and various other anomalies [31] and are subject to different experimental constraints [32][33][34][35]. The severe limits on invisible decays of sub-GeV A s have been obtained from the results of beam dump experiments LSND [24,36] and E137 [37], under assumptions of certain values of the coupling strength, α D , and masses of the DM decay particles. Recently, NA64 [38] and BaBar [39] experiments set new direct limits, 2 10 −6 for m A 100 MeV and 2 10 −6 for m A 1 GeV, respectively, which rule out the A parameter space explaining the muon g-2 anomaly, leaving, however, a significant area that is still unexplored.
In the following we assume that the Dark Matter invisible decay mode is dominant, i.e. Γ(A →χχ)/Γ tot 1, and that the A leptonic decay channel is suppressed, Γ(A →χχ) Γ(A → e − e + ). If such A exists, many crucial questions about its mass scale, coupling constants, decay modes, etc. arise. One possible way to answer these questions, is to search for the invisible A in accelerator experiments. The A s could be produced in a high intensity beam dump experiment and generate a flux of DM particles through their decays, which can be detected through the scattering off electrons in the detector target [24,25,33,36,[40][41][42]. In this case, the signal event rate in the detector scales as 2 y ∝ 4 α D , where the parameter y is defined as which was recently constrained to 10 −9 y 10 −8 , for α D = 0.5 and for dark-matter masses of 0.01 < m χ < 0.3 GeV by the MiniBooNE experiment [43]. Another approach, considered in this work and proposed in Refs. [44,45], is based on the detection of the large missing energy, carried away by the energetic A produced in the interactions of high-energy electrons in the active beam dump target, see also [25]. The advantage of this type of experiment is that its sensitivity is proportional to the mixing strength squared, 2 , associated with the A production in the primary reaction and its subsequent prompt invisible decay, while in the former case it is proportional to 4 α D , with 2 associated with the A production in the beam dump and 2 α D coming from the χ particle interactions in the detector.
In this work we report new results on the search for the A and light DM in the fixed-target experiment NA64 at the CERN SPS. The experimental signature of events from the A → invisible decays is clean and they can be selected with small background due to the excellent capability of NA64 for the precise identification and measurements of the initial electron state.
The rest of the paper is organised as follows. Section II outlines the method of search and theoretical setup for the A production in an electron-nuclei scattering, and the signal simulation. Here, we mainly focus on the experimental signature of the A → invisible decays and A production rate. We also attempt to provide an estimate of the experimental uncertainties associated with the A cross section calculation required for the sensitivity estimate. We revisit here the calculations of Refs. [44][45][46] and clarify the apparent disagreements in the numerical factors in the cross section of A production in the Weizsäcker-Williams framework and exact computations at tree level. We also discuss additional experimental inputs that would be useful to improve the reliability of the calculated sensitivity of the NA64 experiment. The H4 beam line and experimental set-up is presented in Sec. III, followed by a description of the event reconstruction and analysis in Sec. IV. The results on the benchmark process of dimuon production are presented in Sec.V. In Sec. VI and VIII the signal efficiency and background sources are discussed. The final results on the searches for invisible decays of dark photons and light thermal DM are reported in Sec. IX and X, respectively. We present our conclusions in Sec. XI.

II. METHOD OF SEARCH AND THE A PRODUCTION
As seen from the Lagrangian (1), any source of photons will produce all kinematically possible massive A states according to the appropriate mixing strength. If the coupling strength α D and A masses are as discussed above, the A will decay predominantly invisibly.
The method of the search for the A → invisible decay is as follows [44,45]. If the A exists it could be produced via the kinetic mixing with bremsstrahlung photons in the reaction of high-energy electrons absorbed in an active beam dump (target) followed by the prompt A → invisible decay into DM particles in a hermetic detector: see Fig. 1. A fraction f of the primary beam energy E A = f E 0 is carried away by χ particles, which penetrate the target and detector without interactions resulting in zero-energy deposition. The remaining part of the beam energy E e = (1 − f )E 0 is deposited in the target by the scattered electron. The occurrence of the A production via the reaction (7)  above those expected from backgrounds. Here we assume that in order to give a missing energy signature the χs have to traverse the detector without decaying visibly. No other assumptions are made on the nature of the A → invisible decay . In previous work [38,46], the differential cross-section A -production from reaction (1) was calculated with the Weizsäcker-Williams (WW) approximation, see [47,48]. The cross-sections were implemented a Geant4 [49,50] based simulations, and the total number n A of the produced A per single electron on target (EOT), depends in particular on , m A , E 0 and was calculated as where ρ is density of Pb target, N A is the Avogadro's number, A Pb is the Pb atomic mass, n(E 0 , E e , s) is the number of e ± with the energy E e in the e-m shower at the depth s (in radiation lengths) within the target of total thickness T , and σ(E e ) is the cross section of the A production in the kinematically allowed region up to E A E e by an electron with the energy E e in the elementary reaction (7). The energy distribution dn A dE A of the A s was calculated by taking into account the differential cross-section dσ(Ee,E A ) dE A , as described in Ref. [46]. The numerical summation in Eq. (8) was performed with the detailed simulation of e-m showers done by Geant4 over the missing energy spectrum in the target, see Fig. 4. According to the simplified WW approximation [47] the e − N scattering total rate can be written as Here, t min = m 4 A /(4E 2 0 ) and t max = m 2 A are approximated values of the A momentum transfer. For most energies the elastic form-factor G 2,el (t) dominates and can be approximated as where a = 111 Z −1/3 /m e and d =0.164 GeV 2 A −2/3 . Note that for heavy atomic nuclei A one has to take into account the inelastic nuclear form factor. The flux is given by Φ = Z 2 · Log, where the Log value depends weakly on atomic screening, nuclear size effects and kinematics [47,51]. Numerically Log ≈ (5 − 10) for m A ≤ 500 MeV. It has been recently pointed out, that for a certain kinematic region of the parameters m A , E A , the A yield derived in the WW framework could differ significantly from the one obtained with the exact tree-level (ETL) calculations [52,53]. Therefore, it is instructive to perform an accurate calculation of the A cross-section based on precise phase space integration over the final state particles in the reaction e − Z → e − ZA . A reliable theoretical prediction for the A yield is essential for the proper interpretation of the experimental results in order to obtain robust exclusion limits in the A parameter space or the possible observation of the A signal.
In order to derive more accurately the A yield, we have used the A production cross sections of (7) obtained without the WW approximation, but with ETL calculations of Ref. [54]. These cross sections were cross checked with those calculated by Liu et al. [52,53] and were found to be in agreement. The comparison showed that the difference between the two calculations in the wide A mass range does not exceed 10% (see Fig.  1 in Ref. [54]), which will be further used as a systematic uncertainty in the calculation of the A yield. This difference is, presumably, due to the different accuracy of computation programs used for integration over the phase space and, possibly due to the different parameterisation of the form factors used as an input for the cross section calculations.
In order to implement the ETL cross-section formula into Geant4 [50] based NA64 simulation package, we introduce in Eq. (8) a correction k-factor defined by the following ratio Here, the cross-section σ A EXACT takes into account the phase space integration over the final states of the particles and represents the overall uncertainties in the crosssection (9) calculated in simplified WW approach. The A yield was calculated by using (8) with the replacement We refer to this method as k-factor approach throughout the paper. For heavy target nuclei k depends rather weakly on Z and A, i.e. for tungsten and lead the deviation is about 0.5%. In Fig. 2, k values are shown as a function of the electron beam energy E 0 for various m A . One can see that for m A = 100 MeV and E 0 = 100 GeV, the ETL cross sections of the A production are smaller by a factor 1.7 than the corresponding WW crosssections. On the other hand, σ EXACT exceeds σ W W for the masses m A below 5 MeV, so that one can slightly improve the limits on the mixing strength for this mass region. An example of differential A spectra from the electron beam interactions in the thin, X 0 , Pb target (here X 0 is the radiation length) calculated for the mass m A = 100 MeV as a function of x = E A /E e for different electron energies is shown in Fig. 3. Once the A flux (8) was defined, the next step was to simulate the A emission spectrum from the target. The decay electrons and positrons were tracked through the dump medium including bremsstrahlung photons, their conversion and multiple scattering in the target. The A reconstruction efficiency in the target was computed and convoluted with the target details and detector geometrical acceptance (see Sec. III) based on the NA64 Monte Carlo (MC) simulation package used in our previous search [38]. The comparison of the energy distributions of A s emitted from the thick target (t X 0 ) with the energy E A 0.5E 0 calculated for masses m A = 10 and 100 MeV with and without WW approximation for the 100 GeV beam energy is shown in Fig. 4. The spectra have the similar shape and differ mostly in the overall nomalization factor. Note that these distributions represent also the spectra of the missing energy in the detector.

III. H4 BEAM AND NA64 DETECTOR
The experiment employs the optimized 100 GeV electron beam from the H4 beam line at the North Area (NA) of the CERN SPS described in details in Ref. [55]. The H4 provides an essentially pure e − beam for fixedtarget experiments. The beam was designed to transport the electrons with the maximal intensity up to 10 7 per SPS spill of 4.8 s in the momentum range between 50 and 150 GeV/c that could be produced by the primary proton beam of 400 GeV/c with the intensity up to a few 10 12 protons on a beryllium target. The main contribution to the e − yield from the target was the production of π 0 followed by a process π 0 → γγ → e + e − . The short-lived π 0 decays inside the target, and the electrons are produced through the conversion of the decay photons in a separate converter [56]. Protons and charged secondaries that did not interact in the convertor are separated from the neutrals by deflecting them in a mag-netic field to a thick absorber. The electrons produced in the converter are transported to the NA64 detector inside an evacuated beam-line tuned to an adjustable beam momentum. The hadron contamination in the electron beam was π/e − 10 −2 . The beam has the transverse size at the detector position of the order of a few cm 2 and a halo with intensity a few %.
The signal event recognition in NA64 must rely on the detection of the incoming and outgoing electron only, since the decay product for the A → invisible decay are undetectable. The NA64 detector, which is located at about 500 m from the proton target, is schematically shown in Fig. 5. The setup utilized the beam defining scintillator (Sc) counters S1-S3 and veto V1, and the spectrometer consisting of two successive dipole magnets with the integral magnetic field of 7 T·m and low-material-budget tracker. The tracker was a set of two upstream Micromegas chambers (MM1,2) and two downstream MM3,4 and GEM1,2 stations, measuring the beam e − momenta, P e , with the precision δP e /P e 1% [57]. The Straw Tubes chambers (ST) were used for calibration purposes. The in (out-)coming electron azimuthal angle was tuned to be within θ in(out) 1 mrad with respect to the primary beam axis. A small fraction of events with the larger incoming angle in the range θ in(out) 1 − 20 mrad which was typically correlated with the smaller track momentum was rejected by further analysis. The magnets also served as an effective filter rejecting the low energy electrons present in the beam. To improve the high energy electron selection and suppress background from a possible admixture of low energy electrons, a tagging system utilizing the synchrotron radiation (SR) from high energy electrons in the magnetic field was used, as shown schematically in Fig. 5. The basic idea was that, since the SR energy emitted by a particle per revolution with a mass m and energy E 0 is < E SR >∝ E 3 0 /m 4 , the low energy electrons and hadrons in the beam could be effectively rejected by using the cut on the energy deposited in the SR detector (SRD) [44,58]. A 15 m long vacuum vessel was installed between the magnets and the ECAL to minimize absorption of the SR photons detected immediately at the downstream end of the vessel with a SRD, which was an array of PbSc sandwich calorimeter with fine longitudinal segmentation. Compared to the previous measurements [38], the SRD was also segmented transversely by three SRD counters, each 60 × 80mm 2 in lateral size assembled from 80 − 100 µm Pb and 1 mm Sc plates with wave length shifting (WLS) fiber read-out. This allowed to additionally suppress background from hadrons, that could knock off electrons from the output vacuum window of the vessel producing a fake e − SRD tag, by about two orders of magnitude [58]. The detector was also equipped with an active target, which was a hodoscopic electromagnetic calorimeter (ECAL) for the measurement of the electron energy deposition, E ECAL , with the accuracy δE ECAL /E ECAL 0.1/ E ECAL [GeV] as well as the X, Y coordinates of the incoming electrons by us- ing the transverse e-m shower profile. The ECAL was a matrix of 6 × 6 Shashlik-type counters assembled with Pb and Sc plates with WLS fiber read-out. Each module was 40 radiation lengths (X 0 ) and had an initial part 4 X 0 used as a preshower (PS) detector. By requiring the presence of in-time SR signal in all three SRD counters, and using information of the longitudinal and lateral shower development in the ECAL, the initial level of the hadron contamination in the beam π/e − 10 −2 was further suppressed by more than 4 orders of magnitudes, while keeping the electron ID efficiency at the level 95% [58]. The ECAL PMTs were read-out with sampling ADC (MSADC) electronics which consist of shapers and the ADCs themselves [59,60]. The shaper stretched the PMT signal to 100 ns while the MSADC sampled the signal amplitude every 12.5 ns. As described below this allowed using an algorithm to extract precise timing and amplitude from the MSADC information in the presence of pileups at high intensity. The NA64 Data Acquisition system was adapted from the one used in the COM-PASS experiment at CERN [61]. A high-efficiency veto counter V 2 , and a massive, hermetic hadronic calorimeter (HCAL) of 30 nuclear interaction lengths (λ int ) were positioned just after the ECAL. The V 2 was a plane of scintillation counters used to veto charged secondaries incident on the HCAL detector from upstream e − interactions. The HCAL which was an assembly of four modules HCAL1-HCAL4, served as an efficient veto to detect muons or hadronic secondaries produced in the e − A interactions in the ECAL target. Each module was a sandwich of 48 alternating layers of iron and scintillator (Sc) with a thickness of 25 mm and 4 mm, respectively, with a total length of 7λ int , and with a lateral size of 60 × 60 cm 2 . Each Sc layer consisted of 3×3 plates with WLS fiber readout allowing to assemble the whole HCAL module as a matrix of 3 × 3 cells, each of 20 × 20 cm 2 . The number of photoelectrons produced by a minimum ionizing particle (MIP) crossing the single module was in the range 150-200 photoelectrons. The HCAL energy resolution was δE HCAL /E HCAL 0.6/ √ E HCAL [GeV]. The single electron events were collected with the hardware trigger designed to accept events with in-time hits in beamdefining counters S i and clusters in the PS and ECAL with the energy thresholds E th P S 0.3 GeV and E th ECAL 80 GeV, respectively. The missing energy events have the signature with the incoming track momentum P e 100 GeV, and V 2 and HCAL zero-energy deposition, defined as energy below the thresholds E th V2 1 MIP and E th HCAL 1 GeV, respectively.

IV. DATA ANALYSIS AND SELECTION CRITERIA
The search for the A → invisible decay described in this paper uses the full data sample collected during July and October runs in 2016 corresponding to n EOT = 4.3× 10 10 EOT. The results reported here are obtained using three sets of data in which n EOT = 2.3 × 10 10 , 1.1 × 10 10 and 0.9×10 10 EOT were collected with the beam intensities (1.4−2)×10 6 , (3−3.5)×10 6 and (4.5−5)×10 6 e − per spill, respectively. Data of these three runs (hereafter called respectively the run I,II, and III) were processed with selection criteria similar to the one used in our previous paper [38] and finally combined as described in Sec. IX. Compared to the analysis of Ref. [38], a number of improvements in the event reconstruction, e.g., adding the pileup algorithm, were made in order to increase the reconstruction efficiency.
The strategy of the analysis was to identify A → invisible candidates by precise reconstruction of the initial e − state and an isolated low energy e-m shower in the ECAL that are accompanied by no other activity in the V 2 and HCAL detectors. The measured rate of such events was then supposed to be compared to that expected from known sources. The spectra of A s produced in the ECAL target by primary electrons were calculated using the approach reported in Ref. [46]. An example of the distributions of energy deposited in the target calculated for the masses m A = 2, 20 and 200 MeV is shown in Fig. 6. A detailed Geant4 based MC simulation was used to study the detector performance and acceptance losses, to simulate background sources, and to select cuts and estimate the reconstruction efficiency.
The candidate events were pre-selected with the criteria chosen to maximize the acceptance for simulated signal events and to minimize the numbers of events expected from background sources discussed in Sec. VIII. The following selection criteria were applied: • There must be one and only one incoming particle track having a small angle with respect to the beam axis. This cut rejects low momentum electrons as they were typically correlated with a large-angle incoming tracks originating presumably from the upstream e − interactions. The reconstructed momentum of the particle was required to be P e = 100 ± 2 GeV.
• The track should be identified as an electron with the SRD detector. The energy deposited in each of the three SRD modules should be within the SR range emitted by e − s and in time with the trigger. This was the key cut identifying the pure initial e − state, with the pion suppression factor < 10 −5 and electron efficiency > 95% [58].
• The lateral and longitudinal shower shape in the ECAL should be consistent with the one expected for the signal shower [46,62]. It is also used to distinguish hadrons from electrons providing an additional hadron rejection factor of 10 [63].
• There should be no activity in the veto counter V 2 .
In total 7 × 10 4 events passed these criteria from the combined 2016 data sample. The final selection is a cut-based and uses the cuts on the ECAL missing energy E miss = E beam − E ECAL and on the energy deposition in the HCAL. In order to avoid biases in the choice of selection criteria for A events, a blind analysis was performed, with a preliminary definition of the signal box as E miss > 50 GeV and E HCAL < 1 GeV. The HCAL zeroenergy threshold E th HCAL = 1 GeV in (15), see Sec.VI was determined mostly by the noise of the read-out electronics. Events from the preliminary signal box were excluded from the analysis of the data until the validity of the background estimate in this region was established. The cut on E miss was optimized as described in Sec.IX.
In Fig. 7 the left panels show an example of the distributions of events from the reaction e − Z → anything in the (E ECAL ; E HCAL ) plane measured in the runs II(top) and III(bottom) with moderate selection criteria requiring only the presence of the SRD tag identifying the beam electrons. Here, E HCAL is the sum of the energy deposited in the HCAL1 and HCAL2. The distributions of events from the run I with low intensity are similar to the one shown in Ref. [38]. Events from the areas I in Fig. 7 originate from the QED dimuon production, dominated by the the muon pair photoproduction by a hard bremsstrahlung photon conversion on a target nucleus: with some contribution from γγ → µ + µ − fusion process. The µ + µ − pairs were characterised by the HCAL energy deposition of 10 GeV. This rare process whose fraction of events with E ECAL 60 GeV was 10 −5 /EOT served as a benchmark allowing to verify the detector performance and as a reference for the background prediction. The regions II shows the events from the SM hadron electroproduction in the ECAL which satisfy the energy conservation E ECAL + E HCAL 100 GeV within the detector energy resolution. The leak of these events to the signal region mainly due to the HCAL energy resolution was found to be negligible. The fraction of events from the region III was due to pileup of e − and beam hadrons. It was beam rate dependent with a typical value from about a few % up to 20%.
To evaluate the performance of the setup, a cross-check between a clean sample of 10 4 observed and MC simulated µ + µ − events was made. The process (16) was used as a benchmark allowing to verify the reliability of the MC simulation and to estimate the corrections to the signal reconstruction efficiency and possible additional uncertainties in the A yield calculations. Let us first briefly review the description of the gamma conversion into a muon-antimuon pair implemented in Geant4. The dimuon production was also used as a reference for the prediction of background, see Sec. VIII.

A. Simulation of dimuon events
The dimuon production has been simulated with Geant4 [50] and a code developed by NA64 used also for simulation of dark photon production [46]. Here, we report our comparison with data based mostly on Geant4 simulation for decays and propagation of muons through the detectors. However, we anticipate that comparison of dimuon results with the NA64 code will follow in the future, as it will be an important additional cross-check of the A yield calculations reported in this work and in Ref. [54].
The gamma conversion into a muon-antimuon pair on nuclei is a well known reaction in particle physics (Bethe-Heitler process). The simulation of this reaction in Geant4 is based on the differential cross section for electromagnetic creation of muon pairs on nuclei (A, Z) in terms of the energy fraction of muons [64,65]: Eγ , α = 1 137 and r µ = α mµ is the classical radius of muon and where Here, A is the atomic number of the nuclei. The differential cross section is symmetric in x + and x − , and a relation is takes place. The differential cross section (18) can be rewritten in the form Here x = x + or x = x − . The total cross section was obtained by integration of the differential cross section, namely where . Numerically for Pb nuclei σ tot = 30.2; 334; 886 µb for E γ = 1, 10, 100 GeV, respectively.
Note, that formula (18) for the cross section was obtained from the tree level formula for the differential cross section γZ → µ + µ − Z by taking into account both the atomic and nuclear form-factors and without using the WW approximation of equivalent photons. Even though the production mechanisms of the A and µ + µ − pair are different, the number of A and dimuon events, are both proportional to the square of the Pb nuclear form factor F (q 2 ) and are sensitive to its shape. As the mass (m A m µ ) and q 2 (q m 2 A /E A m 2 µ /E µ ) ranges of the final state for both reactions are similar, the observed difference can be considered as due to the accuracy of the dimuon yield calculation for heavy nuclei.

B. Yield of dimuon events
The dimuon events were selected with the trigger (14), which accepted only events with the ECAL energy deposition smaller than 80 GeV. Because only muons can punchthrough the total length of the modules HCAL1-3 ( 21λ int ) without interactions, the selection was based on the requirement of the energy deposited in HCAL1 and HCAL4 modules to be in the range 1 E HCAL1, 4 6 GeV, which is comparable with that expected from a single muon or dimuon pair. In Fig. 8  is shown. Here, the HCAL energy is defined as the total energy deposited in the four HCAL modules.
The dimuon yield was estimated from the observed number of reconstructed dimuon events. The comparison of the number of observed (n data  Table I. One can see, that the reconstruction efficiency of µ + µ − pairs were found to be beam rate dependent. The difference between the number of observed and MC predicted µ + µ − events with E ECAL 60 GeV is the range 8-19%. It can be interpreted as due to the inaccuracy of the dimuon yield determination for heavy nuclei target and, can be conservatively accounted for as an additional reduction factor f of the signal efficiency, which depends on the beam intensity. The uncertainty in this factor includes uncertainty due to the difference of the ECAL energy spectra for dimuon and A events which is taken into account by the reweighting procedure discussed below in Sec.VI.

C. The HCAL and ECAL energy distributions
An example illustrating good agreement between distributions of energy deposited by µ + µ − in the HCAL module 2, for the data and MC is shown in Fig. 9. On the right panel of the plot one can see a small peak at 2.5 GeV from single muons originated from events when one of the muon from the µ + µ − pair did not reach the HCAL3. An additional cross-check was made by comparing the distributions of the energy E ECAL deposited by scattered electrons from the reaction (16) in the ECAL taking into account small corrections due to dimuon energy depositions. This comparison of the data vs MC E ECAL distributions for the high intensity run III is shown in Fig. 10. One can see that the predicted and measured spectra are in a reasonable agreement and are not significantly distorted by pileup events.

VI. SIGNAL EFFICIENCY
Several signal detection efficiency contribute to the value of tot (m A ) in the NA64 detector: where e , A , ECAL , V and HCAL are the efficiency factors for the primary e − selection, which include also the reduction factor f discussed in Sec. V.B, the A acceptance in the signal box range, and the efficiencies for the signal to pass the ECAL, V 2 , and HCAL selection criteria, respectively. The ECAL value includes also the ECAL spectrum reweighting factor discussed below in Sec. VI.A. These factors were determined from the sample obtained with MC simulations and from the data samples of e − and dimuon events. The flux and spectra of the A s produced in the ECAL target by primary electrons were calculated using the approach reported in ref. [46] taking into account the development of the signal e-m shower from reaction (7) in the ECAL target (see, Sec. V).

A. The ECAL signal efficiency
The reconstruction efficiency ECAL for signal events was calculated for different A masses as a function of energy deposited in the ECAL. Compared to the ordinary em shower, the ECAL value for the e-m shower induced by an A event has to be corrected mainly due to difference in the longitudinal e-m showers development at the early stages in the PS detector [46]. This correction depends on the threshold E th P S of the energy deposited in the PS used in the trigger (14) and was typically (5 ± 3)% where the errors came from the E th P S threshold variation during data taking.
The sensitivity of the NA64 experiment is defined by the number of accumulated events which depends on the beam intensity. The intensity is limited by the pulse duration (τ ECAL 100 ns) from the ECAL MSADC shaper resulting in a maximally allowed electron counting rate of 10 6 e − /s in order to avoid significant loss of the signal efficiency due to the pileup effect. To evade this limitation, we have implemented a pileup removal algorithm to allow for high-efficiency reconstruction of the A signal and energy in high electron pileup environments, and run the experiment at the electron beam rate a few 10 6 e − /spill. This is in particular important in the case of signal events, because the shape of the E ECAL spectrum can be used for the A mass evaluation [46]. The shape is in particular sensitive to the mass in the low energy region which is the most affected by the pileup pulses which may occur somewhat earlier or later than the desired pulse and may seriously affect the reconstruction efficiency of signal events.
A simple pileup removal algorithm was used in the analysis of the data and MC samples of events obtained for high beam intensity. All ECAL cells were requested to have a single MSADS peak with a cell-time within ±2 ns of the trigger time if the energy deposited in the cell was more than 1 GeV. If several peaks were found, the one closer to the expected cell-time position was selected, and an attempt was made to remove contributions from the neighbouring pileup peak(s) to the signal area. The efficiency of the pileup removal algorithm as a function of the E ECAL value was estimated by using clean data and MC dimuon samples obtained at different intensities as described below. This method was also used to evaluate the efficiency of the A signal reconstruction in the energy range predicted by the simulations. At high in- tensity the dependency of dimuon events reconstruction efficiency on E ECAL can be important. For this reason using the efficiency corrections directly from the overall ratio data/MC in dimuons can be inaccurate due to the difference in E ECAL spectra. In order to check this the samples of reconstructed µ + µ − events from data and MC simulation were compared with more details. For this purpose the individual correction RR factors, ratios between efficiency measured in data and in Monte Carlo, have been derived and they were applied as an event weights in the MC simulations to obtain a better agreement between simulated and real data samples. These scale factors are used to correct Monte Carlo efficiencies to agree with measurements on data samples. The RR ratios have been formed as a function of E ECAL value: where n EC i and n EC tot is the number of events in i-th bin and the total number of events in the ECAL energy distribution, see Fig. 10. The ratio of the above ratios RR = R Data /R M C is then a measure of any additional differences in the signal reconstruction efficiency between data and MC as a function of the energy deposited in the ECAL. In Fig. 11 the distribution of RR values over the ECAL energy range from 1 to 60 GeV is shown for the beam intensity 5 × 10 6 e − /spill from run III. One can see that the method works well for the energy region, E ECAL 5 GeV, when the distortion of the ECAL spectrum from the pileup effect is relatively small at any distance between the true and pileup pulses in the ECAL yielding a correction factor close to 1. In the low energy region E ECAL 10 GeV, the reconstruction efficiency is more affected by the statistical uncertainties in the shape of the reconstructed pileup pulse, and the true pulse is not identified well anymore. This difference was used to additionally correct the MC efficiency for the signal in order to account for pileup and other effects not present in the MC. For this purpose the signal spectrum shown in Fig. 6 was reweighted by applying bin-by-bin corrections to the signal efficiency for the given energy obtained from dimuon data sample as shown in Figs.10 and 11. This procedure results in the overall correction factor to the signal efficiency 0.93 for the case of highest intensity of the run III, and is 95% for the runs I and II. It is slightly mass dependent.

B. The Veto and HCAL cuts and efficiency
corrections.
The V 2 and HCAL signal acceptance was defined as a fraction of events below the corresponding zero-energy cuts: where n HCAL(V2) (E < E th ), n tot is the number of events below the threshold energy E th and the total number of events, respectively. The shape of the distributions of energy deposited in these detectors from the leak of the signal shower energy, deposited in the ECAL, was simulated for different A masses and cross-checked with measurements at the e − beam of several energies. The HCAL energy cut was chosen by using the distribution of E HCAL in the data from the calibration runs, where the ECAL threshold was removed from the trigger. The admixture of non-electrons in these runs can be neglected. In Fig. 12 the right panel shows the dependence of HCAL on the energy threshold E th HCAL . One can see that for the value E th HCAL 1 GeV the acceptance is HCAL 0.98 for the maximal beam intensity.
The corrections to efficiency V and HCAL in the V 2 and HCAL were determined directly from the data using the calibration runs. The left panel in Fig. 12 shows an example of the measured distribution of the energy in the HCAL in such run. In general, the agreement with MC spectrum is good. The small differences between data and MC distributions are dominated by the pileup effects. The results for the V 2 acceptance look quite similar. Quantitatively, to estimate the V 2 and HCAL signal efficiency we studied the ratios of the number of events above and below the corresponding cuts used for the A event selection.
An example of the summary of corrected efficiencies used for calculations of the limits for the run III data sample obtained with a beam intensity 5 × 10 6 e − /spill is presented in Table II. These efficiencies were slightly different for the data samples from runs I and II, mostly because of the different pileup algorithm efficiencies which were rate dependent and also determined from measurements in calibration runs at different beam rates. The quantities e , ECAL were the most rate-dependent detection efficiencies of the tracker chambers and SRD, and ECAL cluster reconstruction, respectively. The DAQ deadtime was a function of the beam rate and was 7.4% averaged over the full data-taking period.
The total number of collected EOT in 2016 was ob-  tained from the recorded number of events from the em e − Z interactions in the ECAL target by taking into account the trigger suppression factor ( 10 2 ) and DAQ dead time which was beam rate dependent. The e − beam loss due to interactions with the beam line materials was estimated to be small. The trigger and SRD efficiency obtained by using unbiased samples of events that bypass the selection criteria was found to be 0.95 and 0.97 with a small uncertainty 2%. The probability of A events to pass all selection criteria, A was evaluated by processing the simulated signal events through the same reconstruction program as data, with the same cuts. The A yield calculated in accordance with Ref. [46,54] was corrected for the production cross section as described in Sec.II. The overall signal efficiency was in the range tot (0.7 − 0.5) decreasing for the higher intensity run.

VII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties are determined to stem from the overall normalization, signal cross section computations, reweighting the EAL signal energy distribution, description of the dimuon spectrum, and the uncertainty in the signal efficiencies in the Veto and HCAL. Systematic uncertainties are determined by varying cuts and taking the largest change in the calculated rate as the systematic error. Details of the systematic checks are given below.
The 10% additional uncertainty, estimated from the comparison of the cross sections calculated in [53] and [54] as discussed in Sec. II, was taken into account as the systematic error for the A production in the target. Note, that possible contributions from the purity of the target ( 99.9%) and 22% admixture of spin1/2 isotope 207 Pb are estimated to be small. Another contributions are due to the events selection in the ECAL and the reweighting procedure of the ECAL signal spectrum described in Sec. V. The former was estimated as a difference in the signal yield with respect to the nominal value due to the PS energy threshold variation during the run. This contribution increases for large A masses. The systematic uncertainty from the reweighting procedure was estimated by varying the parameters of the empirical fitting functions shown in Fig. 11 and considering differences in the number of obtained signal events. In this case the quoted systematic uncertainty is taken to be the quadratic sum of the observed shift and the statistical error of 4% on the shift. The reweighting correction is slightly A mass dependent and has the total uncertainty 7% for the highest intensity run III. Other contributions to the systematic uncertainty on the A signal efficiency come from the choice of the cut threshold and definition of the signal efficiency for the V 2 and HCAL, see Sec.VI. The uncertainties of the corrections for the V 2 and HCAL signal efficiency were studied by varying the corresponding energy threshold within the range determined from the data. The calculated variations were assigning to the systematic errors, which were estimated to be 3% and 2% for the V 2 and HCAL, respectively. An example of sources and the corresponding magnitudes of the systematic errors for the A masses 10 and 100 MeV, estimated for the run III is shown in Table III.

VIII. BACKGROUND
The search for the A → invisible decays requires particular attention to backgrounds, because every process with a single track and an e-m cluster in the ECAL can potentially fake the signal. In this Section we consider all background sources, which were also partly studied in Refs. [44,45].
There are several backgrounds resulting in the signature of Eq.(15) which can be classified as being due to detector-, physical-and beam-related sources. The selection cuts to reject these backgrounds have been chosen such that they do not affect the shape of the true E miss spectrum.
The estimation of background levels and the calculation of signal acceptance were both based on the MC simulation, as well as direct measurement with the beam. Because of the small A coupling strength value, performing a complete detector simulation in order to investigate these backgrounds down to the level of a single event sensitivity 10 −11 would require a very large amount of computing time. Consequently, we have estimated with MC simulations all known backgrounds to the extent that it is possible. Events from particle interactions or decays in the beam line, pileup activity created from them, hadron punchthrough from the target and the HCAL were included in the simulation of background events. Small event-number backgrounds such as the decays of the beam µ, π, K or µs from the reaction of dimuon production were simulated with the full statistics of the data. Large event-number processes, e.g. from e − interactions in the target or beam line, punchthrough of secondary hadrons were also studied, although simulated samples with statistics comparable to the data were not feasible. To eliminate possible instrumental effects not present in the MC calculations, the uniformity scan of the central part of the ECAL target was performed with e − by using the MM3 and MM4. We also examined the number of events observed in several regions around the signal box, which were statistically consistent with the estimates.
The main detector background sources are related to • Instrumental effects. The leak of energy throughout the possible holes, cracks, etc. in the downstream coverage of the detector which allows secondary particles to pass through without interactions. To study this effect a X − Y scan over the transverse area of the ECAL and HCAL detector has been performed with a particular attention to the boundaries between cells, fibers positions, and dead materials. No significant leak of energy has been observed.
• Detector hermeticity. The fake signature of Eq.(15) could also arise when either: i) a highenergy bremsstrahlung photon from the reaction eZ → eZγ, or ii) leading hadron h from the reaction eZ → eZX + h in the target escape detection due to punchthrough in the HCAL. The reaction i) may occur if an energetic photon induces a photonuclear reaction accompanied by the emission of a leading neutral particle(s), such as e.g. a neutron. The neutron then could be undetected in the rest of the detector. Taking into account the estimated non-hermeticity of the detector, the probability of the reaction is found to be 10 −14 . For the charged secondaries the punchthrough is highly suppressed by the observation of energy deposition in the HCAL modules. As the number of photoelectrons per MIP crossing the HCAL module was measured to be in the range n ph.e.
150 − 200/ MIP the inefficiency of the punchtrough detection is 10 −10 making the overall background negligible.
For the case ii) the punchthrough probability of a leading neutral hadron, such as a neutron and/or K 0 L , is defined by exp(−L tot /λ int ), where L tot is the (ECAL+HCAL) length sum. It has been estimated separately with a pion beam and compared with simulations [45] . It has been found that the overall hadron punchthrough probability is below 10 −12 for the total thickness of the ECAL and HCAL of about 30 λ int . This value should be multiplied by a factor 10 −4 , which is the probability of a leading hadron electroproduction in the ECAL target. Taking this into account the final estimate results to the negligible level of this background per incoming electron. The HCAL non-hermeticity for high energy neutral hadrons was cross-checked with Geant4-based MC simulations [46]. For the energy threshold E th HCAL 1 GeV the non-hermeticity is expected to be at the level 10 −9 . Taking into account the probability to produce a single leading hadron per incoming electron as P h 10 −4 , an overall level of this background of 10 −13 is obtained. This is in agreement with the above rough estimate.
• Large transverse fluctuations. Another possible source of background was caused by the large transverse fluctuations of hadronic showers from the reaction eZ → eZ+ ≥ 2 neutrals induced by electrons in the ECAL. In such events all secondary long-lived neutral particles (such as neutrons and/or K 0 L 's) could be produced in the target at a large angle, the HCAL and escape the detector without depositing energy through the lateral surface, thus resulting in the fake signal event. Taking into account results from the previous study [46,66,67], a conservative estimate for this background gives the level 10 −14 per incoming electron.
The beam backgrounds can be subdivided into two categories: upstream interactions and particle decays.
The main background sources is caused by the upstream beam interactions with beamline materials, such e.g. as entrance windows of the beam lines, residual gas, S1, MM1,2 etc., resulting in an admixture of low energy electrons with a large incoming angle in the beam. Those may fake a missing energy signal as they still could be within acceptance of the the spectrometer, while some or all of the accompanying produced secondaries fall outside the acceptance of the downstream ECAL and HCAL calorimeters. The limited detection acceptance of the secondaries along the downstream beam axis enhances this background. In 2017 run a zero-angle detector, as well as the lead-glass counters are planned to be installed to improve the downstream coverage and detection efficiency. The fraction of upstream scattered events is estimated to be at the level 10 −5 per EOT. An uncertainty arises also from the lack of accurate knowledge of the dead material composition in the beam line and is potentially the largest source of systematic uncertainty for accurate calculations of the fraction and energy distribution of these events. As it is not clear whether such rare large angle scattering could be reproduced with MC simulations the amount of background events from the beam upstream interactions was mainly estimated from the data itself.
In addition to the SRD cut which helps to re-ject low energy electrons and V 2 cut, which rejects most of the charged secondaries up to 35 cm away from the deflected beam axis, two additional cuts were used to study possible contribution from this source. The first one eliminated charged secondaries with multiple hits in the upstream tracker chambers more than expected from a single track event. The second one, used information on lateral reconstructed energy and time spread in the HCAL cells from charged and neutral secondaries. It was used to reject mostly events accompanied by low energy neutrals and charged secondaries with typically more activity in the HCAL than expected from interaction in the ECAL target. In Fig. 7 the comparison of events distribution before (central panels) and after using the HCAL cut (right panels) is shown. One can clearly see that the amount of background events in the vicinity of the masked signal region is substantially reduced. Finally the background level in the signal box was estimated from the extrapolation of the number of data events observed in dedicated control regions to the signal region using the fitting procedure described below. Namely, we looked at the ECAL energy distribution in the control region E ECAL > 50; E HCAL < 1 GeV for the runs I-III and estimated the contamination level in the signal box by fitting the E ECAL distribution with the where p1 is a constant, and p2 is the slope. In order to validate the exponential shape of the extrapolation function from control to signal region, dedicated validation region 1 < E ECAL < 80 GeV, E HCAL > 10 GeV containing a bigger data sample of events from hadronic interactions of the beam electrons with nuclei of the ECAL target was defined, see right panels in Fig. 7. The exponential shape of the distribution of energy deposited by the scattered electrons in the ECAL observed in the validation region was cross checked with MC simulations, see Fig. 5 in Ref. [46], and found to be in agreement for the full energy range. Therefore, the amount of background in the signal region was estimated form the extrapolation procedure assuming that the energy distribution of beam electrons scattered in hadronic reactions the upstream part of the beamline has exponential shape similar to the one observed in hadronic interactions of beam electrons in the ECAL target. Note that background of electrons from the upstream QED bremsstrahlung scattering is strongly suppressed as they typically follow the beam direction after the scattering, and fall outside of the detector acceptance after being deflected in the magnets. The yield of the background events was estimated by extrapolating the fit functions from the side band C to the signal box, see Fig. 7, assessing the systematic uncertainties by varying the back-ground fit functions within the corresponding errors. An example of the fit extrapolation for the side C of the ECAL energy distribution is shown in Fig. 13 for the run II. The slopes in the exponential fitting of the E ECAL distributions for the runs I-III were 0.315 ± 0.0021, 0.396 ± 0.0072, 0.49 ± 0.0026 in unit of 1/(GeV/c), and the number of n b events in the signal region were expected to be 0.043 ± 0.017, 0.041 ± 0.02, 0.01 ± 0.003, respectively. Possible variation of the HCAL zero-energy threshold during data taking were also taken into account. The fit was also performed for both sideband A shown in the right panel of Fig. 7. Events in the region A (E ECAL < 50 GeV ; E HCAL > 1 GeV ) are pure neutral hadronic secondaries produced by electrons in the ECAL target, while events from the region C (E ECAL 50 GeV ; E HCAL < 1 GeV ) are likely from the e − interactions in the upstream part of the beam line. As a result, 0.001 events in total are expected in the signal box from the side band A for all runs, and was further neglected.
Finally, the contamination of backgrounds to the signal region due to beam interactions was estimated to be 0.09 ± 0.03 events. The uncertainty in the background estimate due to upstream scattered events was dominated by the systematic uncertainty of the upstream veto V1 and tracker efficiency, precise knowledge of the material in the beam line, and statistics of the data samples. This systematic uncertainty was estimated by performing measurements on several samples of upstream scattered events tagged by a signature of scattering in the HCAL. The uncertainties in this background estimate are evaluated by considering differences in the estimates of the event number by varying the electron identification probabilities and changing the parameters of the extrapolation functions.
• Particle decays. Other backgrounds were expected from the decays of µ, π, K → e + in flight in the beamline accompanied by emission of an energetic neutrino. These backgrounds were highly suppressed by requiring the presence of the SRD tag. However, there might be cases when, e.g. a pion could knock electrons off the downstream window of the vacuum vessel, which hit the SRD creating a fake tag for a 100 GeV e − . The pion could then decay into eν in the upstream ECAL region thus producing the fake signal.
The main background source in this category was K e3 decays where the electron overlapped with photons from π 0 decay thus producing a single-like e-m shower in the ECAL. In addition to the SRD cut and the probability of decay in the downstream part of the setup, this process was further suppressed by requiring shower energy to be < 50 GeV and the incoming track azimuthal angle to be below 5 mrad. The transverse and longitudinal shower Background source Estimated number of events, n b hermeticity: punchthrough γ's, cracks, holes < 0.001 loss of hadrons from e − Z → e − + hadrons < 0.001 loss of muons from e − Z → e − Zγ; γ → µ + µ − 0.005 ± 0.001 µ → eνν, π, K → eν, Ke3 decays 0.02 ± 0.004 e − interactions in the beam line materials 0.09 ± 0.03 µ, π, K interactions in the target 0.008 ± 0.002 accidental SR tag and e − from µ, π, K decays < 0.001 Total n b 0.12 ± 0.04 shape at the ECAL was also used to distinguish the single electron shower from the overlapped one.
Similar background was caused by a random superposition of uncorrelated low-energy, 50 -70 GeV, electron from the low-energy beam tail and 100 GeV beam µ, π, K occurring during the detector gate-time. The electron could emit the amount of SR energy above the threshold which is detected in the SRD as a tag of 100 GeV e − and then is deflected by the spectrometer magnets out of the detector's acceptance angle. While the accompanying mistakenly tagged µ, π or K could either decays in-flight in front of the ECAL into the e − + X state with the decay electron energy less then the beam energy, or could also interact in the target producing an e-m like cluster below 50 GeV though the µZ → µZγ or π, K charge-exchange reactions, accompanied by the poorly detected scattered µ, or secondary hadrons, thus resulting in both cases to the signal signature of Eq. (15). These background components were simulated with a statistics higher or comparable to the number of events expected from the data and was found to be small.
The remaining physical backgrounds were • Dimuon, τ , charm decays. The process (16) could mimic the signal either i) due to muons decay in flight inside the ECAL target into eνν state, or ii) if the muons escape detection in the V 2 and HCAL modules due to fluctuations of the energy (number of photoelectrons) deposited in these detectors. In the case i) the relatively long muon lifetime results in a small probability to decay inside the ECAL. For the case ii) the background is suppressed by the high-efficiency veto system V 2 +HCAL. The V 2 was a ∼ 4 cm thick high-sensitivity scintillator array whose inefficiency for a single muon detection was estimated to be 10 −4 . Therefore, the level of dimuon background is expected to be < 10 −13 per EOT. The fake signal could also arise from the reactions of τ , e.g., eZ → eZτ + τ − ; τ → eνν, or charm, e.g., eZ → eZ + D s + anything; D S → e + ν + anything, production and their subsequent prompt decays into an electron accompanied by emission of neutrinos. The estimate show that these backgrounds are also expected to be negligible.
• Finally, the electroproduction of a neutrino pair eZ → eZνν resulting in the invisible final state accompanied by energy deposition in the ECAL1 from the recoil electron can occur. An estimate showed that the ratio of the cross sections for this reaction to the bremsstrahlung cross section is well below 10 −13 [44].
In Table IV the contributions from all background processes estimated by using the MC simulations, exept for those from beam interactions in the upstream part of the setup, are summarized. The final number of background events estimated from the combined MC and data events is n b = 0.12 ± 0.04 events for 4.3 × 10 10 EOT. The estimated uncertainty of about 30% was due mostly to the uncertainty in background level from upstream beam interactions. It also includes the uncertainties in the amount of passive material for e − interactions, in the cross sections of the hadron charge-exchange reactions on lead (30%), and systematic errors related to the extrapolation procedure. The total systematic uncertainty was calculated by adding all errors in quadrature.

IX. RESULTS AND CALCULATION OF LIMITS
In the final statistical analysis the three runs I-III were analysed simultaneously using the multi-bin limit setting technique. The corresponding code is based on the RooStats package [68]. First of all, the above obtained background estimates, efficiencies, and their corrections and uncertainties were used to optimize more accurately the main cut defining the signal box by comparing sensitivities, defined as an average expected limit calculated using the profile likelihood method, with uncertainties used as nuisance parameters. Log-normal distribution was assumed for the nuisance parameters [69]. The most important inputs for this optimization were the expected values from the background extrapolation into the signal box for the data samples of the runs I,II,III. The uncertainties for background prediction were estimated by FIG. 14: The sensitivity, defined as an average expected limit, as a function of the ECAL energy cut for the case of the A detection with the mass m A 20 (blue) and 2 (green) MeV.
varying the extrapolation functions, as previously discussed. An example of the optimization curves obtained for the m A = 2 and 20 MeV is shown in Fig. 14. It was found that the optimal cut value depends very weakly on the A mass choice and can be safely set to E ECAL < 50 GeV for the whole mass range.
Overall optimization and improvement of the signal selection and background rejection criteria resulted in roughly more than a factor 10 reduction of the expected backgrounds per EOT and an increase of a factor 2 in the efficiency of A → invisible decay at higher beam rate for the run III compared to those obtained in the analysis reported in Ref. [38]. For the full 2016 exposure, the estimate of the number of background events expected from the sources discussed above per 10 10 EOT was n b = 0.03, while for the run of Ref. [38] it was n b = 0.5.  . For more limits obtained from indirect searches and planned measurements see e.g. Ref. [13,14].
After determining and optimizing all the selection criteria and estimating background levels, we examined the events in the signal box and found no candidates, as shown in Fig. 7. We proceeded then with the calculation of the upper limits on the A production. The combined 90% confidence level (C.L.) upper limits for the corresponding mixing strength were determined from the 90% C.L. upper limit for the expected number of signal events, N 90% A by using the modified frequentist approach for confidence levels (C.L.), taking the profile likelihood as a test statistic in the asymptotic approximation [70][71][72]. The total number of expected signal events in the signal box was the sum of expected events from the three runs: where i tot is the signal efficiency in the run i given by Eq. (23), and the n i A ( , m A , ∆E A ) value is the signal yield per EOT generated by a single 100 GeV electron in the ECAL target in the energy range ∆E e . Each ith entry in this sum was calculated by simulating the signal events for corresponding beam running conditions and processing them through the reconstruction program with the same selection criteria and efficiency corrections as for the data sample from the run-i. The expected backgrounds and estimated systematic errors were also taking into account in the limits calculation. The combined 90% C.L. exclusion limits on the mixing strength as a function of the A mass can be seen in Fig. 15. In Table V [13,[25][26][27] from the results of the LSND [24,36], E137 [37], BaBar [39], MiniBooNE [43] and direct detection [76]experiments. The favoured parameters to account for the observed relic DM density for the scalar, pseudo-Dirac and Majorana type of light thermal DM are shown as the lowest solid line.
different m A values are also shown for comparison. One can see that the corrections are mostly relevant in the higher mass region m A 100 MeV. The derived bounds are the best for the mass range 0.001 m A 0.1 GeV obtained from direct searches of A → invisible decays [15].
The limits were also calculated with a simplified method by merging all three runs into a single run as described previously by Eq. (26). The total error for the each N i A value includes the corresponding systematic uncertainties calculated by adding contributions from all sources in quadrature, see Sec.VII. In accordance with the CL s method [72], for zero number of observed events the 90% C.L. upper limit for the number of signal events is N 90% A (m A ) = 2.3. Taking this and Eq.(26) into account and using the relation N A (m A ) < N 90% A (m A ) resulted in the 90% C.L. limits in the (m A ; ) plane which agreed with the one shown in Fig. 15 within a few %.

X. CONSTRAINTS ON LIGHT THERMAL DARK MATTER
As discussed previously, the possibility of the existence of light thermal Dark Matter (LTDM) has been the subject of intense theoretical activity over the past several years [13,14], see also [73,74]. The LTDM models can be classified by the spins and masses of the DM particles and mediators. The scalar dark matter mediator models are severely restricted or even excluded by non-observation of rare B-meson decays [14,15], so we consider here only the case of a vector mediator. As was discussed in Sec.I, the most popular vector mediator model is the one with additional massive dark photon A which couples with DM particles via interaction L = e D A µ J µ χ . The currents J µ χ =ψ χ γ µ ψ χ and J µ χ = i(φ + χ ∂ µ φ χ − φ χ ∂ µ φ + χ ) for spin 1/2 and 0, respectively. Here, χ denotes both, either scalar or fermion LTDM particle. As discussed in Sec.I, the γ − A mixing leads to nonzero interaction of dark photon A µ with the electrically charged SM particles with the charges e = e. As a result of the mixing the cross-section of DM particle annihilation into SM particles, which determines the relic DM density, is proportional to 2 . Hence using constraints on the cross section of the DM annihilation freeze out (resulting in Eq.(5)), and obtained limits on mixing strength of Fig. 15, one can derive constraints in the (y;m χ ) plane, which can also be used to restrict models predicting existence of LTDM for the masses m χ 1 GeV.
These limits obtained from the full data sample of the 2016 run are shown in the left panel of Fig. 16 together with the favoured parameters for scalar, pseudo-Dirac (with a small splitting) and Majorana scenario of LTDM taking into account the observed relic DM density [14]. The limits are calculated by using Eq.(6) under the conventional assumption α D = 0.5, and m A = 3m χ , here m χ stands for the LTDM particle's masses, either scalars or fermions. The plot shows also the comparison of our results with limits from other experiments. Note, that some of these limits were obtained by using WW approximation for the cross section calculation and therefore might require revision. The choice of α D = 0.5 is compatible with the bounds derived in Ref. [75] based on the running of the dark gauge coupling. However, it should be noted that differently form the results of beam dump experiments, such as LSND [24,36], E137 [37], MiniBooNE [43], the χ-yield in our case scales as 2 , not as 4 α D . Therefore, for sufficiently small values of α D our limits will be much stronger. This is illustrated in the right panel of Fig. 16, where the NA64 limits and bounds from other experiments are shown for α D = 0.005. One can see, that for this, or smaller, values of α D , the direct  [14] from the results of the LSND [24,36], E137 [37], BaBar [39] and MiniBooNE [43] experiments.
search for the A → invisible decay in NA64 excludes model of scalar and Majorana DM production via vector mediator for the remaining mass region m χ 0.05 GeV. While being combined with the BaBar limit [39], the result excludes the model for the entire mass region m χ 1 GeV.
The experimental upper bounds on also allow to obtain lower bounds on coupling constant α D which are shown in Fig. 17 in the (α D ;m χ ) plane. For the case of pseudo-Dirac fermions and small splitting, the limits in the left panel of Fig. 17 were calculated by taking the value f = 0.25 in Eq. (5). For the mass range m χ 0.05 GeV the obtained bounds are more stringent than the limits obtained from the results of LSND [24,36] and E137 [37]. The limits for the Majorana case shown in the right panel of Fig. 17 were calculated by setting f = 3. To cross check our calculations, we also derived limits on α D by using BaBar bounds on [39], see Fig. 15, Eq. (5) and the previous f values for the pseudo-Dirac and Majorana cases. The obtained BaBar limits were found to be in good agreement with those shown in Fig. 17 for the mass region m χ 0.1 GeV. Note, that new constraints for the large pseudo-Dirac fermion splitting can also be derived. They will be more stringent than for the case of the small splitting and similar to the one obtained for the Majorana case.

XI. CONCLUSION
From the analysis of the full 2016 data sample, we found no evidence for the existence of dark photon with the mass in the range 1 GeV which mixes with the ordinary photon and decays dominantly invisibly into light DM particles A → χχ. New limits on the mix-ing strength were derived by taking into account the A production cross sections calculated at the exact treelevel in Ref. [54] without using Weizsäcker-Williams approximation. These cross sections, implemented into NA64 simulation package, were cross checked with the one calculated by Liu et el. following their approach reported in Ref. [52,53] and was found to be in agreement. Good agreement between the data and MC for the rare QED dimuon production in the reaction (16) was also found. This process was used as a robust benchmark for the signal event simulation and analysis. For the mass range m A 0.1 GeV the most stringent upper limit on the mixing strength, were obtained. Using conventional choices for the DM parameters we also set the 90% C.L. limits on the value y, representing the dimensionless DM annihilation cross section parameter, for the pseudo-Dirac and Majorana DM in the χ mass region 0.001 < m χ < 0.1 GeV. With these DM parameter combinations, our result has expanded the search for DM to y values an order of magnitude smaller than MiniBooNE DM experiment [43].
For the vector portal DM model and the chosen parameter constraints, the obtained lower limits α D 10 −3 for pseudo-Dirac Dark Matter in the mass region m χ 0.05 GeV are more stringent than the bounds from beam dump experiments. For values α D 0.005 the combined results from direct searches of the A → invisible decay in NA64 and BaBar experiments exclude model for scalar and Majorana DM production via vector portal for the mass region m χ 1 GeV. The obtained results are used to constrain an interpretation of the A as a mediator of light thermal DM production. The remaining windows of parameter space between our result and the Landau pole bounds obtained from arguments based on the running of the dark gauge coupling [75] for these scenarios can be covered through future searches with the NA64 experiment. This would require an improved sensitivity of about two orders of magnitude which is feasible.