Constraining Non-minimal Dark Sector Scenarios with the COHERENT Neutrino Scattering Data

Abelian dark sector scenarios embedded into the two-Higgs doublet model are scrutinized within the Coherent Elastic Neutrino-Nucleus Scattering (CE$\nu$NS) experiment, which was first measured by the COHERENT Collaboration in 2017 with an ongoing effort to improve it since then and recently released data for the CsI target in 2022. In the theoretical framework, it is assumed that there is a $U(1)$ gauge group in the dark sector with a non-zero kinetic mixing with the hypercharge field. The COHERENT data for the targets CsI and liquid argon (LAr) are treated in both single and multi bin bases to constrain the multi dimensional parameter space, spanned by the dark gauge coupling, kinetic mixing parameter and the dark photon mass, of totally seven different representative scenarios which are also compared and contrasted among each other to find out about the most sensitive one to the data. The effect of refined quenching factor is also addressed.


I. INTRODUCTION
With the discovery of the Standard Model Higgs boson at Large Hadron Collider [1,2], the last missing piece of the Standard Model (SM) was found. While the search for physics beyond the Standard Model at high energies has been continuing, there are low-energy tools which could be used to test physics beyond the SM. Recently measured coherent elastic neutrino-nucleus scattering (CEνNS ) [3][4][5][6] is one of such powerful tools, which has been considered for new physics scenarios [7][8][9][10], for performing precision tests of the Standard Model [10,11], as well as being instrumental for many other implications [12][13][14][15][16][17].
The CEνNS is an elastic scattering process which takes place between low-energy neutrinos (in the MeV range) and atomic nuclei. It is a weak interaction mediated by Z boson with a tiny momentum transfer (in the keV range), which makes the interaction coherent with the nucleus, whose cross section in the SM is scaled as the square of the number of neutrons in the nucleus (∼ N 2 ). This brings an enhancement to the cross section as compared to the other scatterings in the same energy range.
The idea of measuring CEνNS cross section had been first proposed by Freedman [18] and after more than four decades, the COHERENT Collaboration has managed to measure the CEνNS by using the targets Cesium-Iodide (CsI[Na]) first [3,4] and then argon (Ar) [5,6].
There is also a recent update [19] by the COHERENT Collaboration with more data on CsI together with various improvements in the analysis like using an updated energy-dependent quenching model, including energy-smearing effect and time dependent efficiency. There is a very recent analysis [20] combining full CsI data with the available LAr results which considers various implications like testing the SM, studying the electromagnetic properties of neutrinos (see also [21]) and constraining some new physics scenarios by allowing generalized types of neutrino-neucleus interactions while assuming universal couplings.
The discovery of the SM Higgs as the first elementary spin-0 scalar has boosted the interest in extending the scalar sector of the SM, the simplest of which is known as Two Higgs Doublet Models (2HDM). Within the 2HDM scenario it is possible to address both the flavor changing neutral current problem and mass for neutrinos [39] by extending the gauged sector of the SM with a U (1) D group together with right-handed neutrinos added to the spectrum. Various anomaly free models have been constructed and discussed in highand low-energy scales [39]. The use of CEνNS data as a new physics probe to constrain the 2HDM gauged with U (1) D is missing in the literature and the current study aimed to provide a thorough investigation on this matter to fill the gap.
The anomaly free extension of the SM (plus right-handed neutrinos) by coupling it with a U (1) B−L gauge group from a hidden sector is a well-established scenario to study dark photon effects through its kinetic mixing with the SM hypercharge gauge field as well as its gauge coupling with the SM fermions. Such a scenario would be naturally a limiting case of the 2HDM mentioned above. All of the viable scenarios will be compared and contrasted with each other in the light of COHERENT data.
In Section II, the theoretical framework is summarized and the relevant SM vertices which receive corrections or the new vertices are also listed for completeness. The cross section formulas contributing the CEνNS and some basics about the COHERENT data are given in Section III and Section IV, respectively. After having the statistical method briefly explained, the numerical study is conveyed in Section V and we conclude in Section VI.

