The first dual-phase xenon TPC equipped with silicon photomultipliers and characterisation with $^{37}$Ar

For the first time, a small dual-phase (liquid/gas) xenon time projection chamber was equipped with a top array of silicon photomultipliers for light and charge readout. Here we describe the instrument in detail, as well as the data processing and the event position reconstruction algorithms. We obtain a spatial resolution of ~1.5 mm in the horizontal plane. To characterise the detector performance, we show calibration data with internal $^{83\text{m}}$Kr and $^{37}$Ar sources, and we detail the production of the latter as well as its introduction into the system. We finally compare the observed light and charge yields down to electronic recoil energies of 2.82 keV to predictions based on NEST v2.0.


Introduction
Detectors that use liquid xenon are widely employed in the field of astroparticle physics. In particular, dualphase (liquid/gas) time projection chambers (TPCs) are used for direct dark matter detection, searches for the neutrinoless double beta decay of 136 Xe, as well as other rare event searches [1,2].
A defining requirement to observe a rare interaction is an ultra-low background rate, to maximise the signal-to-background ratio and thus the sensitivity to a given search channel. This is achieved via the reduction of the radioactivity levels of detector materials, the definition of central detector regions via fiducialisation to benefit from the self-shielding power of liquid xenon, and the distinction between single and multiple interactions. The latter rely on precise position reconstruction capabilities of the interaction sites.
In a dual-phase TPC, particle interactions are detected by observing prompt scintillation light (S1) and a e-mail: patricia.sanchez@physik.uzh.ch b e-mail: kevin.thieme@physik.uzh.ch a delayed (S2) signal caused by ionisation electrons, drifted and extracted into the gaseous phase where they produce electroluminescence. Both S1 and S2 signals are observed by photosensors. The depth of an interaction is given by the drift time between S1 and S2, while the (x, y)-position is derived from the light distribution in the photosensor plane placed in the gas phase above the liquid.
Many past dark matter searches based on liquid xenon TPCs [3][4][5] and experiments under construction [6,7] employ VUV-sensitive photomultiplier tubes (PMTs) which have been optimised for low levels of radioactivity and are up to 3 inches in diameter [3,[8][9][10]. However, they are among the dominant sources of electronic (ER) and nuclear recoil (NR) backgrounds, and they limit the (x, y)-position resolution, which is correlated to their size. Next-generation detectors such as DARWIN [11] and nEXO [12,13] aim to reduce the backgrounds further, and consider the use of silicon photomultipliers (SiPMs) for one or both photosensor planes.
Motivated by these future experiments, and by our previous studies of VUV-sensitive SiPM arrays in single-phase detectors [14], we have built and characterised the first dual-phase xenon TPC equipped with a top photosensor array containing 16 SiPM channels read out individually. The new detector is based on our previously described Xurich II TPC [15], with several modifications as detailed below. The calibration of the detector, i.e. the determination of its charge and light response parameters was performed with two sources, 83m Kr and 37 Ar, which allow us to investigate ER interactions at low energies.
This article is structured as follows: We describe the detector in Section 2 and the 37 Ar calibration source in Section 3. We detail the data taking and process-ing tools in Section 4, and show the analysis of data acquired with both internal calibration sources in Section 5. We present and discuss the results in Section 6, and summarise the main findings, addressing some challenges facing future detectors in Section 7.

The Xurich II detector with SiPM readout
The current TPC is an upgrade of the Xurich II detector operated at the University of Zurich (UZH) and described in [15]. The upgrade focuses on two aspects: the replacement of the top PMT with an array of SiPMs and their read-out electronics, and the design of a stand-alone setup to introduce the gaseous 37 Ar calibration source into the xenon recirculation loop.
The detector is placed inside a vacuum-insulated stainless steel vessel which is coupled to a liquid nitrogen bath via a copper cold finger. The TPC drift region is defined by a 31 mm×31 mm (diameter × height) polytetrafluoroethylene (PTFE) cylindrical shell with a cathode mesh at the bottom and a gate mesh at the top, shown in Figure 1. The anode is placed 4 mm above the gate. The liquid surface was kept in the middle of the two top electrodes, at ∼ 2 mm above the gate for the presented data, optimised for S2 amplification. Seven copper field-shaping rings, separated by PTFE spacers and connected via a resistor chain, ensure a uniform drift field, which is maintained between the negatively biased cathode and the gate at ground potential. By varying the cathode voltage, we can acquire data at different drift fields. For the 83m Kr source, we acquired data in the range 194 − 968 V/cm while we probed the interval 80 − 968 V/cm for 37 Ar. The electron extraction field is created between the gate and the positively biased anode and was kept constant at 10 kV/cm by setting the anode to a nominal voltage of +4 kV.
The xenon is continuously purified of electronegative impurities by circulating the gas phase through a hot metal getter, an integral part of the gas handling system. The gas system features a source chamber which allows for the introduction of 83m Kr from the decay of 83 Rb with T 1/2 = 86.2(1) d [15,16]. With an additional gas mixing setup, described in Section 3, that can be connected in parallel to the 83 Rb source chamber, we also deploy a gaseous 37 Ar source.
The S1 and S2 signals are detected by a 2-inch diameter PMT (R6041-06 MOD, Hamamatsu Photonics), placed in the liquid phase at the bottom of the TPC, and sixteen 6 × 6 mm 2 VUV4 SiPMs (2 × 2 array of S13371, Hamamatsu Photonics) placed in the gas phase. To install the SiPMs in the TPC, a PTFE holder which exactly replaced the area occupied by the former Fig. 1 Rendering of the upgraded Xurich II TPC. The active liquid xenon (LXe) is contained within a PTFE cylinder and surrounded by copper field-shaping rings. A 2-inch PMT is placed in the liquid, and an array of 16 SiPMs is placed in the gaseous phase at the top. For the sake of visualisation, two field-shaping rings have been cut. Legend: 1 -PCB with ×10 preamplifier, 2 -SiPM array, 3 -Anode and gate mesh, 4 -Level meters, 5 -PTFE reflector wall, 6 -Weir for liquid level control, 7 -Copper field-shaping rings, 8 -Cathode,  top PMT was designed and built. The SiPMs were arranged in the holder to maximise the light collection efficiency (LCE) as studied in MC simulations [17]. The distance from the surface of the SiPMs to the anode is 5 mm, optimised to obtain good position resolution.
The new configuration of the TPC requires 17 readout channels, 15 more than in the previous design. To handle this number of cables in the small detector, we installed a new CF-40 feedthrough for 20 doubleshielded coaxial signal cables (RG 196) for the photosensors and 7 Kapton ® insulated wires for voltage biasing. To independently read out the 16 SiPM channels, a printed circuit board (PCB) was designed. Each SiPM, on its cathode side, is connected in series with a 10 kΩ resistor to the bias voltage to limit the current and avoid possible damage to the device. On its anode side, each SiPM is grounded with a 50 Ω resistor and coupled to a custom-made ×10 voltage amplifier circuit [17]. The main component of the circuit is the ultra-low noise, non-inverting voltage feedback operational amplifier OPA847 from Texas Instruments. This amplifier has a nominal gain of 20 when connected to an open circuit at the output and is able to reduce the electronic pick-up noise to a negligible level. It operates from direct current to 250 MHz and has a 50 Ω output impedance that matches the input impedance of the ADC. These two resistances form a voltage divider at the amplifier output. Because of this loading effect, the effective final signal amplification with a 50 Ω load from the ADC is tenfold. The heat dissipation of the entire amplifier board is ∼ 3 W, whereas it was negligible for the former 2-inch PMT base. This additional heat load can safely be handled by the cooling system of the detector. Located in the warmer gas phase at ∼ 190 K, the SiPMs have shown an error-weighted mean dark count rate of (8.05±0.03) Hz/mm 2 at a bias voltage of 51.5 V, which is consistent with the temperature behaviour we reported [14]. The dark count rate is the uncorrelated noise rate above a threshold of 0.5 photoelectrons (PE) in scintillation-free data acquired with nitrogen gas in the TPC under very similar thermodynamic conditions to those observed when using xenon.

The 37 Ar calibration source
The 37 Ar isotope decays with T 1/2 = 35.01(2) d via electron capture into stable 37 Cl, with a Q-value of 813.9(2) keV [18]. Depending on the inner atomic shell from which the electron is captured, 2.82 keV, 0.27 keV and 0.0175 keV X-rays and Auger electrons are released (K-, L-, and M-shell, respectively, see Table 1 for the corresponding branching ratios). These offer mono- Ar from 38 Ar and 41 Ar from 40 Ar. Of less relevance are 43 Ar from 42 Ar, of which only traces are present in natural argon, and 42 Ar itself from secondary activation of 41 Ar.
As the natural abundance of 38 Ar is more than five times lower and its cross section for thermal neutrons is more than six times lower, we expect to have produced 39 Ar with an activity lower than 0.23 Bq per ampule. It decays with a half-life of 268 y into its stable daughter 39 K solely via beta decay with an endpoint of 565 keV [23].
The high abundance of 99.6 % of 40 Ar yields a high activity of 41 Ar shortly after the irradiation. This required that the activation happened in two steps with a cool down time in between. 41 Ar decays into stable 41 K with a half-life of 109.6 min [24]. Hence, its activity decreased below 1 Bq after 2 days and was negligible after ∼ 160 days when the source was introduced into the detector. Similarly, if 43 Ar (T 1/2 = 5.4 min) had been produced, it would have been negligible after a short time as the subsequent decay of its daughter isotope 43 K into stable 43 Ca takes place with a halflife of 22.3 h [25]. The secondary activation of 41 Ar to 42 Ar (T 1/2 = 32.9 y), which decays via 42 K into stable 42 Ca [26], is negligible as well since the initial activity can be calculated to be ∼ 50 µBq.

Source introduction setup
To introduce the 37 Ar source into the TPC, the quartz ampules must be placed and broken inside the gas system. We designed and built a dedicated system that contains a mixing chamber and an ampule breaking mechanism similar to the one described in [27] and deployed in XENON1T at the end of its lifetime [28]. We refer to Figure 2 for a schematic view.
The fragile ampule is firmly held in place by a custom-made holder welded onto a CF-40 blank mounted on one of the sides of a CF-Tee. By actuating a guillotine through a vacuum bellow, the neck of the ampule can be cut off in vacuum. Opening a valve to the evacuated mixing volume transfers ∼ 94 % of the argon through a filter. After isolating the ampule chamber again, xenon from a gas bottle can be added until the TPC working pressure of ∼ 2 bar is reached. Bypassing the recirculation flow through the mixing chamber and integrating the flow, the introduced 37 Ar activity can be estimated. The argon gas of multiple ampules can be injected by stopping the recirculation, isolating the TPC and recuperating the xenon/argon mixture remaining in the gas system into a cold trap cooled by liquid nitrogen. The gas system can be pumped after inserting a new ampule from a port on the mixing chamber (cf. Figure 2). The use of the cold trap ensures that the amount of xenon and, hence, the pressure and liquid level remain unchanged during ampule changes. For the data described in the subsequent sections, three ampules with an estimated total activity of 2.7 kBq were broken mid-May 2019 to release the activated argon into the xenon.

Data acquisition
The signal sizes are energy dependent, thus deploying different calibration sources requires an adjustment of the settings in the data acquisition (DAQ). The S1 signals of the 37 Ar lines, for instance, are much smaller than those of the 83m Kr lines and can even be missed by the trigger or absent in the waveform. At the high end in signal size, the PMT channel could, if not limited in voltage, saturate the electronics with S2 signals from 83m Kr. Below, we discuss the DAQ settings and the deployed hardware in detail.
The raw data is digitised by three 8-channel CAEN V1724 modules connected in daisy chain, each with a 100 MHz sampling rate (10 ns samples) and 14-bit resolution at 2.25 V dynamic range. The preamplified SiPM signals are fed directly to the Flash ADCs. The PMT signal is attenuated by a factor of 10 for the 83m Kr data, while not attenuated for 37 Ar, then sent to the digitisers. After passing through a CAEN 625 fan-in/fan-out module the trigger is generated by a CAEN N840 leading edge discriminator based on the PMT signal. For 83m Kr, we require the attenuated PMT signal to exceed 2 mV, which corresponds to a peak with integrated charge of > 15 − 21 PE. For 37 Ar, we trigger on 7 mV of the raw PMT signal corresponding to > 4 − 5 PE. Only for the data at 80 V/cm drift field we use 6 mV. The event window for the 83m Kr data was chosen to be 40 µs with 30 µs post-trigger time, which is sufficiently large to contain the S2 even for the lowest field of 194 V/cm when triggered on the S1. Triggers on 83m Kr S2 signals are expected to be rare and are removed in the analysis phase to avoid a bias of events towards the upper half of the TPC. For 37 Ar, the event window was chosen to be 60 µs with the trigger placed in the middle at 30 µs to contain the entire waveform, regardless of whether the trigger was issued by an S1 or S2. The maximum drift time at the lowest applied field of 80 V/cm inside the TPC volume is ∼ 22 µs.
Depending on the user settings, entire waveforms can be written to disk or, to save memory, only the samples exceeding a certain threshold (henceforth referred to as "good" as distinguished from "skipped") by means of Zero Length Encoding (ZLE), as shown in Figure 3. The first half of the 83m Kr data was taken without ZLE while the rest, including the 37 Ar data, was acquired with ZLE. The ZLE threshold for 83m Kr, which only comes into play if a trigger is issued, was chosen to be 15 bins (∼ 2.06 mV at ADC input) for the PMT channel and 10 bins (∼ 1.37 mV at ADC input) for the SiPM channels. To catch the S1s of 37 Ar at the single PE level, we have chosen 20 bins (∼ 2.75 mV at ADC input) for the PMT channel and 6 bins (∼ 0.82 mV at ADC input) for the SiPM channels as ZLE thresholds. These values were determined empirically and, in the case of 37 Ar, set just above the baseline. In addition to the over-threshold samples, a certain number of backward-and forward-samples can be specified to include the rise and the tail of the peak. To keep timing information, the counts of good samples and skipped samples are recorded as so-called control words [29]. We recorded 150 (200) backward-samples and 150 (100) forward-samples for 83m Kr ( 37 Ar).

Data processing
The data processing consists of three stages. The first is devoted to event alignment and merging of the data files from the three ADC modules into one. The second stage is the main raw data processing of the waveforms on an event-by-event basis. It retrieves basic information about the peaks in the individual channels, identifies pulses as coincidences of these peaks and classifies their pulse type (S1/S2/noise). The third step is the post-processing tailored to the analysis phase: events are built by combining S1s and S2s based on measured charges and physical drift times. Geometry or drift time related corrections are applied later in the analysis, as detailed in Section 5.

Pre-processing
The first stage ensures the correct matching of the events recorded by the three ADC modules by means of a comparison of the trigger time tag changes among triggers. Misalignment among the modules can be caused by an incorrectly propagated busy state. However, an efficient offline realigning-algorithm typically restores 99.9 % of the raw events. Once all the events are aligned, events of one dataset are merged into one ROOT file.

Main processing
Prior to the actual processing, in the ZLE case, the waveforms are reconstructed from the control words. The second stage starts with the baseline calculation of the raw waveforms and the inversion of the negative PMT signal. In the ZLE case, we account for baseline shifts by calculating the baseline for each good-region individually, shown in Figure 3. The threshold for peakdetection (not to be confused with the ZLE threshold for the DAQ) was chosen to be two standard deviations of the baseline distribution around zero. The peak integration limits are dynamically defined based on the variation of the baseline. Integration windows smaller than 3 samples are discarded for noise suppression.
The processor can, to a certain extent, distinguish multiple peaks that overlap. This is necessary to resolve the event topology of 83m Kr with two S2 signals arriving shortly after one another (see Section 5.2) and is helpful to separate afterpulses from the main peak. To reduce the sensitivity to random fluctuations, the method requires over-threshold regions of at least 6 samples. We split two peaks based on their moving average of order of 4 samples (order 2 samples for over-threshold regions ≤ 25 samples) whenever both have fallen below half of their individual maximum value. Thus, for a Gaussian-shaped peak almost 90% of its charge is integrated when cut on one side. The moving average limits the algorithm to a minimum peak separation of 40 ns or 20 ns, respectively. This time-scale is sufficiently short to split even close S1 signals efficiently and is irrelevant for S2 signals for which the requirement of the contained charge dominates.
After a peak is found, separated from others and its integration window fixed, the peak properties are extracted. The area is determined from the summed bin content of each sample, the peak position is defined as the position of the maximum sample. We also determine the width of the peak which serves as discriminator for the S1 and S2 identification. S1 and S2 signals are distinguished by width-based filters, where the filters at the i-th bin are defined as: A j is the baseline-subtracted signal at the j-th bin and w 1 , w 2 are the two widths. While a maximum summation width of w 1 = 20 samples contains the entire S1, w 2 was chosen to be maximum 100 samples. The sums are bounded by the peak integration window to guarantee that it is not summed over a distinct close-by peak. Both filters are evaluated at i being the centre of the full width at half maximum of the peak to account for their asymmetric shape. The S2 filter will be zero for an isolated S1-like signal and much larger for an S2. We choose the ratio S2/S1 as discriminator and identify a peak with S2 if this ratio is > 0.2, i.e. if > 20 % more charge is contained within w 2 than within w 1 . This reveals the importance of good S1 splitting as we rely on the fact that S1s are contained inside w 1 .
For the event building stage, it is essential to form physical pulses from the detected peaks. We define a pulse as a set of peaks of at least two channels in time Fig. 3 Recorded SiPM waveform of a 83m Kr event acquired with ZLE. For visualisation, we present the relevant region of the 40 µs DAQ window of an event from a SiPM channel located at an outer corner of the array that exhibits small S2 signals for events in the horizontal center of the TPC. The peaks are located in two good-regions surrounded by 150 backward-and forward-samples with skipped-regions in between. The integration windows are shaded and enclosed by lines. The baseline, from which the threshold is defined, is calculated individually for each good-region. While for disjointed good-regions, like here, it is based on the first half of the backward-samples, for two consecutive good-regions, the baseline is calculated within the [60, 90] % interval between the two adjacent over-threshold regions. The S2 signals are overlaid with the moving average on which the splitting algorithm is based. While it splits the two S2s because they have both fallen below half of their individual maximum value, it is insensitive to the visible fluctuations of the S2 light. The widths w 1 , w 2 used for the filters are limited to the integration window and, hence, are equal for the S1 signals. In the zoom, it is visible that over-threshold regions below 3 samples are not integrated.
coincidence. The coincidence window is different for S1 and S2 signals and set based on the classification of the PMT peak.

Post-processing
The post-processing script creates high-level variables based on the described output. It selects the two highest-charge S1s and S2s of an event with the corresponding peak properties from those pulses in which the PMT was involved. In addition, it stores the photosensor gains, reconstructs the (x, y)-position of each event and calculates the drift times from the delay between the S1 and S2 signals. We require a positive drift time, i.e. the S1 must happen before the corresponding S2 and their distance must not exceed a maximum value. To reduce the contribution of high-energetic events and ensure correct S1/S2 pairing, we remove saturated signals just below the maximum input voltage of the limiting fan-in/fan-out module of 1.6 V at 11600 ADC bins or higher. The calculation of higher-level variables and the conversion from ADC bins to PE are performed at the analysis stage.

Photosensor gain calibration
The gain of the photosensors was monitored and calibrated about once a week with light from an external blue LED as described in [15]. The gain values are calculated following the model-independent approach proposed in [30]. For the bottom PMT, the measured gain during the 37 Ar data taking phase is (3.76 ± 0.06) × 10 6 at an operating voltage of 940 V. The gain of the SiPMs is found to be uniform across the array and stable during the entire data taking. The average gain of the 16 SiPM readout channels over the data-taking period is (3.12 ± 0.01) × 10 7 at a bias voltage of 51.5 V including the tenfold amplification of the pre-amplifier. The given uncertainties correspond to the standard deviation of the time averaged gains. The stability of the averaged SiPM gain is shown in Figure 4.

Position reconstruction and fiducial volume
The top photosensor array allows for the reconstruction of the (x, y)-position of interactions. We use a simple centre-of-gravity algorithm to determine the (x, y)- Fig. 4 Error-weighted mean gain (tenfold amplified) of the 16 SiPM channels over the entire data taking period in 2019. We measured (3.12 ± 0.01) × 10 7 . The ±1σ and ±2σ uncertainty bands around the mean are shown in red. Data with 83m Kr ( 37 Ar) was acquired before (after) the time-axis break.
position of an event from the light distribution among the photosensors caused by the largest S2 signal. We denote the physical position of the i-th SiPM in the array by (X, Y ) i , with the origin of the coordinate system being the centre-of-gravity of the array. We then obtain the (uncorrected) event position in x from a geometryand gain-weighted sum of the integrated S2 charge of each sensor 1 , where the normalisation Q tot S2 := 16 i=1 Q i S2 /g i is the total collected S2 charge in the top array with the gain g i of the i-th SiPM channel. While this approach neglects the fact that the light intensity follows an inverse-square law with the distance from the light production site, we find that the algorithm is an adequate choice for our small-scale TPC and leads to excellent reconstruction in the TPC centre. Because of the square arrangement of the SiPMs and the nature of the centre-of-gravity approach, the resulting positions will not represent the actual circular cross-section of the TPC. In fact, while the reconstruction in the TPC centre is very good, events at the TPC boundary feature a square-shaped bias. For this reason, the square is scaled down to unit side length with a constant, x scal = x/c, and then mapped onto a unit circle that is scaled back, While these mapped positions are centred around the origin of the coordinate system due to the gain correction in Equation 2, they still feature a radial bias because of the solid angle seen by the individual photosensors. Using a data-driven approach, we correct for this effect by comparing the positions of the spots visible in the (x, y)-distribution with the known gate mesh junctions. In Figure 5, left, we show the (x, y)distribution from 37 Ar data overlaid by the drawing of the meshes. The spots occur due to focusing of drifting electrons to the high-field regions around the junctions of the grounded hexagonal gate mesh. This limits the accuracy of the (x, y) event position reconstruction to roughly half the distance between the junctions, which is ∼ 1.5 mm.
Up to a radius of ∼ 10 mm from the TPC centre, the correction function is fairly linear, i.e. the correction can be performed with a scaling factor, whereas for larger radii it diverges. This can be described by a projection of a sphere onto a plane as: where r map = x 2 map + y 2 map is the mapped radius of the interaction site with r map,max being the boundary of the TPC after the mapping and d TPC = 31 mm the true diameter of the TPC. The linear region of this projection defines our radial fiducial volume cut to be r < 10 mm. The z-position of an interaction is obtained from the time difference between the S1 and S2 signals. We extract the electron drift speed for a given electric field from the clearly visible position of the cathode and the gate in the drift time histogram and their spatial separation. The procedure, together with a drift speed measurement at different fields in Xurich II, are detailed in [15]. The drift speeds that we find here are in agreement with this former measurement and are in the range 1.5 − 2.0 mm/µs.
To avoid high and potentially non-uniform electric field regions around the electrodes, we select only events within z ∈ [−29, −2] mm, where 0 mm is the position of the gate and −31 mm is the position of the cathode. The fiducial region defined by these boundaries is shown in Figure 5 and contains a xenon mass of ∼ 24.5 g.

Data analysis
In this section, we describe the data selection, the applied quality cuts and the signal corrections. We require correctly paired S1 and S2 events for both calibration sources. The focus is on the 37 Ar source, since it is a rather new approach to calibrate the low-energy region, while calibrations with 83m Kr are well established. In this work, we show 37 Ar data from the K-shell capture, with a line at 2.82 keV only.

Data selection and quality cuts
The data used for the calibration was taken between late-May and early-July 2019 while the campaign for the half-life measurement lasted until mid-October 2019.
As described at the end of Section 4.2, the postprocessing accounts for physical drift times and removes saturated signals. Because of the low 37 Ar signal rate in the detector volume, O(10 Hz), we expect a negligible pile-up rate for a DAQ window of 60 µs. This includes pile-up with background events with a rate of the same order of magnitude. Furthermore, we select only single S2s, for we do not expect multiple scatters from a single 37 Ar decay. To ensure correct S1 and S2 matching, we require the S2 signal to contain more charge than the S1. In addition, we apply the fiducial volume cut described in Section 4.4 to remove wall and gas events, regions with an inhomogeneous electric field, and those outside of the linear region of the (x, y)-position reconstruction algorithm. To efficiently remove events in the gas phase that yield a large S2 in the top array compared to the total S2, we apply a so-called area fraction top cut shown in Figure 6, left. To remove the lower energy lines from the analysis, we apply a cut on the energy, namely in the centre of the 5σ boundaries of the L-and the K-shell populations in the total S1+S2 space. Since accidental coincidences have by definition a ran-dom drift time, they can be efficiently removed with a cut on a variable with a strong drift time dependence. The longitudinal diffusion of the electron cloud during the drift leads to such a dependence of the S2 width, as shown in Figure 6, right. This distribution follows an empirical function with three parameters: Here, the width corresponds to the full width at tenth maximum. The upper and lower cuts are defined as follows: the histogram is sliced in drift time bins, and the fit is performed through the 97.7 % and 2.3 % quantiles of those bins. We thus select the ±2σ region around the mean. For 37 Ar data, 5.6 % of the recorded events that generate an S2 signal pass all data quality cuts.

Corrections
As a result of the recombination of the drifting free electrons with electronegative impurities in LXe, the S2 charge degrades with drift time. The free electron lifetime τ e , defined as the decay constant in an exponential decay law, is a measure of the LXe purity. Being a property of the LXe, the drift time dependence of the S2 signal of any channel can be used to correct for this systematic effect. We base our correction on the bottom PMT as it collects more S2 light than the top array. To this end, we fit the 50 % quantiles of the 2-dimensional (a) Area fraction top cut to remove gas events with a fraction of 0.3−0.5 . Besides the K-shell population on the right we also select parts of the L-shell population below 600 PE S2b and the single photoelectron region that are both mostly located in S2-only space. These are removed by an energy cut.
(b) S2 width cut to remove accidental coincidences and remaining gas events (see text). The drift region defined by the cathode and the gate are marked by solid black lines. S2b 2 versus drift time histogram with an exponentially decreasing function (see Figure 7, left). We thus obtain the corrected S2 charge by multiplying with its inverse, i.e. we scale according to the S2 at zero drift time: The stable detector conditions during the time in which the calibration data was taken resulted in an electron lifetime of 85 − 187 µs without showing any abrupt changes. While we correct for the electron lifetime on a daily basis, in Figure 7, left, we show for simplicity the combined data at 968 V/cm drift field with an averaged electron life time. The LCE depends on the geometry and the optical properties of the TPC. In particular, the collected S1 light varies with the depth of the event. We use a data-driven approach and fit the 50 % quantiles of the 2-dimensional S1b versus drift time histogram with a 4th order polynomial, as shown in Figure 7, right. The corrected signal cS1b is obtained by scaling around the drift time corresponding to the centre of the TPC. Thanks to the excellent single photoelectron resolution of SiPMs, the S1t histogram features a band-like structure, as visible in Figure 8. Since the small light signal seen by the SiPMs does not show any drift time dependency within a single band, no data-based drifttime dependent correction could be performed for S1t. Because of this discreteness, we would need to know which events changed to the next lower band after a certain drift time. However, these bands help us to fix the charge scale in PE. To that end, Gaussian fits on the seven lowest bands are performed in the S1t histogram and then linearly mapped onto the expected charge in PE. Figure 8 shows the corrected result. All the described corrections are performed on data within the fiducial volume only. For simplicity, we do not apply (x, y)-dependent corrections.

Energy calibration with 83m Kr
A common technique for the energy calibration of a xenon dual-phase TPC is to introduce 83m Kr into the system [4,[31][32][33]. This source with T 1/2 = 1.83(2) h offers a homogeneous distribution of events in the target volume with a line at 32.15 keV followed by a line at 9.41 keV with T 1/2 = 155.1(12) ns for the intermediate state [16].

Data selection and quality cuts
The selection of 83m Kr events is performed by exploiting the double-S1 and double-S2 signal topology as a result of the intermediate 9.41 keV level. In the following, we only consider split S1 and S2-pairs and not the combined 41.56 keV line. That is, in the PMT channel, we require exactly two S1s with S2s and the S1s must contain less charge than their paired S2s. In a next step, a cut on the delay time of the PMT signals is applied. We expect 83m Kr events to have the same positive delay ∆t between the two S1 signals as between the two S2 signals and select ±2σ around the centred mean in the ∆t S1 − ∆t S2 histogram. As for 37 Ar, we use the fiducial volume cut from Section 4.4, an area fraction (a) Electron lifetime correction fit with relative residuals. The electron lifetime from the combined data at 968 V/cm is (124.2 ± 0.7) µs. The error is the error on the fit parameter and purely statistical.
(b) S1b correction fit with relative residuals.  top cut, an energy cut to separate the lines, as well as the S2 width cut. For 83m Kr data, 1.9 % of the recorded events that generate an S2 signal pass all data quality cuts.

Corrections
We apply the same corrections on the 83m Kr lines as described above for 37 Ar. Due to the higher energy of the source, and thus the higher S1 signal seen by the top SiPMs, we are able to correct S1t with a 4th order polynomial showing qualitatively the opposite drift time dependency compared to the S1b signal (cf. Figure 7, right).

Results and Discussion
After applying the described corrections, we show in Figure 9 the localised populations from the 2.82 keV line of 37 Ar (left) together with the two lines of 83m Kr (right) in total S2 versus total S1 space. The anticorrelation between ionisation and scintillation signals is evident, however less pronounced at the low energy of 37 Ar due to the small S1 signals. From the mean of the 37 Ar distribution, we deduce a light yield of (4.985 ± 0.003) PE/keV and a charge yield of (1504.0 ± 0.2) PE/keV for 2.82 keV at 968 V/cm drift field. For 83m Kr, at the same drift field, we obtain a light yield of (7.941 ± 0.003) PE/keV and a charge yield of (702.0 ± 0.3) PE/keV at 9.4 keV, and (6.284 ± 0.002) PE/keV and (1080.1 ± 0.3) PE/keV, respectively at 32.1 keV. The uncertainties reflect the statistical errors on the mean of the yields. The discussion of the systematics, in particular the limitation of the 83m Kr charge yield precision, follows below. In Figure 10, we present the 2.82 keV population of 37 Ar for the S2 light of the top array only. From the comparison with Figure 9 (left), we deduce that 23.5% of the total detected S2 light is collected by the SiPMs.
Calculating the total light and charge yield for data at different drift fields yields the so-called Doke plot [34], shown in Figure 11. The light and charge yield measurements of the three lines from 37 Ar and 83m Kr feature a clear anti-correlated linear dependence. The 37 Ar points are less spread than those from 83m Kr because the drift field dependence of the ionisation and scintillation yield decreases rapidly below ∼ 10 keV ER energy [35]. The errors include statistical uncertainties,   as well as systematic uncertainties of the PMT gain and of the peak-splitting method described in Section 4.2. The condition on the contained charge explains the systematic bias of the 32.15 keV points towards lower charge yield and is the origin of the asymmetric vertical error bars assigned to 83m Kr. However, this is not relevant for the light yield as the S1 delay time is sufficiently long to provide complete S1 separation for any split S2 signal.
The detector-specific gains g1 and g2 are defined as g1 := S1/N γ and g2 := S2/N e − , where N γ is the total number of photons produced by de-excitation or recombination and N e − is the total number of electrons producing the S2 signal. With W = (13.7 ± 0.2) eV being the mean energy required to produce an excited or ionised xenon atom [36], the combined energy scale of an event is where we substituted the number of quanta for the definitions in Equation 7. From the axis intersections of the linear regression in Figure 11, we obtain g1 = (0.150 ± 0.014) PE/γ and g2 = (36.9 ± 2.1) PE/e − . The combined energy resolution σ/E from S1 and S2 signals of the 2.82 keV 37 Ar line is (17.2 ± 0.1) % at 968 V/cm and stays constant down to 600 V/cm. Below this value, the resolution degrades slowly to (18.0 ± 0.1) % at the lowest measured field. The resolution at the highest energy of 32.15 keV is (5.80 ± 0.01) % which is consistent with the one we previously measured with the TPC with two PMTs [15]. The yields of all lines at 968 V/cm are compared to the Noble Element Simulation Technique (NEST) v2.0 [37] model predictions in Figure 12. Both charge and light yield at 2.82 keV and 32.15 keV agree well with the gamma model. In addition, we show the predictions for the two lines of 83m Kr.
While we see good agreement at 32.15 keV, we find a lower charge yield than predicted at 9.41 keV. Due to anti-correlation, the light yield at the latter energy is higher than predicted but consistent with NEST v2.0 83m Kr within the error. Data with the 37 Ar source was acquired over a period of more than 4 half-lives. Selecting the data as described in Section 5 in daily slices, we show in Figure 13 the evolution of the activity of the 2.82 keV line. For simplicity, we assume that all the systematics and applied cuts are stable over the time of the measurement. The exponential fit yields a half-life of (35.7 ± 0.4) d, where the most recent literature value is (35.011 ± 0.019) d from a weighted average of four measurements [18]. The given error is the error of the decay parameter of the least squares fit and does not include systematic uncertainties. However, within the scope of this work, we can be confident that we indeed selected events from 37 Ar decays.
For 83m Kr, exponential fits of the S1 and S2 delay histograms yield (161 ± 5) ns for the half-life of the intermediate state, in agreement with the literature value stated in Section 5.2. Fig. 13 Activity evolution of the selected 2.82 keV data of 37 Ar over more than 4 half-lives with a one-day binning starting from the time at which the first ampule was broken. The fit was performed on 968 V/cm data only, because at other drift fields the data features systematic deviations caused by the varying efficiency of the applied data-driven cuts and corrections. The errors on the individual data points are Poissonlike. We obtain a half-life of (35.7 ± 0.4) d.

Summary
In this work, we have presented for the first time a dualphase xenon TPC equipped with SiPMs in the top photosensor array. We operated the detector under stable conditions for more than one year and have demonstrated its performance with data from internal 83m Kr and 37 Ar calibration sources. Based on the centre of gravity of ionisation signals detected by the SiPM top array, we have reconstructed the (x, y)-position with a resolution of 1.5 mm within the fiducial radius. The low energy threshold of the detector and its energy resolution have allowed for an S1 and S2 analysis of the 2.82 keV line of 37 Ar. The light and charge yields were compared with NEST v2.0 predictions, and provide valuable input for the software. Ongoing studies include an S2-only analysis of the clearly visible 0.27 keV line and an attempt to look for the lower energetic Mshell capture at 17.5 eV. Our study has demonstrated that SiPMs are indeed a promising alternative to PMTs in small detectors. The scaleability to much larger arrays, such as required by the DARWIN multi-ton liquid xenon detector [11] is under investigation. A major challenge, apart from reducing the dark count rate, is the readout and the signal pre-amplification. The main shortcoming of the described amplifier board for 16 channels is its heat dissipation of ∼ 3 W. To that end, we are currently investigating the lowest allowable pre-amplification factor and testing various low-power pre-amplifiers at cryogenic temperatures.