Constraints on the kinetic mixing parameter $\epsilon^2$ for the light dark photons from dilepton production in heavy-ion collisions in the few-GeV energy range

The vector $U$-bosons, or so called 'dark photons', are one of the possible candidates for the dark matter mediators. They are supposed to interact with the standard matter via a 'vector portal' due to the $U(1)-U(1)^\prime$ symmetry group mixing which might make them visible in particle and heavy-ion experiments. While there is no confirmed observation of dark photons, the detailed analysis of different experimental data allows to estimate the upper limit for the kinetic mixing parameter $\epsilon^2$ depending on the mass $M_U$ of $U$-bosons which is also unknown. In this study we present theoretical constraints on the upper limit of $\epsilon^2(M_U)$ in the mass range $M_U \le 0.6$ GeV from the comparison of the calculated dilepton spectra with the experimental data from the HADES Collaboration at SIS18 energies where the dark photons are not observed. Our analysis is based on the microscopic Parton-Hadron-String Dynamics (PHSD) transport approach which reproduces well the measured dilepton spectra in $p+p$, $p+A$ and $A+A$ collisions. Additionally to the different dilepton channels originating from interactions and decays of ordinary matter particles (mesons and baryons), we incorporate the decay of hypothetical $U$-bosons to dileptons, $U\to e^+e^-$, where the $U$-bosons themselves are produced by the Dalitz decay of pions $\pi^0\to \gamma U$, $\eta$-mesons $\eta \to \gamma U$ and Delta resonances $\Delta \to N U$. Our analysis can help to estimate the requested accuracy for future experimental searches of 'light' dark photons by dilepton experiments.


I. INTRODUCTION
An understanding of the structure of our Universe is one of the intriguing topics of modern physics. According to the present knowledge, the standard matter represents less then 5% of our Universe, while about 27% of it consists of so-called 'dark matter' (DM) and about 68% is the 'dark energy' [1]. The dark matter is supposed to be a relic from the Big Bang, which makes itself noticeable by its gravitational action on the large-scale cosmic structures. It was advocated that the dark matter mediators can interact with the Standard Model (SM) particles by four possible 'portals' -vector, Higgs, neutrino and axion (cf. the reviews [2][3][4] and references therein).
The 'vector' portal assumes the existence of a U (1) − U (1) gauge symmetry group mixing [5], i.e. the corresponding Lagrangian is defined by the hypercharge fieldstrength tensor of the SM photon field and the DM vector boson field: L ∼ 2 /2 F µν F µν . The mediators in this case are vector U -bosons, which often are called 'dark' or 'hidden' photons or A , with a mass M U remaining presently unknown. Here 2 is a kinetic mixing parameter, which characterizes the strength of the interaction of SM and DM particles [6][7][8][9][10][11]. This mixing allows the decay of U -bosons to a pair of leptons -e + e − or µ + µ − . The 'light' U -bosons can be produced by the decay of SM particles, e.g. by Dalitz decays of pseudoscalar mesons -pions and η-mesons, as well as by the Dalitz decay of * Electronic address: E.Bratkovskaya@gsi.de baryonic resonances such as ∆'s. This provides the possibility to observe dark photons in dilepton experiments, stimulating a lot of experimental as well as theoretical activities [2,3,12,13]. A recent measurement of the excess electronic recoil events by the XENON1T Collaboration might be interpreted also in favour of dark matter sources and in particular dark photons are possible candidates [14].
The HADES Collaboration at GSI, Darmstadt, performed an experimental search for dark photons in dilepton experiments at SIS18 1 with both proton and heavyion beams [15]. The HADES experiment presented an upper limit for the kinetic mixing parameter 2 in the mass range of M U = 0.02 − 0.55 GeV based on the experimental measurements of e + e − pairs from p + p and p+N b collisions at 3.5 GeV as well as Ar+KCl collisions at 1.76 AGeV. Later the HADES result has been superseded by the A1 [16], the NA48/2 [17] and the BaBar [18,19] experiments which further lowered the limit for 2 in this mass range. The NA48/2 experiment investigated a large sample of π 0 Dalitz decays obtained from in-flight weak kaon decays, the BaBar collider experiment used their cumulated luminosity of e + e − reactions to survey a very large mass range up to M U = 8 GeV, and the A1 experiment [16] at MAMI investigated electron scattering off a 181 T a target at energies between 180 and 855 MeV to search for a dark photon signal. In the mass range discussed here, i.e. 20 -500 MeV, the limit on 2 has thus been pushed down to about 10 −6 . The compilation of the world data collected by various experiments can be found in the review [3]. The goal of this study is to estimate the upper limit for the kinetic mixing parameter 2 (M U ) depending on the mass of the hypothetical U -boson from the theoretically calculated dilepton spectra using the microscopic Parton-Hadron-String Dynamics (PHSD) transport approach which describes the whole evolution of heavy-ion collisions based on microscopic transport theory by solving the equations-of-motion for each degree-of-freedom (hadronic and partonic) and their interactions. The PHSD provides a consistent description of hadron production in p+p, p+A and A+A collisions as well as electromagnetic probes -dileptons and photons, from SIS18 to LHC energies (cf. [20][21][22]). In particular, the PHSD describes very well the dilepton data of the HADES experiment [23]. Having the SM particle production under control, we incorporated in the PHSD the dynamical production of U -bosons by π 0 , η and ∆ Dalitz decays as well as the U -boson decay to e + e − pairs. We compare our results with the experimental data from the HADES Collaboration and provide constraints on the mass dependence of the kinetic mixing parameter 2 for light dark photons (M U ≤ 0.6 GeV).
Our paper is organised as follows: In Section II we recall the basic ideas of the PHSD approach and the SM production. In Section III we describe the production of dark photons and in Section IV we show the results for the mixing parameter 2 (M U ). We summarise our findings in Section V.

