Search for Majoron-emitting modes of $^{136}$Xe double beta decay with the complete EXO-200 dataset

A search for Majoron-emitting modes of the neutrinoless double-beta decay of $^{136}$Xe is performed with the full EXO-200 dataset. This dataset consists of a total $^{136}$Xe exposure of 234.1 kg$\cdot$yr, and includes data with detector upgrades that have improved the energy threshold relative to previous searches. A lower limit of T$_{1/2}^{\rm{^{136}Xe}}>$4.3$\cdot$10$^{24}$ yr at 90\% C.L. on the half-life of the spectral index $n=1$ Majoron decay was obtained, a factor of 3.6 more stringent than the previous limit from EXO-200, corresponding to a constraint on the Majoron-neutrino coupling constant of $|\langle g_{ee}^{M}\rangle|$$<(0.4$-$0.9)\cdot10^{-5}$. The lower threshold and the additional data taken resulted in a factor 8.4 improvement for the $n=7$ mode compared to the previous EXO search. This search provides the most stringent limits to-date on the Majoron-emitting decays of $^{136}$Xe with spectral indices $n=1,2,3,$ and 7.

A search for Majoron-emitting modes of the neutrinoless double-beta decay of 136 Xe is performed with the full EXO-200 dataset. This dataset consists of a total 136 Xe exposure of 234.1 kg·yr, and includes data with detector upgrades that have improved the energy threshold relative to previous searches. A lower limit of T 136 Xe 1/2 >4.3·10 24 yr at 90% C.L. on the half-life of the spectral index n = 1 Majoron decay was obtained, a factor of 3.6 more stringent than the previous limit from EXO-200 and a factor of 1.6 more stringent than the previous best limit from KamLAND-Zen. This limit corresponds to a constraint on the Majoron-neutrino coupling constant of | g M ee |< (0.4-0.9) · 10 −5 . The lower threshold and the additional data taken resulted in a factor 8.4 improvement for the n = 7 mode compared to the previous EXO-200 search. This search provides the most stringent limits to-date on the Majoron-emitting decays of 136 Xe with spectral indices n = 1, 2, 3, and 7. Double-beta (ββ) decay is a rare weak transition between two nuclei with the same mass number A and nuclear charges Z that differ by two units. The process is only observable if single-beta (β) decay is highly suppressed or forbidden by energy conservation. Decays in which two neutrinos are emitted (2νββ) are an allowed process in the Standard Model, and have been observed in a number of nuclides [1] including 136 Xe with a half-life of T 2νββ =[2.165±0.016 (stat) ±0.059 (syst)]×10 21 yr [2]. However, if neutrinos are massive Majorana fermions, ββ decays can also proceed without emission of neutrinos, violating lepton number conservation [3]. The simplest of such modes, the 0νββ decay with the emission of two electrons and nothing else, is a subject of an intense experimental search [4]. The most recent measurements have set stringent lower limits on the half-life for 0νββ decay of several isotopes, including 136 Xe (EXO-200 [5] and KamLAND-Zen [6]), 76 Ge (GERDA [7]), and 130 Te (CUORE [8]).
In this paper we present results of a search for neutrinoless ββ decay modes of 136 Xe in which one or two additional bosons, denoted as χ 0 here, are emitted together with the electrons, i.e.: Any bosons emitted in the 0νββχ 0 (0νββχ 0 χ 0 ) modes are usually referred to as "Majorons". Majorons were originally proposed as a Goldstone boson associated with spontaneous lepton number symmetry breaking [9][10][11]. Majorons are possible dark matter candidates [12], and may be involved in cosmological and astrophysical processes [13,14]. Precise measurement of the width of the Z boson decay to invisible channels [15] has disfavored the original Majoron models. However, other analogous models have been proposed, free of this constraint, in which Majorons more generally refer to massless or light bosons that might be neither Goldstone bosons, nor be required to carry a lepton charge [16]. The spectral index n is used to characterize different Majoron-emitting modes that are experimentally recognizable by the shape of the sum electron spectrum [17,18]. A novel model of 0νββ decay with a emission of a light Majoron-like scalar particle φ was also proposed in [19], where the Majoronlike particle couples via an effective seven-dimensional operator with a right-handed lepton current, and with right-handed ( φ RR ) or left-handed ( φ RL ) quark current. The normalized spectra for various Majoron-emitting decay modes of 136 Xe are shown in Fig. 1. The calculation of the spectra uses the Fermi function suggested in [20] that fully includes the nuclear finite size and electron screening [21] and evaluates the value of the Fermi function at the nuclear radius R as recommended in [22]. The single-state dominance model of 2νββ is used, while the higher-state dominance models that yield slightly different spectral shapes [23,24] are not considered in the analysis.
Recent sensitive searches for Majoron-emitting ββ decays have been carried out in 76 Ge (GERDA [25]), 130 Te (CUORE [26]), and 136 Xe (KamLAND-Zen [6] and EXO-200 [21]). EXO-200 reported a lower limit of T 136 Xe 1/2 > 1.2 × 10 24 yr at 90% C.L. on the half-life of the spectral index n = 1 Majoron decay mode based on 100 kg·yr exposure of 136 Xe [21], compared with a lower limit of 2.6 × 10 24 yr at 90% C.L. reported by KamLAND-Zen [6]. Following that analysis, several detector upgrades were made to EXO-200 permitting a lower analysis threshold, and additional data were acquired in "Phase-II" of EXO-200 operations from May 2016 to December 2018 utilizing these technical improvements. The total Phase-II exposure collected was similar to that of the first run ("Phase-I", September 2011 to February 2014) from which the previous searches for Majoron-emitting modes were reported. This paper reports a search for Majoron-emitting modes of ββ decay using the full EXO-200 dataset that totals 1181.3 days of live time after data quality cuts described in [27], corresponding to a 134% increase in 136 Xe exposure relative to the previous search [21].