II. A NONMINIMAL DARK SECTOR FRAMEWORK
One way of realizing interactions between the hidden (dark) sector and the visible (SM) one is through various portals which are dimension-4 operators and hence are free of suppression. A so-called dark vector boson from the hidden sector coupling with the SM is one popular scenario, known as the vector portal.
In a minimal scenario, such a coupling includes only a kinetic mixing between the dark vector boson and the weak hypercharge field, which allows the dark vector boson indirect access to the SM fermions. If the vector boson of the dark sector is a gauge field of a dark group, say U (1) D (again in the minimalistic approach), having the SM fermions charged under U (1) D would allow them to get a direct coupling with the dark vector boson. This will require checking additional gauge-anomaly conditions which restrict the dark quantum charges of the SM fermions. B − L is one popular choice for quantum charges which makes the overall scenario anomaly free when the SM is extended with the right-handed neutrinos.
In the above framework, both the visible and dark sectors are assumed to be minimal for the sake of taking advantage of mainly the predictive power of the scenario. Either or both sectors could be enlarged, when the SM being an effective theory is especially considered, which leads to nonminimal scenarios.
One popular way of extending the SM while keeping the dark sector minimal (specifically, assuming an Abelian dark-gauge group U (1) D ) is through its scalar part, simply by adding another SU (2) doublet 1 , known as the two-Higgs doublet models (2HDM) [40] which have motivations from supersymmetry [41], axions [42], baryogenesis [43][44][45][46] etc. On the other hand, adding right-handed neutrinos to the particle content of the construct would not only address the issue of neutrino masses but also usually needed for satisfying anomaly equations due to gauging the 2HDM under the additional dark group U (1) D . The details of the scenarios have been worked out in other studies [39,47,48] (see also [49][50][51][52][53] for earlier theoretical studies including different aspects of the mixing among the gauge bosons) and here we only reproduce the relevant part of the model whose Lagrangian terms are given below Here X 0 µ is the gauge field of U (1) D while Y µ and W 3µ are the usual weak hypercharge field of U (1) Y and the third of weak gauge fields, respectively. ϕ 1 and ϕ 2 in L KE Scalar are the usual scalar doublets under SU (2) L . sin ϵ is the strength of the kinetic mixing between U (1) Y and U (1) D . m X represents the Stueckelberg mass parameter for the dark gauge field X 0 µ , which will be explained briefly.
Unlike in the non-Abelian case, the mass generation for the Abelian gauge bosons does not necessarily require the existence of a Higgs mechanism since there is no issue of unitarity or renormalizability for Abelian gauge theories. Hence alternative ways of mass generation other than the Higgs mechanism can be employed for such sectors like U (1) D while in the SM part, masses are acquired through the Higgs mechanism. One of the popular ways is the so-called Stueckelberg mechanism [54] where the U (1) gauge boson couples with the derivative of an axionic (dark) scalar field in a gauge invariant manner. We assume that the scalar field from the dark sector is charged under U (1) D only so that the m X mass term in Eqn. (2) is obtained. The Stueckelberg extension of the SM with [55] and without the kinetic mixing [56] have been discussed in detail. The two-Higgs doublet model extended with an additional U (1) gauge symmetry has also been studied in [57] by considering the Stueckelberg contribution to the mass terms.
The Lagrangian L Fermion contains the kinetic energies and electroweak interactions of the fermions. At this stage, the relevant part of the covariant derivative in the (Y µ , W 3µ , X 0 µ ) basis would look like where t 3 , Y , and Q ′ D are the SU (2) L generator, weak hypercharge and dark charges of the field that the covariant derivative is acting on, respectively. Here g W , g Y and g D are the corresponding gauge coupling constants.
The details of the scalar sector as well as the discussion of anomaly cancellations of the model are given in [39]. Some of the relevant parts are going to be repeated here. For example, the dark charge assignment of the SM fields satisfying the anomaly conditions are listed in Table I. The mass terms for the gauge fields that will originate from Eqn.
Having said that the original gauge basis in the neutral sector, i.e., (Y µ , W 3µ , X 0 µ ), is not diagonal, the physical basis, say (A µ , Z µ , A ′ µ ), can be obtained by making three successive rotations. The first two of these are the ones to eliminate the kinetic mixing term and  the usual Weinberg angle rotations, which results in the rotated gauge basis (A µ , W 3µ , X µ ), given by in terms of (Y µ , W 3µ , X 0 µ ), At this stage the mass term becomes where the parameters a, b, and c are given Here the functions ∆ i are defined as follows and is the ratio of the VEVs, usually defined in the 2HDMs. Q ϕ i D is the dark charge of the scalar doublet ϕ i , i = 1, 2.
The third and last rotation is needed due to b ̸ = 0 and the rotation angle ξ should satisfy tan 2ξ = 2b/(a − c) and the final form of the overall rotation matrix from the ( The corresponding eigenvalues are These are the mass squares of the photon (A µ ), dark photon (A ′ µ ) and electroweak neutral boson (Z µ ), respectively. The mass-squared expressions for the gauge bosons in the minimal B − L model can be read off from Eqns. (11) and (12) by taking the limit g D → 0.
The SM Z boson mass, M Z , receives corrections due to being mixed with the initial U (1) D gauge field, X 0 µ , and coupling with the Higgs doublets having nonzero dark charges. Even though the analytical expression of M Z in Eqn. (12) differs from the SM mass, , it is observed that the numerical value of M Z lies on the SM prediction for almost all the parameter space as long as m X is not greater than ∼ O(GeV) together with small enough kinetic mixing parameter sin ϵ. the dark photon mass is only determined by the terms proportional to g D . Therefore, for a fixed sin ϵ or g D in the region m X ≤ m critical X , there will be a nonzero minimum value for the dark photon mass. This feature will be useful in the numerical discussion section.
At the end of the section, we prefer to give the list of vertex factors relevant to the CEνNS process. They are all listed in Table II, based on the parametrization given in the following Lagrangian terms, The vertex factors of the models, listed in Table I, namely Model C, D, E, F, G and B − L (including the minimal B − L as well), can be read off from the entries in Table II by plugging the values of the corresponding charges Q ′ u and Q ′ d , given in Table I. The vertex factors listed in Table II can also be reduced to the SM ones by taking the limiting values; Q ′ u,d → 0 and ϵ → 0 (ξ → 0). In this limit, one can see that any dark photon vertex vanishes Table II. The relevant vertex factors contributing to the CEνNS in the two-Higgs doublet models extended with a dark U (1) D group. A shorthand notation is used for the trigonometric expressions.
For example, (s ξ , t ϵ ) stand for (sin ξ , tanϵ) and similar for the others.
as expected.

III. COHERENT ELASTIC NEUTRINO NUCLEUS SCATTERING (CEνNS)
Coherent elastic neutrino-nucleus scattering (CEνNS) includes low-energy neutrino coupling and provides an accessible window to many areas of physics research especially for physics beyond the Standard Model. CEνNS is suitable for analyzing the nonminimal dark sector framework. When scattering of neutrinos from nuclei at rest is considered, if incident neutrino energy is below 50 MeV, the interaction occurs coherently, meaning that the neutrinos interact with the nuclei as a whole rather than with its constituents individually [18].
In the SM, the cross section for elastic scattering is two orders of magnitude larger than inelastic scattering, which makes CEνNS viable for observation. If a new neutral current interaction mediated by a light vector boson exists, it would not be suppressed by SM interactions in this region. In spite of the fact that the earliest experimental proposal to measure CEνNS was rather old [58]; it took lmost four decades to be able to make significant progress on the way of measuring CEνNS cross section, eventually succeeded by the COHERENT Collaboration in 2017 [3]. Hence CEνNS has become one of the important probes for physics beyond the SM since then.
The differential cross section for CEνNS is well established in the literature. The SM predicts a coherent elastic scattering cross section proportional to the weak nuclear charge, For spin-0 and spin-1/2 targets, the differential CEνNS cross section with respect to the nuclear recoil energy T , in the coherent limit where the form factor approaches unity, is given by [7] where G F is the Fermi coupling constant, M is the nucleus mass, J N = 0, 1/2 is the spin of the nucleus, E ν is the incident neutrino energy and Q W is the weak nuclear charge [59,60] given by where g u V and g d V are vector couplings of u and d quarks, Z and N are the atomic and neutron number respectively, and s W is the sine of the weak mixing angle. The differential cross section is expressed in nuclear recoil energy because detectors measure this quantity.
The effect of the extra T 2 /E 2 ν term appearing in the spin-1 2 case is negligible. Therefore it is convenient to work on the spin-1 2 case since more developed computational tools are available for fermionic particles. The same applies to the dark sector extended models.
The differential cross section with respect to the nuclear recoil energy T in the minimal B − L model is given by where g D = g BL and A = N + Z is the so-called atomic mass number of the nucleus. F (q 2 ) is the Helm-type nuclear form factor [61] and q 2 denotes the squared momentum transfer given by q 2 = 2M T . The Helm form factor is where j 1 (x) is the spherical Bessel function and q = √ q µ q µ with q 2 defined above, s ≈ 0.9 fm is the nuclear skin thickness and R 1 is the effective nuclear radius where The differential cross section expression for the 2HDM case is rather long and is not illuminating to present here.

IV. COHERENT DATA FOR CEνNS
COHERENT is a neutrino-based fixed target experiment. Neutrino beams are produced by striking proton beams of pulse ∼ 1 µs to a mercury target at 60 Hz, creating ≈ 5 × 10 20 collisions per day. These collisions produce π − and π + as byproducts which come to rest inside the target. Then, while negatively-charged pions are mostly absorbed by mercury, positively-charged ones decaying through the channel π + → µ + +ν µ , which happens when the pion at rest, results in monochromatic energy of ≈ 29.7 MeV for ν µ 's; they are called prompt neutrinos. The pion decay is then followed by the decay of the antimuon, µ + → e + +ν µ + ν e , which occurs about ∼ 2.2 µs after the pion decay. Hence,ν µ and ν e are called delayed neutrinos. Energy spectra for neutrinos stemming from this decay are continuous up to 52.8 MeV. Each proton collision on target (POT) produces on average 0.08 neutrinos per flavor.
Before proceeding further, we would like to comment on the following possibility 2 . So far in this study, we have assumed that the impact of new physics has been considered only on the scattering of neutrinos from the nucleus while the production channels of neutrinos 2 We thank the anonymous referee for mentioning this possibility. for incoming neutrino flux f να for prompt ν µ and delayed ν e andν µ beams are given by [62,63]: Time-dependent flux density f να t (t) is presented in [64] and after normalizing f να E and f να t , the total neutrino flux can be expressed as where N = rN P OT /4πL 2 . Here L is the distance between the detector and the neutrino source, r is the number of neutrinos per flavor created in each POT collision, and N P OT is the total number of POT collisions throughout the entire experiment, which is 1.76 × 10 23 for the CsI 2017 data, 3.198 × 10 23 for the CsI 2022 data and 1.37 × 10 23 for the LAr-based experiment.
Nuclear recoil energy is picked up by scintillation detectors and converted into photoelectrons (PE). This conversion is parametrized with light yield L Y , which is the amount of PE produced per unit energy. Values provided with the data releases are L Y = 13.348 PE/keVee for CsI and ∼ 4.5 PE/keVee for LAr-based experiment. The number of produced PE can be expressed as Nuclear recoil energy is primarily dissipated through secondary nuclear recoils, and only a small amount of energy is turned into scintillation (or ionization). The Quenching factor QF is the ratio of energy turned into scintillation where the subscript "ee" stands for electron equivalent and E ee is the equivalent energy of a recoiling electron in which the energy is dissipated through only scintillation. There have been different approaches for quenching factor in each data release by the COHER-ENT Collaboration. In the 2017 release, a constant quenching factor QF CsI = 8.78 ± 1.66% was suggested. Later, an energy-dependent model for the quenching factor was proposed in [65], which increases the accuracy of the SM expectation for the event rate. The energy -dependent model describes the scintillation in CsI by slow ions with low-energy approximation to Birks' scintillation model [66] multiplied by an adiabatic factor to account for the behavior of scintillation production cutoff at low nuclear recoil energy limit. The quenching factor takes the form where kB = 3.311 ± 0.075 × 10 −3 g MeV −1 cm −2 and E 0 = 12.97 ± 0.61 keV. And (dE/dr) i is the total stopping power for ions, extracted from the software SRIM-2013 [67]. Both QF proposals for CsI will be considered in the analysis in the following section.
In the recent study by the COHERENT Collaboration [19], the previously proposed quenching factor for CsI [3] was reassessed by including a new scintillation response to nuclear recoil measurement on CsI[Na] crystal. Quenching in the region of interest is modeled as a fourth-degree polynomial fit to the available data [68] E ee = g(T ) = 0.05546 T + 4.307 T 2 − 111.7 T 3 + 840.4 T 4 (24) where the detector response E ee is in MeVee and the nuclear recoil energy T is in MeVnr.
For the LAr-based experiment, the suggested quenching factor is a linear fit to the world data on argon which is given in Fig. 7 of [5].
Collision events are counted in n PE and time bins. For the CsI 2017 and LAr releases, the expected number of events in the ith PE and jth time bin can be written as where N β targ is the number of nuclei in the target, E min ν and E max ν are energy limits of the incident neutrinos with E max ν ≈ 52.8 MeV, and T i min and T i max give the boundary values in energy for the ith bin. Lastly A is the signal acceptance function [4],  (27) and Θ(n P E ) is a modified Heaviside step function defined as In the CsI 2022 analysis of the COHERENT Collaboration, functions for energy smearing and acceptance for both energy and time spectra are provided [19]. Acceptance in the energy − 0.333322 (29) and the time-dependent part of the acceptance function is given by where a = 0.52µs b = 0.0494/µs, (31) and the energy smearing is parametrized in the following form by using the gamma function where "true" stands for the true spectrum which would be the one obtained without the smearing effect, and "reco" stands for the measured spectrum. Here the parameters a = 0.0749/E true ee and b = 9.56 × E true ee depend on the quenched-energy deposition. The energy smearing is normalized by using the following condition Ω reco dT reco P n reco PE (T reco )|E true ee = 1 .
With these effects taken into account, the number of events in the ith PE and jth time bin is given by At the end of the section, our findings for the total number of events in the SM are listed in Table III and some of the results available in the literature are also added for comparison.

V. STATISTICAL ANALYSIS
In this section, we will adopt a χ 2 -fit to study the sensitivity of the COHERENT data to phenomenological parameters in the framework of new physics interactions. For the current where N i meas and N i exp are measured and expected number of events in the ith bin respectively, B i on is the estimated number of background events when the beam is on, α and β are the systematic parameters for the signal rate and B on , respectively. The χ 2 function is minimized over α and β. Here, the statistical uncertainty of the measurement is given where B i SS denotes the steady-state backgrounds. The so-called pull terms presented by the COHERENT Collaboration [3] have the form where σ α and σ β are the fractional uncertainty of α and β, corresponding to 1-sigma variation. The pull term in the above form is observed to lead to unphysical results around its limiting values. This behavior is also noted in [69] and instead, an asymmetric pull term of the form is suggested to use. For completeness, we use both forms of the pull terms in the χ 2 calculation and comment on their effect.
Even though the scattering data are obtained in a multibinned detector, the earlier analyses provided by the COHERENT Collaboration combined all events in a single bin [3] for CsI 2017 and LAr data. Later it is suggested to adopt rather a multibin analysis and indeed in the latest study by the COHERENT Collaboration more than one multibin options have been performed. The χ 2 function for the single bin analysis can be obtained from Eqn. (35) by simply using the total number of events for the signal and the background.

A. Numerical Discussion
In this study it is promised to analyze various two-Higgs-doublet models extended by a        Table I. The results are depicted in Fig. 6. In the upper left graph, the light blue shaded region is allowed 90% C.L. by the COHERENT data in the (Q ′ u , Q ′ d ) parameter space. The wider shaded strip represents the allowed region by the dark photon-mass constraint. The models are also marked on the graph and, for example, for the chosen inputs which are also indicated on the graph, Models E and F are ruled out by the data. The sensitivity of the allowed region on tan β is shown on the upper right graph for tan β = 2 and tan β = 20 and the allowed region is slightly wider for tan β = 2. In the lower row, the allowed regions are indicated for various g D values on the left panel and for various M A ′ values on the right panel. The size of the allowed region has effected significantly by varying g D , which is somehow less pronounced for the variation of the dark photon mass.
Last but not least, in Fig. 7, the exclusion regions from the neutrino-neucleus scattering data of COHERENT Collaboration for the two most promising scenarios, the Model F and Minimal B −L model, are displayed with other laboratory bounds relevant in the MeV-GeV range. The figure is adapted from Fig. 6 in Ref. [70] where the details of the laboratory experiments are given (see Table III there). As seen from Fig. 7 that there is a new region in 30 MeV-1 GeV range which is now probed and excluded by the COHERENT data depending on the value of the kinetic mixing parameter. most part of the tan β interval and indeed the sensitivity entirely disappears for tan β ≳ 5.

VI. CONCLUSION
More than four decades after its first proposal, coherent elastic neutrino nucleus scattering was successfully measured by the COHERENT Collaboration in 2017, which has initiated After briefly explaining the theoretical framework and listing seven representative models in Table I where the scalar doublets are allowed to have nonzero dark charges under the gauge group U (1) D , we go on to provide analytical expressions for the differential cross sections for CEνNS both in the SM and in the 2HDM extended with U (1) D . The analyses of the COHERENT data for CEνNS from the year 2017 to 2022 have been modified in various ways and some of these details including statistical analysis have been explained together with our basics to carry out the numerical study.
In this study we aim to find out the constraints on the parameter space of the 2HDMs extended with U (1) D by using the COHERENT neutrino scattering data. This has been achieved by looking at the effects of different factors and parameters which might play some role in the analysis. These are the dark photon mass M A ′ , kinetic mixing parameter sin ϵ, the dark gauge coupling g D , the ratio of the VEVs tan β and the free dark charges Q ′ u , Q ′ Indeed, in the 2HDMs extended with U (1) D , the dark photon mass receives contributions from the scalar sector being proportional to the parameter g D . No matter how small the kinetic mixing is, there exists g D proportional contribution to the mass in addition to the m X term which brings the behavior at low dark photon mass tail (see for example the right panel of Fig. 5). Therefore the minimal B − L model gives the best bound for M A ′ ≥ 1MeV while the other models do much better for the lighter dark photon region. This might be taken as a way to distinguish them. It is also observed that the best bounds are obtained for the single bin case (1PE-1t), for the CsI 2022 data and for the QF taken to be constant (even though energy-dependent QF has been used throughout the numerical analysis as proposed by the COHERENT Collaboration in their latest analysis [19]). In a set of plots, Fig. 6, the allowed bands on the (Q ′ u , Q ′ d ) plane have been shown and this could be used as a reference for better assessment of the scenarios beyond the chosen representative ones listed in Table I. In a final plot, our results are overlaid on a global picture and a new region could further be excluded in the 30 MeV-1GeV dark-photon mass range by the COHERENT data depending on the value of the kinetic mixing parameter. Further data available at low energies may either probe new physics scenarios like 2HDM better or even point out a smoking gun signal