II. STANDARD MATTER (SM) PRODUCTION IN THE PHSD
In this section we recall the basic ideas of the PHSD approach and treatment of the dilepton production from SM matter.

A. The PHSD transport approach
The Parton-Hadron-String Dynamics is a nonequilibrium microscopic transport approach that incorporates hadronic as well as partonic degrees-of-freedom [20,[24][25][26][27][28][29]. The PHSD describes the full evolution of a relativistic heavy-ion collision, from the initial hard NN collisions out-of equilibrium, to the formation of the quark-gluon plasma (QGP), its partonic interactions, up to the hadronisation and final-state interactions of the resulting hadrons. The dynamical description of the time evolution of the interacting system is based on the solution of the Cassing-Juchem generalised off-shell transport equations in test-particle representation [30,31] on the basis of the Kadanoff-Baym equations [32] in first-order gradient expansion, which are applicable for the dynamical description of strongly interacting degrees-of-freedom [25,33].
The hadronic sector follows the early developments from the HSD transport approach [34,35] which includes as explicit hadronic degrees-of-freedom the baryon octet and decouplet, the 0 − and 1 − meson nonets and higher resonances. The description of multi-particle production in elementary baryon-baryon (BB), meson-baryon (mB) and meson-meson (mm) reactions is incorporated based on the Lund model [36]. This is realized in terms of the event generators FRITIOF 7.02 [36,37] and PYTHIA 6.4 [38] which are "tuned" (cf. Ref. [39] for details) for a better description of elementary reactions at low and intermediate energies as well as for the incorporation of the in-medium effects related to chiral symmetry restoration and modification of the properties of the formed hadronic degrees-of-freedom, while the string formation and decay occurs in a dense hadronic medium [40][41][42][43][44]. The strings are built from primary nucleonnucleon (N N ) collisions and secondary energetic BB, mB and mm interactions (the corresponding thresholds are taken as (s th BB ) 1/2 = 2.65 GeV, (s th mB ) 1/2 = 2.4 GeV and (s th mm ) 1/2 = 1.3 GeV). They decay to the leading hadrons (the energetic ends of the strings) and 'prehadrons', i.e. the newly produced mesons and baryons which are considered under formation time and v is the velocity of the particle in the calculational frame which is chosen to be the initial N N center-of-mass frame).
The transition from hadronic to partonic degrees-offreedon and vice versa (i.e. hadronization) occurs when the local energy density is above (below) the critical energy density of C ∼ 0.5 GeV/f m 3 in line with lattice quantum chromodynamics (lQCD) [45]. If the energy density is below critical the 'pre-hadrons' evolve into asymptotic hadronic states after the formation time t F and interact with hadronic cross sections.
The description of the partonic degrees-of-freedom and their interactions during the QGP phase is based on the Dynamical Quasi-Particle Model (DQPM) [46,47] which describes the thermodynamic properties of QCD in equilibrium in terms of massive strongly-interacting quasiparticles whose masses are distributed according to spectral functions (imaginary parts of the complex propagators). The widths and pole masses of the spectral functions are defined by the real and imaginary parts of the parton self-energies and the effective coupling strength in the DQPM; both depend on the local temperature T (or local energy density) and the baryon chemical potential µ B . They are fixed by adjusting the DQPM entropy density to the respective lQCD results from Refs. [48,49]. The QGP phase is evolved by the off-shell transport equations with self-energies and cross sections from the DQPM. When the fireball expands the probability of the partons for hadronization increases close to the phase boundary between hadronic and partonic phases (which is a crossover at all RHIC energies); the hadroni-sation takes place using covariant transition rates. The resulting hadronic system is further-on governed by the off-shell HSD dynamics incorporating (optionally) selfenergies for the hadronic degrees-of-freedom [43].
Since in this study we concentrate on low energy heavyion collisions, only hadronic sources are relevant for the considerations here. For the description of dilepton production from the QGP sources we refer the reader to the review [20] and to the recent study on the relative contribution of the QGP dileptons and those from the correlated charm decays [22].
The dilepton production by a (baryonic or mesonic) resonance R decay can be schematically presented in the following way: i.e. in a first step a resonance R might be produced in baryon-baryon (1), meson-baryon (2) or meson-meson (3) collisions (or at high energy collisions the resonance can be formed in the hadronization process). Then this resonance can decay to dileptons directly, e.g. a direct decay of vector mesons (ρ, ω, φ). The decay of resonances to dileptons can be accompanied by the production of other particles (5) (e.g., Dalitz decay of the ∆ resonance: ∆ → e + e − N ). It can decay firstly to a meson m (+ hadrons) and later to dileptons (6) by e.g. the Dalitz decay (π 0 , η, ω). The resonance R might also decay into another resonance R (7) which later produces dileptons via Dalitz decay. The electromagnetic part of all conventional dilepton sources -π 0 , η, ω, ∆ Dalitz decays as well as direct decay of vector mesons ρ, ω and φ -are calculated following our early work [51]. However, here we adopt the "Wolf model" for the differential electromagnetic decay width of the ∆ resonance [52] which is a default setting in the PHSD 4.0 used for this study. The details of the evaluation of the ∆ Dalitz decays are given in Ref. [23].
For the bremsstrahlung from pp and pn reactions (N N → N N γ * → N N e + e − ) as well as from πN 'quasielastic' scattering (πN → πN γ * → πN e + e − ) we adopt the results from the OBE model calculations by Kaptari and Kämpfer in Ref. [53] as implemented in the PHSD in Ref. [42] and used for the dilepton study at SIS18 energies in Ref. [23].
We note that we account for the in-medium effects for the vector meson dynamics such as a 'collisional broadening' scenario for the spectral functions [42]. The vector meson production and propagation follow the off-shell transport description [33,47] where their spectral functions change dynamically during the propagation through the medium and evolve towards on-shell spectral functions in the vacuum. The effect of collisional broadening of the vector-meson spectral functions is incorporated by a modification of the vector meson width in the dense baryonic medium: Here Γ V (M ) is the total width of the vector mesons (V = ρ, ω, φ) in the vacuum. The collisional width in (8) is approximated as Here v = | p|/E; p, E are the velocity, 3-momentum and energy of the vector meson in the rest frame of the nucleon current and γ 2 = 1/(1 − v 2 ); ρ N is the nuclear density and σ tot V N the meson-nucleon total cross section. We use the 'broadening coefficients' α coll ≈ 150 MeV for the ρ and α coll ≈ 70 MeV for ω mesons as obtained in Ref. [42]. We note that the collisional broadening scenario is supported by a comparison of transport model calculations with dilepton data from low energies [23,42,[54][55][56] to ultra-relativistic energies (cf. e.g. [20,[57][58][59]).
The dilepton yield is calculated using the time integration (or 'shining') method, i.e. the virtual emission of dileptons by different sources is accommodated in the dilepton rates at each time step. In this way, the vector mesons and baryon resonances continously 'emit' dileptons from their production ('birth') up to their absorption or hadronic decay ('death'). This is especially important for the study of in-medium effects because this method takes the full in-medium dynamics into account.

III. U-BOSON PRODUCTION AND THEIR DILEPTON DECAY IN THE PHSD
Now we step to the description of e + e − pair production from U -boson decays U → e + e − in the PHSD approach. We note that at SIS18 energies, which are in the focus of our present study, the dominant production channels of dileptons for M < 0.6 GeV are the Dalitz decays of π 0 , η, and the ∆ resonance [23,42]. Thus, one expects that if the hypothetical dark photons have a mass M U < 0.6 GeV, they might stem from the Dalitz decays of π 0 , η, and the ∆ resonance, too. The evaluation of the corresponding partial decay widths of pseudoscalar mesons and ∆ baryons to U -bosons due to the U (1) − U (1) mixing has been performed in Refs. [10,11] and employed also in the HADES experimental search for dark photons in Ref. [15].
In this study we follow the strategy of Ref. [15] and consider three production channels of dark photons which are dominant for low energy heavy-ion collisions as measured by the HADES Collaboration [15], i.e. dark photons from Dalitz decay of 1) neutral pions: π 0 → γ+U , 2) and η-mesons: The dilepton yield from a U -boson decay of mass M U can be evaluated as the sum of all possible contributions (for a given mass M U ): where Br U →e + e − is the branching ratio for the decay of U -bosons to e + e − . In this study we assume that the width of the U -boson is zero (or very small), i.e. it contributes only to a single dM bin of dilepton spectra from SM sources. Accordingly, if M U > m η , only the ∆ channel is kinematically possible in Eq. (10). On the other hand, the yield of U -bosons of mass M U themselves can be estimated from the coupling to π 0 , η and ∆ decays to the virtual photons [15]: Following Refs. [10,11] the ratio of the partial widths for the Dalitz decays of π 0 , η mesons to U -bosons and real photons can be evaluated as follows: for m = π 0 , η. Here 2 is the kinetic mixing parameter and λ is the triangle function (λ(x, y, z) = (x − y − z) 2 − 4yz) from the expression of particle 3-momentum. Since m γ = 0, one obtains: Here, M U is the U -boson mass, F m are the electromagnetic transition formfactors for π 0 and η; they are taken as in our previous studies [23,42] and in the experimental HADES study [15] as well: with Λ = 0.72 GeV.
The U -boson production by the ∆ Dalitz decay ∆ → N U has been proposed in Ref. [10]. For the evaluation of the partial decay widths of a broad ∆ resonance, one has to take into account the ∆ spectral function A(M ∆ ) as used also in the HADES study [15]: where M ∆ is the mass of the ∆ resonance distributed according to the spectral function A(M ∆ ), m N the mass of the remaining nucleon. Following Ref. [15] we adopted |F ∆ (M 2 U )| = 1 for this study since an experimental formfactor is unknown.
In the PHSD the spectral function of a ∆ resonance of mass M ∆ is taken in the relativistic Breit-Wigner form [23]: with M ∆0 being the pole mass of the ∆. The factor C 1 is fixed by the normalization condition: where M lim = 2 GeV is chosen as an upper limit for the numerical integration. The lower limit for the vacuum spectral function corresponds to the nucleon-pion decay, M min = m π + m N . In N N collisions the ∆'s can be populated up to the M max = √ s − m N and hence the available part of the spectral function depends on the collision energy. We recall that for the total decay width of the ∆ resonance Γ tot ∆ (M ∆ ) in the PHSD we adopt the "Monitz model" [60] (cf. also Ref. [52]): We note that when accounting for the mass-dependent total width of the ∆ resonance our calculation for Uboson production by the ∆ Dalitz decay differs from the evaluation in Ref. [15] where a constant total width of the ∆ has been used. As discussed in Ref. [23] (see Section VI), the shape of the spectral function strongly depends on the actual form of Γ tot ∆ (M ∆ ). The branching ratio for the decay of U -bosons to e + e − , entering Eq. (10), is adopted from Ref. [11] and used also in Ref. [15]: = 1 .
Here m µ is the muon mass. The total decay width of a Uboson is the sum of the partial decay widths to hadrons, e + e − and µ + µ − pairs: Γ U tot = Γ U →hadr + Γ U →e + e − + Γ U →µ + µ − . The expression (22) has been evaluated using that Γ U →µ + µ − = Γ U →e + e − due to lepton universality for M U 2m µ . The hadronic decay widths of U -bosons is chosen such that Γ U →hadr = R( √ s = M U )Γ U →µ + µ − , where the factor R( √ s) = σ e + e − →hadrons /σ e + e − →µ + µ − is taken from Ref. [61]. In this Section we present the numerical results within the PHSD approach for the dilepton spectra including the contribution from the dilepton decay of U -bosons.
Since the kinetic mixing parameter 2 is unknown as well as the mass of the U -boson, we invent the following procedure to obtain the constraints on 2 (M U ): for each bin in dilepton mass dM , which is taken to be 10 MeV in our simulations, we calculate the integrated yield of dileptons from U -bosons of masses [M U , M U + dM ] according to Eq. (10) and divide by the bin size dM . The resulting dilepton yield per bin dM we denote as dN sumU /dM , which is the sum of all contributions (kinematically possible in the mass bin) from the dilepton decay of U -bosons produced by the Dalitz decays of π 0 , η and the ∆ resonance. Assuming that 2 is a constant in dM we can write that dN sumU /dM = 2 dN sumU =1 /dM where the notation dN sumU =1 /dM is the dilepton yield calculated without 2 or formally with = 1.
Thus, the total sum of all possible sources of dileptons, from the SM channels and from U -boson decays, can be written as Now we can obtain constraints on 2 (M U ) by requesting that the total sum dN total /dM cannot surplus the sum of SM channels by more than a fraction C U in each bin dM , i.e. C U controls the additionally "allowed" dielectron yield resulting from dark photons on top of the total SM yield (e.g. C U = 0.1 indicating that the dark photons add 10% extra yield to the SM yield, C U = 0.2 meaning 20% extra, etc.). We then express this as Combining Eqs. (24) and (23), one obtains that the kinetic mixing parameter 2 for M U can be evaluated as Eq. (25) allows to compute 2 for each bin [M U , M U + dM ] and presents the properly weighted dilepton yield from dark photons relative to the SM contributions. Moreover, now we can apply the experimental acceptance for e + e − pairs from U -boson decays in the same way as for the SM channels and compare our results to the experimental data from the HADES Collaboration. The latter will allow to explore the possible range of the factor C U that controls the additional yield from dark photons to the SM contributions. Since the dark photons have not been observed in any dilepton experiments, one can require that this enhancement should be still in the acceptable agreement with experimental data, i.e. within the experimental error bars (under condition that the SM yield agrees well with experimental data).
In Fig. 1 we present the compilation of the PHSD results for the differential cross section dσ/dM for e + e − production in p + p (upper, left) and p + N b reactions (upper, right) at 3.5 GeV beam energy and for the mass differential dilepton spectra dN/dM -normalized to the π 0 multiplicity -for Ar + KCl collisions at 1.76 A GeV (lower, left) and for Au + Au collisions at 1.23 A GeV (lower, right) in comparison to the experimental measurements by the HADES Collaboration. The solid dots present the HADES data for p + p [62], for p + N b [63,64] for Ar + KCl [65] and for Au + Au [66]. We note that the theoretical calculations are passed through the corresponding HADES acceptance filter and mass/momentum resolution.
The contributions from the various SM channels of dilepton production in the PHSD calculations are presented as individual colored lines described in the legend: the black dashed lines show the contribution from π 0 Dalitz decays, the red dot-dashed lines -from η Dalitz decays, the blue dot-dot-dashed lines -from ∆ Dalitz decays, the green short dashed lines -from ω Dalitz decays, the light blue short dashed lines -from direct decay ω → e + e − , the magenta dotted lines -from direct decay ρ → e + e − , the grey dotted lines with open dots -from direct decay φ → e + e − , the navy solid lines correspond to the sum of all SM contributions ("Sum SM" in legend). We note that for low beam energies (E kin < 2 A GeV) we account also for the Bremsstrahlung contributions (as described in Section II.B): the dark cyan dashed lines show the N N = pn+pp Bremsstrahlung and orange wave lines indicate the πN Bremsstrahlung. However, we do not include the Bremmsstrahlung contributions for larger energies since the OBE calculations used for low energies can not be easily extrapolated to the high energy regime. As has been shown in Ref. [23], such an extrapolation leads to a slight overestimation of the experimental data. We stress that we include the collisional broadening scenario for the vector meson (ρ, ω, φ) spectral functions which leads to a smearing of the peaks in the case of heavyion collisions and to an approximate exponential shape of the dilepton yields for heavy systems such as Au+Au collisions compared to the clear peak structure for p + p collisions. As seen from Fig. 1, the PHSD results for the SM sources are in a good agreement with the HADES results for all four systems as well as with the previous PHSD results in Ref. [23].
The individual contributions from U -boson decay to dileptons are shown in Fig. 1 as symbols: the contributions from U → e + e − produced by Dalitz decays of π 0 are shown by black triangles, from η -as red squared, from ∆ -as blue rhombus. The sum of all DM channels (allowed kinematically in a given M = M U bin) is shown by brown stars (indicated as "Sum U" in the legend) and the sum of dileptons from all SM channels and U-decays is displayed as yellow dots ("Sum SM+U" in legend). The calculation of dileptons from U -bosons have been done with the factor C U = 0.1, i.e. by allowing 10% surplus of the theoretical results from the total sum of all SM channels.
In Fig. 2 we show the PHSD result for the mass differential dilepton spectra dN/dM -normalized to the π 0 multiplicity -for the SM and DM channels for minimum bias Au+Au collisions at 4.0 A GeV which is a prediction for future FAIR and NICA experiments. The description of the individual lines are the same as in Fig. 1. Here we do not apply any acceptance cuts. Also we consider for the dark photons C U = 0.1, i.e. 10% allowed surplus of the total SM yield. Additionally, we show the sum of all U -boson contributions with 50% allowed surplus of the SM yield by the grey open dots.
As follows from Figs. 1 and 2 the total contributions from DM sources follow the shape of the corresponding SM contributions. This follows from the constraints on 2 (M U ) given by Eq. (25) which "scale" the DM yield by the same factor C U in each M U bin. One can see that with increasing beam energies the contribution of vector mesons becomes more and more important for M > 0.6 GeV. As advocated in Refs. [10,11] the U -bosons can be also produced by vector meson conversion due to the U (1) − U (1) mixing; these channels are not accounted for in the present study, however, are a subject of future investigation. Finally, in Fig. 3 we show the results for the kinetic mixing parameter 2 versus M U extracted from the PHSD dilepton spectra for p + p at 3.5 GeV (red line), p + N b at 3.5 GeV (magenta line), Ar + KCl at 1.76 A GeV (green line), Au + Au at 1.23 A GeV (blue line), Au + Au at 4.0 A GeV in comparison with the combined HADES results (orange line) from Ref. [15]. The PHSD results are shown with 10% allowed surplus of the Uboson contributions over the total SM yield. The grey dashed vertical line limits the region of applicability of our estimates of 2 , i.e. M U ≤ 0.6 GeV, since at larger dilepton masses the contributions from vector mesons become important. However, we still show our calculations for larger M U to demonstrate that our theoretical method can provide useful constraints on 2 for any M U . Note also that the shape of our theoretically obtained 2 (M U ) curve is not affected by the experimental detector acceptance since the latter acts in the same manner on the SM and DM contributions at a given M = M U . One can see that the extracted 2 for small (M U ) < m π 0 depends only slightly on the size of the collision system and energy since, the π 0 Dalitz decay is the dominant channel. With increasing M U more channels are open and the result depends on the fraction of the U -boson production channels relative to all other dilepton channels.
As follows from Fig. 3, our 2 (M U ) results calculated with C U = 0.1 are close to the HADES experimental result [15] for 0.15 < M U < 0.4 GeV. We point out that this HADES finding is above the upper limit provided in the compilation of world wide experimental results [4]: the measurements by the NA48/2 [17], BaBar [18,19] and A1 [16] experiments set the present world limit at about 10 −6 in this mass range.
In order to explore the theoretical uncertainties in setting the upper limit we present in Fig. 4 (similar to Fig. 1) -the PHSD results for the differential cross section dσ/dM for e + e − production in p + p reactions at 3.5 GeV beam energy calculated for different 2 scenarios in comparison to the experimental measurements by the HADES Collaboration [62] : i) C U corresponding to  1: The PHSD results for the differential cross section dσ/dM for e + e − production in p + p (upper, left) and p + N b reactions (upper, right) at 3.5 GeV beam energy and for the mass differential dilepton spectra dN/dM -normalized to the π 0 multiplicity -for Ar + KCl collisions at 1.76 A GeV (lower, left) and for Au + Au collisions at 1.23 A GeV (lower, right) in comparison to the experimental measurements by the HADES Collaboration. The solid dots present the HADES data for p + p [62], for p + N b [63,64], for Ar + KCl [65] and for Au + Au [66], respectively. The individual colored lines display the contributions from the various SM channels of dilepton production in the PHSD calculations (cf. color coding in the legend). The contributions from U → e + e − (with 10% allowed surplus of the total SM yield) produced by Dalitz decays of π 0 are shown as black triangles, of η as red squared, of ∆-resonance as blue rhombus, their sum -as brown stars and the sum of dileptons from all SM channels and U-decays -as yellow dots. The theoretical calculations are passed through the corresponding HADES acceptance filter and mass/momentum resolution. 0.1% and 10% (10% as in Fig. 1) allowed surplus of total SM yield; ii) constant 2 = 10 −5 and 2 = 10 −6 (cf. color coding in the legend). We note that the shape of the dσ/dM dileptons from U -bosons decays shown in Fig. 4 for 2 = 10 −6 is different from the shape of the same contributions presented in Fig. 1 for C U = 0.1 (i.e. 10% surplus over the total SM channels). One can clearly see the 'step' behaviour at the mass ∼ 0.21 GeV which reflects the mass dependence of the branching ratio for the decay of U -bosons to e + e − (Eq. (22)) and is related to the opening of the µ + µ − channel (cf. Fig. 1 of Ref. [15]). We find that the dN/dM spectra calculated with constant 2 = 10 −6 are coincident with those calculated with C U = 0.001 which is 0.1% surplus in the 0.15 < M < 0.6 GeV mass range. On the other hand, the 10% surplus -shown in Figs. 1 and 4 and discussed above -corresponds to a constant 2 = 10 −5 in the 0.15 < M < 0.6 GeV mass range. We mention that both scenarios for the extraction of the upper limit of the mixing parameters are based on the comparison of theoretical resultsassuming that the theoretical model reproduces the SM contributions "ideally" which is not exactly the case. The underlying uncertainties come from the calculation of the Dalitz decays of mesons and ∆ resonances due to the lack of knowledge of the corresponding formfactors as well as a relative contribution of individual channels (cf. discussions in Section II).
For the illustration of the discussion above we show in Fig. 5 the kinetic mixing parameter 2 extracted from the PHSD dilepton spectra for p + p at 3.5 GeV calculated for different 2 scenarios in comparison with the combined HADES results (orange line) [15] (as in approximately corresponds to the present knowledge on the upper limit based on the compilation of the world wide experiments [4]. We have to stress that the leverages applied to the extraction of an upper limit for 2 (M U ) in our calculations or in the HADES experiment are rather different: while in our theoretical approach, the excess caused by a hypothetical U -boson is shown as a function of the (yet unkown) U -boson mass (for a given fixed probability taken here as 10%), the experimental approach is to search for an unexpected peak structure of given width (defined by The PHSD results for the differential cross section dσ/dM for e + e − production in p + p reactions at 3.5 GeV beam energy calculated for different 2 scenarios in comparison to the experimental measurements by the HADES Collaboration [62]: i) CU corresponding to 0.1% and 10% (as in Fig.  1) allowed surplus of total SM yield; ii) constant 2 = 10 −5 and 2 = 10 −6 (cf. color coding in the legend). FIG. 5: The kinetic mixing parameter 2 extracted from the PHSD dilepton spectra for p + p at 3.5 GeV calculated for different 2 scenarios in comparison with the combined HADES results (orange line) [15] (as in Fig. 3). The PHSD results are shown with 0.1, 5, 10, 15% allowed surplus of the Uboson contributions over the total SM yield (cf. color coding in the legend). The dotted magenta line shows the constant 2 = 10 −6 .
the detector resolution) on top of a smooth continuum. This also means that in the experiment the sensitivity to find such a peak is basically limited by the statistical fluctuations of the measured total dilepton yield. Except for the π 0 Dalitz region the measured yield is dominated by the irreducible combinatorial background, i.e. the backgorund of uncorrelated lepton pairs. Ultimately, the pair statistics determines the uncertainties at higher invariant masses and the only way to reduce the impact of these fluctuations is by either increasing the detection efficiency, extending the experimental run time and/or increase the event rates, i.e. run with larger beam inten-sities, this is planned by the CBM experiment at FAIR [67].