II. DETECTOR DESCRIPTION, DATA AND MONTE CARLO SIMULATION
The EXO-200 detector consisted of a cylindrical liquid xenon (LXe) time projection chamber (TPC) filled with LXe enriched to 80.6% in 136 Xe. A cathode split the TPC into two drift regions, each with a radius of ∼18 cm and a drift length of ∼20 cm. The TPC was enclosed by a radiopure thin-walled copper vessel. The electric field in the drift region was raised from 380 V/cm in Phase-I to 567 V/cm in Phase-II to improve the energy resolution. The ionization produced from interactions in each drift region was read out after being drifted to crossed-wire planes at each anode, inducing signals on the front-most wire plane (V-wires), after which it was collected by the second wire plane (U-wires). The scintillation light was collected by arrays of large area avalanche photo-diodes (LAAPDs) [28] located behind the wire planes. A more detailed description of the detector can be found in [29,30].
The detector was located inside a clean room at the Waste Isolation Pilot Plant (WIPP) in Carlsbad, New Mexico, under an overburden of 1623 +22 −21 meters water equivalent [31]. An active muon veto system consisting of plastic scintillator panels surrounding the clean room on four sides allowed prompt identification of >94% of the cosmic ray muons passing through the setup, and allowed rejection of cosmogenic backgrounds [31].
Each TPC event is reconstructed by grouping charge and light signals into individual energy deposits. Ionization signals measured by the two wire planes provide information about coordinates x and y perpendicular to the drift field. The time difference between the light signal and the associated charge signal, together with the measured electron drift velocity [32], provides the z position. Events with a single reconstructed charge deposit are referred to as "single site" (SS), and include most β or ββ decays with characteristic spatial extent of 2-3 mm. Events with multiple reconstructed deposits are referred to as "multisite" (MS), and arise mostly from multiple interactions of MeV-energy γ-rays. Additionally, internally generated ββ-like events in the fiducial volume (FV) are uniformly distributed in the LXe, compared to the spatial distribution of background events arising from γ-rays entering the TPC, which tend to be concentrated nearer to the vessel walls. This difference is captured in the analysis by the standoff-distance (SD) variable, defined as the shortest distance between any reconstructed charge deposit and the closest material surface excluding the cathode. The total energy of an event is determined by combining the charge and scintillation signals. This combination achieves better energy resolution than possible from each individual channel alone due to the anticorrelation between them [33,34].
The detector response to ββ decays and background interactions is modeled by a detailed Monte Carlo (MC) simulation based on GEANT4 [35]. Radioactive γ sources are deployed at several positions near the TPC to characterize the detector response and validate the MC simulation. The scintillation and ionization yields were determined with γ interactions from calibration sources over a range of electric fields [34]. The energy scale and resolution are simultaneously determined by fitting the expected energy spectra generated by MC to the corresponding calibration γ sources. The absolute ββ energy scale is found to be consistent with the calibration γ sources at the sub-percent level [5,34].
The dataset and event selection criteria used in this work is the same as in the search for 0νββ decay [5], except that a reduced energy threshold is used here. The 136 Xe exposure of the entire dataset after data quality cuts and accounting for live time loss due to vetoing events coincident with the muon detector is 234.1 kg·yr, or 1727.5 mol·yr, with 117.4 (116.7) kg·yr in Phase-I (Phase-II).

III. ANALYSIS PROCEDURE
The energy threshold used in this work is lowered from the 1000 keV threshold used in the most recent 0νββ search by EXO-200 [5] to 600 keV, to improve sensitivity to Majoron-emitting modes at low energy. As shown in Fig. 2, the energy measurement shows good spectral agreement between the data and the simulation for SS and MS events above 600 keV using 228 Th, 226 Ra, and 60 Co calibration sources. The SD is also observed to have good agreement between the data and the simulation at energies above 600 keV, as shown in Fig. 2.
Probability density functions (PDFs) for signal and background components are generated using the Monte Carlo (MC) simulation. The PDFs are functions of two observables: energy and SD. Studies were performed to estimate the sensitivity improvement possible with additional multivariate discriminators similar to those used in previous 0νββ searches [5,36], but the minimal set of variables (energy and SD) was chosen as it had com-  parable sensitivity while minimizing systematic errors at low energy. The systematic errors in the analysis are described in detail in Sec. IV. To avoid any possible bias in analysis criteria, all selection cuts and the choice of fitting variables were determined from MC-based sensitivity studies alone prior to performing any fits to the data itself. The components of the overall fit model are similar to that in [5] with the 0νββ signal replaced by a Majoron-emitting decay. Since this work used a reduced energy threshold, two components are added to the background model: • 85 Kr dissolved in the LXe that produces β decays with an end point energy of 687.0 keV [37] • 137 Cs in the materials near the LXe, which emits γ-rays with energy of 661.7 keV [37] The simulation of 85 Kr includes the two β decay modes with branching ratios of 99.56% and 0.434% to the ground and excited states of 85 Rb followed by the release of a 514.0 keV γ-ray, respectively [37]. A shape correction accounting for the forbidden nature of the first unique β decay was calculated using the method described in [38] to be between -15% and 80% depending on its energy. This correction was applied as an event weighting in the MC simulation. The PDF model is parametrized by the event counts and SS fractions (f SS = SS/(SS+MS)) of the individual components, as well as two normalization parameters to account for the effects of systematic errors [39]. The search was performed using a maximum-likelihood (ML) function to fit simultaneously both SS and MS events with their corresponding PDFs generated by MC, in a similar approach as [5]. Systematic errors, described in Sec. IV are included in the ML fit as nuisance parameters constrained by normal distributions. The median 90% CL sensitivity is estimated using toy datasets generated from the PDFs of the background model. An energy threshold of 600 keV is used in the fit, which provides near optimal sensitivity for all Majoron-emitting modes considered here. The reduced energy threshold relative to previous analyses [5,21] results in higher signal detection efficiencies, especially for the n = 7 mode, which has a peak around 628 keV in the energy spectrum. Lower energy thresholds do not further improve sensitivity because the increase in signal efficiency is outweighed by the presence of backgrounds, including 85 Kr.

IV. SYSTEMATIC ERRORS
Systematic errors were accounted for by the same technique as described in [5,39]. The five Gaussian constraints added to the ML fit, which are used to propagate systematic errors into the results, correspond to: • uncertainty in the activity of radon in the LXe as determined in stand-alone studies via measurement of time correlated 214 Bi-214 Po decays [40]; • uncertainty in the relative fractions of neutron capture-related PDF components generated by dedicated simulations [31]; • uncertainty in SS fractions obtained by comparisons between calibration data and MC; • uncertainty in the overall event detection efficiency, referred to as normalization, caused by event reconstruction and selection efficiencies; • uncertainty in the signal detection efficiency, referred to as signal-specific normalization, caused by discrepancies in the shape distributions between data and MC and background model uncertainties.
The first two uncertainties were constrained by relative errors of 10% and 20%, respectively, as evaluated in [41]. The uncertainty in the SS fractions is determined by comparisons between the data and MC for calibration sources, as shown in Fig. 3. Taking into account different calibration sources at various positions, these systematics are evaluated to be 3.7% (3.6%) for Phase-I (Phase-II), averaged over the energy range considered here. The uncertainty on the overall efficiency was evaluated to be 3.1% (3.0%) for Phase-I (Phase-II), and differs slightly from that estimated in [5] after accounting for the larger energy range considered here. The top four uncertainties are presented in Table I. Discrepancies in the shape distributions between data and MC are propagated into the signal rate through a normalization parameter that only scales the coefficient of the signal PDFs. This signal-specific normalization parameter is constrained to unity within the errors arising from spectral shape agreement and background model uncertainties as described below. To estimate the effect of spectral shape errors, the ratio between data and MC of the projections onto energy and SD (Fig. 2) were used to re-weight all PDF components by the observed differences (referred to below as un-skewing). 60 Co and 238 Urelated PDFs were weighted by ratios from 60 Co and 226 Ra calibration sources respectively, while the other γlike PDFs were weighted by ratios from the 228 Th source that has the most data. Toy datasets were generated from these un-skewed background PDFs, along with a given number of signal events. These toy datasets were then fit with the standard background and signal PDFs used in the primary analysis. The average difference between the true number of signal events added to the toy datasets and the best-fit signal counts is used to quantify the impact of the spectral discrepancy. To evaluate   the uncertainty associated with the background model, decays of Th, U, and Co were simulated in different locations than that in the default model. For example, all far 238 U are represented by the decays in the air gap between the cryostat and the lead shielding in the background model. To evaluate the errors introduced by this approximation, 238 U simulated in the cryostats is used to represent all 238 U from remote locations. This is taken to represent an extreme deviation from the more realistic case used in this analysis. Toy datasets generated with the default model along with a given number of signal events were fitted with this alternative model, and the resulting difference between the true number of signal events added to the toy datasets and the best-fit signal counts is taken as the systematic error of the background model. The background model error also includes the effects of perturbations to the 2νββ spectrum due to corrections to the Fermi function arising from the finite nuclear size and electron screening effects [20,22]. The 2νββ PDF integrals differed by 1.5% in the case of a differing Fermi function. The signal-specific normalization error is a function of the number of signal counts to account for possible existence of signal events, that is presented in Table II.