V. SUMMARY
In this study we presented the first microscopic transport calculations, based on the PHSD approach, for the dilepton yield from the decay of hypothetical dark photons (or U -bosons), U → e + e − , from p + p, p + A and heavy-ion collisions at SIS18 energies. For that we incorporated in the PHSD the production of U -bosons by the Dalitz decay of pions π 0 → γU , η-mesons η → γU and Delta resonances ∆ → N U decays based on the theoretical model by Batell, Pospelov and Ritz from in Refs. [10,11] which describes the interaction of DM and SM particles by the U (1) − U (1) mixing. The strength of these interactions is defined by the kinetic mixing parameter 2 which is a parameter in the model. Moreover, the mass of the U -boson M U is also unknown, i.e. 2 = 2 (M U ). We introduce a procedure to define theoretical constraints on the upper limit of 2 (M U ): i) Since the dark photons are not observed in dilepton experiments so far (e.g. they were not found in the HADES dark photon search [15]), we can require that their contribution can not exceed some limit which would make them visible in experimental data. ii) While the dilepton spectra from the SM channels are in a good agreement with the HADES experimental data for p + p, p + N b at 3.5 GeV, Ar + KCl at 1.76 A GeV and Au + Au at 1.23 A GeV, we can impose a constraint that the extra contribution from the U -boson decay in each bin of invariant mass M = M U can not exceed the SM yield by more than some small factor C U which we can vary in our calculations. iii) By setting C U to e.g. 0.1 (which implies that we "allow" 10% of extra dark photon yield at each M U ), we obtained constraints on 2 (M U ) which can be calculated according to Eq. (25) as a ratio of the total SM yield over the total dark photon yield when setting formally = 1. This procedure allowed to estimate the 2 (M U ) mixing parameter in the mass interval M U < 0.6 GeV based on pure theoretical results for dilepton spectra from SM sources and dark photon decay with any 'requested' accuracy. We have demonstrated how the result will change by reducing the surplus of DM contribution from 0.1% to 50% over the SM spectra.
We found that the upper limit of 2 (M U ) extracted for 10% surplus of DM contributions over SM sources is consistent with the experimental results of the HADES experiment [15], for 0.15 < M U < 0.4 GeV, i.e. the calculations are on a level of the experimental upper errorbars.
However, the HADES limit on 2 (M U ) has been superseded by the NA48/2, BaBar and A1 experiments: the modern world wide experimental compilation [4] shows that 2 ∼ 10 −6 in this mass range. In fact, the meaning of the quoted upper limit is that 2 ≤ 10 −6 at a confidence level of 90%. This would lead to only 0.1% allowed surplus of the DM over SM contributions. Indeed, such an accuracy is very difficult to achieve in theoretical calculations due to the model uncertainties as well as statistical fluctuations during the simulations. Still, we believe that our model allows to run very realistic simulations, in particular of dark photon production in the environment of a heavy-ion collision, useful for optimizing the analysis of existing data and for planning future, more sensitive measurements.
We note that in the present study we focussed only on light dark photons with M U < 0.6 GeV since the considered Dalitz decay channels for U -boson production are dominant in this mass range, especially for the low energy reactions as measured by the HADES experiment. However, in future our study can be extended to larger M U masses which would require to account for the production of U -bosons by other channels such as direct decays of vector mesons and Drell-Yan type of processes which are dominant for dilepton production at intermediate and high dilepton mass regions.
Furthermore, we predicted limits for the dark photon contribution to the dilepton yield for Au+Au collisions at 4.0 A GeV, which is of interest for the upcoming FAIR [67] and NICA [68] experiments. Since the experimental procedure of searching for U -boson dilepton decays is related to the observation of a sharp peak on top of a smooth continuum, it is also constrained by statistical fluctuations of the combinatorial background. In this respect, it may be more promising to search for peak structures in dilepton experiments optimized for a maximal signal/background ratio, i.e. by using p+A reactions or collisions of light nuclei to reduce the total π 0 multiplicity per event.