V. RESULTS AND CONCLUSION
To determine confidence intervals on each of the Majoron decay modes considered here, the datasets from Phase-I and Phase-II are fit independently with efficiency and livetime from each phase taken into account. The observed data, and an example of the best fit spectrum for the n = 1 Majoron mode, are shown in Fig. 4. No sta- tistically significant evidence for Majoron-emitting 0νββ decays is observed for any mode considered. The fits are initially performed individually for each Majoronemitting decay mode separately for Phase-I and Phase-II, and the combined limits are determined after summing the profile likelihood obtained from each dataset. While the fits are performed for each mode independently, Fig. 4 also overlays the corresponding spectra from all fits at the 90% C.L. upper limits on the number of decays for each mode. The lower limits on the Majoron-emitting 0νββ half-lives are summarized in Table III. EXO-200 has better detection efficiency at low energies and improved spectral agreement between data and MC simulation in Phase-II with upgraded electronics, which result in more stringent lower limits on most Majoron-emitting 0νββ half-lives than in Phase-I. The improvement is particularly significant for n = 7 mode that has a spectrum peaking at lower energy than the other modes.
The lower limit on the Majoron-emitting 0νββ halflives for the models with n =1, 3, and 7 can be translated into limits on the effective neutrino-Majoron coupling constants g M ee using: 1 where M =M ( g A 1.25 ) 2 , M is the nuclear matrix element, g A is the axial coupling constant, m = 2 (4) for the emission of one (two) Majorons, and G 0νM (Z, E 0 ) is the unnormalized phase space integral that depends on the model type and fundamental constants [21]. Table III summarizes the 90% C.L. upper limits on the effective Majoron-neutrino coupling constants. The phase space factors are taken from [21], while the matrix elements are taken from [42,43] for the Majoron decay mode with n = 1, and from [44] for other modes. The spread in the limits on the coupling constants is due to ambiguity in the matrix elements. This is the most stringent limit on g M ee for the n = 1 Majoron among all ββ decay nuclei [25,26,45,46]. The previous best limit from a laboratory experiment comes from KamLAND-Zen, which reported a limit of g M ee <(0.8-1.6)·10 −5 [45]. The phase-space integral for the n = 1 Majoron used by KamLAND-Zen is about a factor two smaller than the most up to date value used here [21]. KamLAND-Zen's half-life limit would translate to a limit on the coupling constant of g M ee <(0.6-1.2)·10 −5 with the same phase space factor, and our new limit corresponds to a factor of 1.3 improvement over this previous result. Our limit on the coupling constant of g M ee is two orders of magnitude more stringent than the recent result obtained in the measurement of pion decays [47].
In conclusion, we report results from a search for Majoron-emitting double-beta decay modes of 136 Xe with the complete EXO-200 dataset. No statistically significant evidence for this process is found, and we obtain limits on half-lives and effective coupling constants of Majoron-emitting 0νββ decay modes that are more stringent than results of previous EXO-200 [21] and KamLAND-Zen [45] searches.