Giovanni Pireddu

Irina Svir

Christian Amatore

Alexander Oleinick

Modelling Diffusion at Random Arrays of Electrodes: Revisiting the Voronoi Tessellation Concept

Keywords: Random array of electrodes, Voronoi tessellation, amperometry, modelling, simulations

In order to assess the precision of simulations of diffusion at arrays based on the Voronoi tessellation approach we investigated two representative types of random arrays involving bands or disk electroactive sites [J. Electroanal. Chem. 147 (1983) 39] and modelled their diffusional patterns when the solution contains only one electroactive species undergoing a simple electron transfer reaction at the active sites under chronoamperometric conditions. On the one hand, the ensuing results establish that in both cases the Voronoi approach produces reasonably good predictions of the total current intensities flowing through the elementary cells of each array. Indeed, the relative error introduced by the Voronoi-based approach, being less than 5% in each case, is acceptable from an experimental point of view owing to many other sources of uncertainties involved for random arrays. On the other hand, this work demonstrates that the current intensities predicted for the individual sites within an elementary cell based on simulations using Voronoi approaches are excessively wrong since they totally neglect the excessive redistributions of flux lines that happen when the diffusion layers expand with time. This hints to possible severe complication when Voronoi approaches are applied to encompass complex kinetics or to predict the outcome of electroanalytical methods that rely on a fine coupling between local concentrations and local fluxes.

Introduction

Many practical micro-and nanoscale systems used presently in heterogeneous catalysis, electrolyzers or electroanalytical sensing can be represented as arrays of active sites distributed randomly on inert surfaces [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Amatore | Electrochemistry at Ultramicroelectrodes[END_REF][START_REF] Nesselberger | The effect of particle proximity on the oxygen reduction rate of size-selected platinum clusters[END_REF][START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF] since this allows to decrease the amount of expensive materials and often provides even better performances than solid materials owing to enhanced micro-and nanoscale diffusion rates and/or more reactive surface atoms. Those may include partially blocked electrodes, membranes or non-(electro)chemically active substrates with active/catalytic nanosites dispersed on their surface to name a few. As shown previously, these systems can be efficiently addressed theoretically by using a mathematical concept introduced by Voronoi [START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Aurenhammer | Voronoi Diagrams -A Surver of a Fundamental Geometric Data Structure[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF]. Indeed, taking advantage of the local symmetries between pairs of active sites scattered on an inert surface, one may adapt the Voronoi tessellation principle to define independent juxtaposed diffusional domains (named unit cells hereafter) around each active site whose full paving of the substrate and of the solution permits evaluating the global current responses of random arrays [START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF].

Briefly, these approaches consist in identifying around each active site a series of planar surfaces perpendicular to the inert surface positioned in such a way that each point of such plane is equidistant from the considered site and an immediately neighboring one. Thus, a collection of juxtaposed polygonal prisms constructed over each active site paves the whole diffusion domain accessible to the global array. This evidently facilitates the treatment of the random array by allowing to perform a series of parallel independent calculations of the diffusional contribution due to each site through considering individually its unit Voronoi cell followed by a simple summation of their individual current intensities at any given time [START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF][START_REF] Oleinick | Reconstruction of Distributions of Nanoparticles or Electroactive Nano-Components in Electrochemical Arrays Based on Chronoamperometric Data[END_REF]. Interestingly, each of these unit Voronoi cells can be additionally approximated by volumes of simplified regular geometries such as cylinders or blades [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Amatore | Electrochemistry at Ultramicroelectrodes[END_REF][START_REF] Sliusarenko | Validating a central approximation in theories of regular electrode electrochemical arrays of various common geometries[END_REF] with the further huge benefit of considerably decreasing the simulations difficulties, hence their time durations and memory occupancies [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Amatore | Electrochemistry at Ultramicroelectrodes[END_REF][START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF] and even possibly allowing analytical empirical formulations [START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF][START_REF] Oleinick | Reconstruction of Distributions of Nanoparticles or Electroactive Nano-Components in Electrochemical Arrays Based on Chronoamperometric Data[END_REF].

This explains why such approaches based on Voronoi tessellation concept is tempting owing to its simplicity and computational efficiency [START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF] although they evidently idealize the real diffusional contributions to and from each active site. Nonetheless, in a series of previous works we established that they allow evaluating with a reasonable accuracy the arrays global currents under chronoamperometric conditions provided that the substrate concentrations at the surface of each site are identical [START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF][START_REF] Oleinick | Reconstruction of Distributions of Nanoparticles or Electroactive Nano-Components in Electrochemical Arrays Based on Chronoamperometric Data[END_REF]. Hence, in this work we wish to examine this issue in more detail, particularly through investigating if the Voronoi tessellation constructed on the surface of the system remains valid at distances from it that exceed the average distance between active sites. For this we will rely on modelling the diffusional patterns created over two typical ideal arrays consisting of band or disk electrodes [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Amatore | Electrochemistry at Ultramicroelectrodes[END_REF] since they will adequately serve our purpose of evidencing the qualities and the limits of the Voronoi approaches when applied to more genuine experimental situations.

Computational details

Simulations were performed with Comsol Multiphysics 5.5 [START_REF]COMSOL Multiphysics Reference Manual, version 5[END_REF] to solve the system of partial derivative equation and boundary conditions described in Section A of the Supplementary Material (SM). However, in the case of the band array, analytical derivations based on conformal mapping (see Section B of the SM) were performed in order to check the accuracy of Comsol results and to identify more accurately the trajectories of flux lines in the elementary cell (see Figure 3 and Figure SM2). For global simulations of the disk array elementary cell behavior, the flux lines (see Figure 7 or Figure SM2) were obtained through using a specific feature of Comsol Multiphysics 5.5.

Results and Discussion

Hereafter, we consider two representative systems involving active sites of identical widths (array of parallel bands) or radii (array of disks) pseudo-randomly embedded on a conductive non-electroactive flat surface. Each array is assumed to operate under amperometric conditions in a solution containing a single electroactive species, A, prone to undergo a single electron transfer (ET) at each active site surface and diffusion being assumed to be a single transport mode. To simplify the presentation and focus on the important outcomes, the physico-mathematical models used to derive the theoretical responses analyzed hereafter are fully described in section A of the SM.

Pseudo random array composed of parallel band active sites

To simplify the following analysis and the discussion of its results we envision a pseudo-random array consisting of a periodic infinite lattice composed of pairs of parallel bands arranged asymmetrically (i.e., in a mirror-like fashion) with a period P like the one whose cross section is shown in Figure 1a.

Considering an infinite array extending to infinity in each direction allows to neglect any edge effect. Therefore, the symmetry planes between two adjacent couples of bands whose cross-sections are represented by the dashed lines perpendicular to the array surface in Figure 1a delimit individual array components that account for any other one (after a mirror inversion as needed). Each of these components, termed elementary cells hereafter, may then be treated independently to predict the whole array current response after a simple scaling taking into account the number of elementary cells.

To evaluate the performance and possible problems associated with the Voronoi tessellation model represented on Figure 1c, we decided to perform two series of simulations, both using with Comsol ® to solve the proper set of partial derivative equations and boundary conditions (see section A of the SM), either (i) within the whole solution area associated to one elementary cell without considering any physical limit inside it as sketched in Figure 1b, or (ii) independently within each Voronoi unit cell composing the elementary one as shown in Figure 1c. This allowed comparing the two series of results as a function of time in terms of (a) the global elementary cell currents or of (b) those passing through each electrode. For these comparisons we decided to opt for representations of the results on log-log plots since they allow identifying the relevant limiting regimes much better than classical chronoamperograms [START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF].

Figures 2a,b compare the results of simulations performed for (Fig. 2a) the whole elementary cell (shown in Fig. 1b) to those (Fig. 2b) obtained based on two independent simulations performed in each of the Voronoi unit cells displayed in (shown in Fig. 1c), as well as to the different limiting Cottrellian behaviors expected in each case at short and long times. A first conclusion is that the two simulations produced almost identical time-dependence behaviors for the overall currents (green curves in Figs 2a,b). Both reached the expected Cottrell limits at short time (proportional to the sum of each band surface area) when the diffusion layers thicknesses are much smaller than the band widths (i.e., Dt << w 2 [START_REF] Deakin | Electrochemical kinetics at microelectrodes: Part II. Cyclic voltammetry at band electrodes[END_REF]) and long time when the two individual diffusion layers generated at each band fully overlap to generate a planar diffusion layer (in the elementary cell case, Fig. 2a) or when each diffusion layer developing in each Voronoi unit cell became planar (Fig. 2b). However, since the diffusion layers transition towards their planar limits [1-2, 6, 8] in each cell (e.g., the elementary one in Fig. 1b, and the two independent ones in Fig. 1c) did not occur exactly in phase owing to the different ratios between the electrode(s) and solution domain widths in each case, the two global currents exhibited a small discrepancy (at maximum close to 2%, see Figure 2c) over the corresponding time range.

Although such discrepancy between the predictions of two models (Fig. 1b vs. 1c) can be considered tolerable given many other uncertainties related to the position of individual sites in any random array it highlighted as a possible issue of concern related to a systematic use of approaches based on the

Voronoi tessellation concept (see below).

A simple comparison of Figures 2a and2b relating to the individual currents flowing through each electrode according to global or Voronoi based simulations perfectly documents the origin and the importance of this point. Indeed, the two sets of predicted variations differ significantly at moderate and long times as can be seen clearly by comparing the respective Cottrellian limits achieved at large times for each band electrode according to the simulation procedure used (compare the blue and orange curves relative to each case in Figures 2a and2b). These large differences arise from the fact that the solution domains delimiting the two concentration fields generated by each band within the elementary double electrode cell varied with time, which cannot be explained by the constant solution domains corresponding to a Voronoi tessellation pattern. These deviations are negligible at short times because the diffusion layers are then limited to the environment close to each electrode, and they exactly compensate at long times since the global planar diffusion layers which progress into the solution far from the surface of the array are identical in each case. However, the two sets of concentration fields are not progressing in phase towards the solution bulk at intermediate times.

Therefore, a mutual compensation of the error associated with the independent Voronoi-based currents cannot occur as observed at long times. Nonetheless, a maximal error of only ca. 3% is obtained for the predictions based on the Voronoi approach relative to the more exact ones performed for the whole elementary cell (Fig. 2c).

This singular behavior becomes definitely apparent by observing the patterns described by the flux lines issued from each opposite edge of adjacent electrodes at different experimental times as done in Figures 3a-c. To simplify the presentation of these results and avoid complex but less relevant issues we decided to modulate the extension of the global diffusion layer by imposing a convection-like limit at  conv = 10, 40 and 60 µm respectively instead of using a classical  diff (t) = (Dt) 1/2 time-dependent expansion (note that the boundaries of the diffusion layers shown in Figure 3a-c would correspond to classical chronoamperometric experiments of durations of about 10, 160 and 360 ms, respectively, for D = 10 -5 cm 2 s 1 ). Indeed, this approach presented the great advantage of eliminating the effects of residual curvatures around the end of diffusion layers [START_REF] Amatore | The real meaning of Nernst's steady diffusion layer concept under non-forced hydrodynamic conditions. A simple model based on Levich's seminal view of convection[END_REF][START_REF] Pebay | Mass transport at infinite regular arrays of microband electrodes submitted to natural convection: theory and experiments[END_REF][START_REF] Amatore | Monitoring concentration profiles in situ with an ultramicroelectrode probe[END_REF][START_REF] Amatore | Cyclic voltammetry at microelectrodes. Influence of natural convection on diffusion layers as characterized by in situ mapping of concentration profiles[END_REF].

Figure 3 clearly evidences that upon increasing the global diffusion layer thickness the left electrode is progressively "eating" its way across the solution domain that pertains to the right one within a Voronoi framework. In other words, the current density flowing through the left electrode increases with time compared to its Voronoi-predicted time-dependence. The opposite evidently occurs for the right electrode so the two phenomena almost exactly compensate at short or long times (Figure 2c). However, one notice in Fig. 3b that for  conv = 40 µm the concentration fields produced by the left electrode has almost reached a planar behavior near the bulk solution, while that generated by the right electrode still displays marked cylindrical characteristics [START_REF] Deakin | Electrochemical kinetics at microelectrodes: Part II. Cyclic voltammetry at band electrodes[END_REF].

Interestingly, at long diffusion times the concentration gradients near the end of the diffusion layer become linear with the same slope independently of the abscissa. Hence, the lateral driving force that curves the flux lines vanishes so these later ones become perpendicular to the array surface although at a position different from that predicted by the Voronoi approximation. Consequently, the ratio of the current intensities flowing through each electrode becomes constant and each electrode may reach a Cottrellian limit as evidenced in Figure 2a. Nonetheless, these limits differ from those predicted upon using a Voronoi-based approach because they correspond to the surface areas delimited by the edge flux lines at infinite distance from the electrodes (compare Fig. 3c) and not to those delimited by the Voronoi-tessellation of the array surface.

Everything appears as if the diffusional domains featuring the range of action of each band at long times corresponded to unit cells with surface areas that differ from those retained within a Voronoitessellation framework. Figure 4 emphasizes this outcome in a compact form for different positions of the midpoint between the facing edges of the two band electrodes (x m values), by plotting the relative magnitudes of the virtual-surface areas sustaining the action of each band at long times. Since at long times the current at each electrode is proportional to these virtual-surface areas Fig. 4 reflects also relative magnitudes of the individual electrode currents and could thus be measured experimentally if the band electrodes could be connected independently.

The above results (Figs. 3,4) concern systems composed of identical band electrodes, however the main conclusions remain valid even when the bands have different widths as reported in the Supplementary Material (Figure SM2). The important outcome is that the Voronoi approach provides in all cases reasonably correct predictions of the total current time-dependence when the concentration of the electroactive species at each active band surface is the same. Then, Voronoi tessellation approaches offer a great benefit by considerably decreasing simulation costs associated with predicting the behavior of arrays involving much higher numbers of electrodes. This advantage was employed for investigating the global current responses of random nanodisk arrays of electrodes [START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF] and even used for identifying the type of array at hand based on chronoamperometric currents [START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Oleinick | Reconstruction of Distributions of Nanoparticles or Electroactive Nano-Components in Electrochemical Arrays Based on Chronoamperometric Data[END_REF]. However, the present analyses established that Voronoi-based predictions concerning the individual band electrodes current responses are not acceptable at intermediate and long times because the diffusion fields developing far from the electrode network surface do not respect the limits imposed by a Voronoi tessellation framework. Although the ensuing consequences of this intrinsic prejudice are trifling under the conditions investigated here, i.e., when the concentration of the electroactive species is identical at the surface of each electrode, this issue suggests the existence of a potential problem whenever this constraint cannot be maintained as, for example, when ohmic drop or migration phenomena are involved [START_REF] Amatore | Electrochemical kinetics at microelectrodes: Part IV. Electrochemistry in media of low ionic strength[END_REF][START_REF] Amatore | Capacitive and Solution Resistance Effects on Voltammetric Responses at a Disk Microelectrode Covered with a Self-Assembled Monolayer in the Presence of Electron Hopping[END_REF] or when electrodes surface concentrations depend on the local current flux as occurs for slow or quasi-reversible electron transfers or many more complex mechanisms.

Random arrays with disk-like active sites

Random arrays of disk-type active sites are more commonly used in actual experiments so they have been the focus of most theoretical works (see e.g. [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Davies | The cyclic and linear sweep voltammetry of regular and random arrays of microdisc electrodes: Theory[END_REF][START_REF] Belding | Diffusional Cyclic Voltammetry at Electrodes Modified with Random Distributions of Electrocatalytic Nanoparticles: Theory[END_REF][START_REF] Oleinick | Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization[END_REF][START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF]). Therefore, the electrochemical behavior of an infinite pseudo-random 2D-array of disk-electrodes whose elementary cell is represented in Figure 5b was investigated. One may construct the full array through an infinite series of recursive mirror-symmetries performed along each side of the squared elementary cell as evidenced in Fig. 5a. Indeed, because of the geometrically exact periodical mirror-symmetries, the 4 planes elevated over each of the elementary cell sides perpendicularly to the array surface into the solution define a squared-based prism in which zero-flux boundary conditions strictly apply at each vertical face. Note that this conclusion is totally independent of any Voronoi-type consideration because it is an intrinsic property of the symmetries and anti-symmetries that define the full array in Figure 5a. Furthermore, since we consider that the array in Fig. 5 is infinite any edge effects are negligible. Hence, the whole array performance can be rigorously anticipated by a simple scaling of the current responses predicted for the elementary cell in Fig. 5b equipped with its vertical insulating walls without the need to consider any of its adjacent elementary cells.

Figure 6a presents the time dependence of the simulated total current response for the whole elementary cell when all three electrodes 1 to 3 are simultaneously active under chronoamperometric conditions, as well as the time-variations of the current intensities evaluated at each of the electrodes.

Conversely, Figure 6b presents the time dependent simulated currents obtained under the same conditions as in Fig. 6a for each Voronoi unit cell 1 to 3 as well as their sum. Comparison of the two sets of data establishes that provided that the concentration of the electroactive species is identical on all disks surfaces the global currents for the elementary cell evaluated through each approach are almost exactly superimposable except over the [10 s, 100 s] time range in which the Voronoi approach yield a maximal error close to 4% relative to more correct simulations considering the elementary cell as a whole (Fig. 6c). Note that over this time range, the three-discs system has almost already reached its longtime limit [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF]. This validates the use of the Voronoi concept when applied to disk electrodes for chronoamperometric conditions whenever this range of error is acceptable. Again, this may only occur because of almost complete mutual compensations of the errors associated at intermediate times to the current intensities evaluated independently for each disc electrode in a Voronoi approach context (see below and compare Fig. 2a-c and Fig. 3).

In strong contrast with this excellent performance of the Voronoi framework for predicting the global time-dependent current intensity of the whole array, the chronoamperometric behavior predicted for each electrode through simulation of the elementary cell as a whole was very different from that simulated in each Voronoi unit cell except at very short times (t < 0.2 s) when they did not yet interact (compare Figures 6a and6b). As observed above for the band electrodes array, at long experimental times, each individual disk reaches a Cottrellian current limit but corresponding to a virtual surface area extremely different from that predicted for the same electrode within a Voronoi-tessellation framework (Fig 5c). In Fig 5c, the surface areas of the Voronoi unit cells 1-3 cross-sections represent respectively 17%, 29% and 54% of the global surface area of the elementary cell. Accordingly [START_REF] Amatore | Charge transfer at partially blocked surfaces. A model for the case of microscopic active and inactive sites[END_REF][START_REF] Amatore | Construction of optimal quasi-conformal mappings for the 2D numerical simulation of diffusion at microelectrodes.: Part 2. Application to recessed or protruding electrodes and their arrays[END_REF][START_REF] Amatore | Numerical simulation of diffusion processes at recessed disk microelectrode arrays using the quasi-conformal mapping approach[END_REF], the long time Cottrellian limits of the three currents obey the same proportions when predicted using a Voronoi framework (Figure 6b). Contrariwise, when the diffusion problem is treated globally within the whole elementary cell, the cross-sections of the diffusion volumes developing from electrodes 1 and 2 located in the smallest Voronoi areas tend to expand well beyond the limits defined by the Voronoi tessellation of the array basal surface with a tendency to equalize the current densities at each electrode. This is readily noticeable in Figure 6a and Table 1 (compare 1 st and 2 nd columns) since long-time Cottrellian limits of the three electrode current intensities equal approximately one third of the global one, viz. contributing to 29%, 34% and 37%.

The origin of this behavior is perfectly understood upon visualizing the development of the diffusion patterns as in Figure 7 by representing the flux lines issued from electrodes 1 and 3 at long time (those issued from electrode 2 are omitted for clarity of the 3D representation). Would the Voronoi approach be valid, each group of flux lines ought to remain confined within the Voronoi boundaries of each unit cells 1 and 3 featured by the green dashed polygon drawn on the bottom and top cross-section of the prismatic diffusion volume of the elementary cell. Contrariwise, one observes that although all flux lines obey this prediction at small distances from the array surface, they rapidly bend farther into the solution, so that those issued from electrode 1 (i.e., the disk electrode associated with the smallest Voronoi domain in Fig. 5c) promptly expand to "invade" the solution volume "reserved" to electrode 3 within a Voronoi-based framework. Conversely, those issued from electrode 3 avoid the solution ranges "invaded" by electrode 1 flux lines. Although not shown for clarity on Figure 7 the same occur for electrode 2 whose flux lines are "compressed" by those of electrode 1 while simultaneously expanding into the Voronoi domain pertaining to electrode 3. Hence, diffusion leads to a progressive homogenization of the cross-section surface areas of the solution domains sustaining the activity of the three electrodes. Accordingly, the three current densities tend to become equal while the experimental time increases. This phenomenology is similar to that noted above in Figure 3 for band arrays (see also Figure SM2 in SM). Nonetheless, in the present case as for the random band array (compare Figs. 2a,b) even if the Voronoi approach does not represent correctly the individual current intensities at each electrode, the global current is almost correctly predicted over the time range that matters [START_REF] Sliusarenko | Development and validation of an analytical model for predicting chronoamperometric responses of random arrays of micro-and nanodisk electrodes[END_REF] because of self-compensation of errors. Indeed, at long diffusion times, flux lines that do not hit the electrode predicted according to the Voronoi-based approach end up anyway onto an adjacent electrode. Hence, the global current intensity is not significantly altered. However, we must point out that this may occur only because diffusion is essentially a zero-order process, since pure-diffusion second Fick's laws are linear and chronoamperometric boundary conditions under ideal conditions on the plateau of a wave are also zero-order for a Nernstian electron transfer. It is expected that this may not be the case when the boundaries at the electrodes surface involve first-order or higher-order conditions coupling fluxes and surfaces concentrations either because of electron transfer kinetics or ohmic drop, or when the Fick's laws involve nonlinear contributions such as migrational transport [START_REF] Amatore | Electrochemical kinetics at microelectrodes: Part IV. Electrochemistry in media of low ionic strength[END_REF][START_REF] Amatore | Capacitive and Solution Resistance Effects on Voltammetric Responses at a Disk Microelectrode Covered with a Self-Assembled Monolayer in the Presence of Electron Hopping[END_REF], or nonlinear chemical kinetics [see, e.g., ref . 21].

To examine this point while remaining within the scope of this work, a Butler-Volmer electron transfer kinetics was considered to apply at the disk electrodes surfaces and the electron transfer rate constant k 0 was decreased while keeping the three electrodes potential poised at a constant value E disk = E 0 + 0.3 V. The corresponding results presented in Table 1 evidence that upon introducing an increasingly severe kinetic constraint on the surface concentrations through decreasing k 0 , the three electrodes perform with practically identical current intensities at long times, passing approximately 1/3 of the overall current flowing through each electrode. This is of course only possible if, in the bulk solution, the surface areas of the cross-sections of the diffusional volumes associated with each electrode are also one third of that of the elementary cell surface area.

Conclusion

The results presented in this work establish that Voronoi approaches produce reasonably good predictions of the current intensities flowing through the elementary cells of random band and disk electrode arrays performing under chronoamperometric conditions for a simple electron transfer mechanism. The relative error introduced by Voronoi-based approaches is less than 5% in each case.

However, this work demonstrated that the current intensities featuring the electrochemical activity of individual sites obtained based on Voronoi approaches are overly wrong as soon as the diffusion layers generated by each electrode activity start to overlap. Fortunately, the ensuing errors nearly compensate each other due to redistributions of the flux lines in the solution far from the array surface. Nonetheless, this discrepancy suggests that Voronoi approaches may afford erroneous results when applied to complex kinetics or for predicting the outcome of electroanalytical methods that rely on a fine coupling between surface concentrations and local fluxes.

Acknowledgements

This work was supported in parts by CNRS, ENS, Sorbonne University and PSL University (UMR 8640 PASTEUR). GP acknowledges CNRS Institute of Chemistry (INC) for the "Emergence" post-doctoral support. AO thanks CNRS INC for the award of the "Emergence@INC2020" grant. CA acknowledges the Sino-French LIA CNRS NanoBioCatEchem for its support, as well as Xiamen University for his Distinguished Visiting Professor position.

Table 1. Dependence of the longtime limits of the relative current intensities experienced by each disk electrode 1,2 and 3 of the disk array in Figure 5b on the magnitude of the electron transfer rate constant for a simple ET oxidation mechanism and a Butler-Volmer kinetic law when the simulations are performed globally in the elementary cell. a

Electrode b

Fractions of the total array current intensity for each electrode

k 0 = 10 2 cm.s -1 k 0 = 10 -3 cm.s -1 k 0 = 10 -5 cm.s -1 1 (17%) 29% 29% 33% 2 (29%) 34% 34% 33%
3 (54%) 37% 37% 34%

(a) E disk = E 0 + 0.3 V;  = 0.5; C bulk = 1 mM; r d = 10 µm; D = 10 -5 cm 2 s -1 ; other geometrical parameters of the disk electrode array were identical to those given in Fig. 5 caption. (b) The percentages indicated between parentheses are the values of the k 0 -independent current intensity fractions evaluated upon using a Voronoi approach (i.e., as in Fig. 5c).

Figure captions:Figure 1 .

 1 Figure captions:

Figure 2 .

 2 Figure 2. Simulated currents (solid curves) for (a) the elementary cell shown in (Fig.1b) or (b) the two Voronoi unit cells shown in (Fig.1c) with P = 140 µm, 𝑥 𝛾 𝑃/2 20 μm, w = 10 µm, L = 1 mm (common lengths of the bands) for 𝐶 1 mM and 𝐷 10 cm /s. Green: total current; orange: right band current; blue: left band; dashed lines represent predicted Cottrellian limits: 1,2single and two bands short times limits, correspondingly; 3 -long time limit for the total current in the elementary cell; 4,5 -long time limits for each electrode according to Voronoi unit cells approximation. In (c) the predicted global currents shown in (a) and (b) are superimposed to stress their small difference (orange: current simulated for the global elementary cell; green: sum of the currents of the three Voronoi unit cells); the inset represents the relative error resulting from the use of a Voronoi approach.

Figure 3 .

 3 Figure 3. Concentration fields enforced by the simultaneous operation of the two band electrodes poised at the same potential in the elementary cell (Fig.1b) under different natural convection intensities imposing different limits of the convection layers at: 𝛿 10 (a), 40 (b) and 60 µm (c). The white curves shown in each panel represent marginal flux lines originating at the facing edges of each band, i.e. the lines across which no flux may occur in the perpendicular to them directions, while the dotted black line perpendicular to the substrate surface at the middle point between these edges (𝑥 20 μm) features the constant Voronoi limit. All other parameters are identical to those reported in the Fig. 2 caption.

Figure 4 .

 4 Figure 4. Relative variation of the equivalent (symbols-solid lines) surface areas associated to each band as a function of the position of the middle point 𝑥 of the insulating gap between the two electrodes at long diffusional times. All other parameters are identical to those reported in the Fig. 2

Figure 5 .

 5 Figure 5. (a) 2D pseudo-random array of disks of identical radius (rd) arranged asymmetrically used in this work and (b) its representative elementary cell containing 3 electrode labelled 1, 2 and 3. In (c) is shown the three Voronoi unit cells 1, 2 and 3 created by splitting the elementary cell (b) according to the Voronoi tessellation concept (the unit cells are labelled according to the labels of their constitutive electrodes in (b)). In all this work the squared elementary cell has a 100 µm side length and the three electrodes have common radius (𝑟 10 μm), being centered at [x = 20 µm, y = 20 µm] (1), [x = 60 µm, y = 20 µm] (2), and [x = 40 µm, y = 50 µm] (2). The areas of the Voronoi unit cells in (c) are: 1.66 10 (1), 2.90 10 (2) and 5.43 10 (3) µm 2 , i.e., occupy respectively 17% (1), 29% (2) and 54% (3)of the total 10 4 µm 2 surface area of the elementary cell cross-section.

Figure 6 .

 6 Figure 6. Red solid curves: currents intensities simulated as a function of time under chronoamperometric conditions for (a) the elementary cell shown in Fig. 5b considered as a whole or (b) obtained through summation of the individual currents intensities simulated for the three Voronoi unit cells 1-3 shown in Fig. 5c. Other solid curves provide the time variations of the current intensities for each of the three electrodes, 1 (blue), 2 (yellow) and 3 (green), obtained in each simulation approach. Dashed lines represent predicted Cottrellian limits: 1,2 -single and three disks short times limits, correspondingly; 3 -long time limit for the total current in the elementary cell; 4,5,6 -long time limits for each electrode according to Voronoi unit cells approximation. In (c) the predicted global currents shown in (a) and (b) are superimposed to stress their small difference (orange: current simulated for the global elementary cell; green: sum of the currents of the three Voronoi unit cells); the inset represents the relative error resulting from the use of a Voronoi approach. The chronoamperometric current intensities were evaluated considering that each electrode potential was poised on the plateau of the redox wave of an electroactive species at a bulk concentration 𝐶 1 mM and a diffusion coefficient of 𝐷 10 cm /s. See Figure 5 caption for the geometry of the array.

Figure 7 .FiguresFigure 1 .

 71 Figure 7. Flux lines issued from electrodes 1 and 3 of the disk-array elementary cell (compare to Fig.5b,c) at t = 1 s. For the sake of clarity of the perspective the flux lines issued from electrode 2 are not represented. The location of the 3 disk electrodes on the elementary cell surface are indicated on the

  Figure 1. (a) Cross-sections of the pseudo-random array of parallel bands of identical widths (w) arranged asymmetrically used is this work and (b) of its representative elementary cell whose width is equal to the halfperiod (𝑃/2) of the array in (a). In (c) is shown the cross-section of the two independent unit cells created by splitting the elementary cell displayed in (b) according to the Voronoi tessellation concept. In (b,c) 𝑥 𝛾𝑃/2is the abscissa of the middle point between the two facing edges of the band electrodes.

Figure 2 .

 2 Figure 2. Simulated currents (solid curves) for (a) the elementary cell shown in (Fig.1b) or (b) the two Voronoi unit cells shown in (Fig.1c) with P = 140 µm, 𝑥 𝛾 𝑃/2 20 μm, w = 10 µm, L = 1 mm (common lengths of the bands) for 𝐶 1 mM and 𝐷 10 cm /s. Green: total current; orange: right band current; blue: left band; dashed lines represent predicted Cottrellian limits: 1,2 -single and two bands short times limits, correspondingly; 3 -long time limit for the total current in the elementary cell; 4,5 -long time limits for each electrode according to Voronoi unit cells approximation. In (c) the predicted global currents shown in (a) and (b) are superimposed to stress their small difference (orange: current simulated for the global elementary cell; green: sum of the currents of the three Voronoi unit cells); the inset represents the relative error resulting from the use of a Voronoi approach.

Figure 3 .

 3 Figure 3. Concentration fields enforced by the simultaneous operation of the two band electrodes poised at the same potential in the elementary cell (Fig.1b) under different natural convection intensities imposing different limits of the convection layers at: 𝛿 10 (a), 40 (b) and 60 µm (c). The white curves shown in each panel represent marginal flux lines originating at the facing edges of each band, i.e. the lines across which no flux may occur in the perpendicular to them directions, while the dotted black line perpendicular to the substrate surface at the middle point between these edges (𝑥 20 μm) features the constant Voronoi limit. All other parameters are identical to those reported in the Fig. 2 caption.

Figure 4 .

 4 Figure 4. Relative variation of the equivalent (symbols-solid lines) surface areas associated to each band as a function of the position of the middle point 𝑥 of the insulating gap between the two electrodes at long diffusional times. All other parameters are identical to those reported in the Fig. 2 caption. The corresponding variations for the Voronoi unit cells are represented by dashed lines. 𝑥 35 μm corresponds to symmetrical electrode configuration within elementary cell and thus to equal relative contributions (50%).

Figure 5 .

 5 Figure 5. (a) 2D pseudo-random array of disks of identical radius (r d ) arranged asymmetrically used in this work and (b) its representative elementary cell containing 3 electrode labelled 1, 2 and 3. In (c) is shown the three Voronoi unit cells 1, 2 and 3 created by splitting the elementary cell (b) according to the Voronoi tessellation concept (the unit cells are labelled according to the labels of their constitutive electrodes in (b)). In all this work the squared elementary cell has a 100 µm side length and the three electrodes have common radius (𝑟 10 μm), being centered at [x = 20 µm, y = 20 µm] (1), [x = 60 µm, y = 20 µm] (2), and [x = 40 µm, y = 50 µm] (2). The areas of the Voronoi unit cells in (c) are: 1.66 10 (1), 2.90 10 (2) and 5.43 10 (3) µm 2 , i.e., occupy respectively 17% (1), 29% (2) and 54% (3) of the total 10 4 µm 2 surface area of the elementary cell cross-section.

Figure 6 .

 6 Figure 6. Red solid curves: currents intensities simulated as a function of time under chronoamperometric conditions for (a) the elementary cell shown in Fig. 5b considered as a whole or (b) obtained through summation of the individual currents intensities simulated for the three Voronoi unit cells 1-3 shown in Fig. 5c. Other solid curves provide the time variations of the current intensities for each of the three electrodes, 1 (blue), 2 (yellow) and 3 (green), obtained in each simulation approach. Dashed lines represent predicted Cottrellian limits: 1,2 -single and three disks short times limits, correspondingly; 3 -long time limit for the total current in the elementary cell; 4,5,6 -long time limits for each electrode according to Voronoi unit cells approximation. In (c) the predicted global currents shown in (a) and (b) are superimposed to stress their small difference (orange: current simulated for the global elementary cell; green: sum of the currents of the three Voronoi unit cells); the inset represents the relative error resulting from the use of a Voronoi approach. The chronoamperometric current intensities were evaluated considering that each electrode potential was poised on the plateau of the redox wave of an electroactive species at a bulk concentration 𝐶 1 mM and a diffusion coefficient of 𝐷 10 cm /s. See Figure 5 caption for the geometry of the array.

Figure 7 .

 7 Figure 7. Flux lines issued from electrodes 1 and 3 of the disk-array elementary cell (compare to Fig.

  5b,c) at t = 1 s. For the sake of clarity of the perspective the flux lines issued from electrode 2 are not represented. The location of the 3 disk electrodes on the elementary cell surface are indicated on the bottom face of the prismatic diffusion volume. Note that the scales are identical for all three directions. The dashed lines in green represented on the top and bottom faces of the prism represent the boundaries separating the three Voronoi unit cells in Fig. 5c.

Supplementary Material

Modelling Diffusion at Random Arrays of Electrodes:

Revisiting the Voronoi Tessellation Concept

Giovanni Pireddu, a Irina Svir, a,b Christian Amatore, a,c Alexander Oleinick a, * a PASTEUR, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 24 rue Lhomond, 75005 Paris, France, b Design Automation Department, Kharkiv National University of Radio Electronics, Nauky Avenue, 14, Kharkiv 61166, Ukraine, c State Key Laboratory of Physical Chemistry of Solid Surfaces, College of Chemistry and Chemical Engineering, Xiamen University, Xiamen 361005, China.

___________________

___________________

A. Description of physico-mathematical equations

We consider here a simple electron transfer (E) mechanism: 𝐴 → 𝐵 𝑒 under amperometric conditions with sufficiently high overpotentials at the working electrode to ensure diffusion limited conditions (except for the simulations whose results are presented in Table 1; see Table 1 caption). It is assumed that only species A is initially present initially in the solution at a concentration C bulk , diffusion coefficients of species A and B are supposed to be equal and supporting electrolyte is present in excess in order to ignore migration and ohmic drop contributions. In this case, the mass transport problem only requires the solution of the second Fick law: 𝜕 𝐴 𝜕𝑡 𝐷Δ 𝐴 associated to the boundary conditions given below; 𝐴 is the concentration of species A, 𝑡 is the time, 𝐷 is the diffusion coefficient, and Δ is the Laplacian operator.

In the main text two kinds of systems are considered so the Laplacian operator is expressed as:

Δ

, for the 2D systems (band-type active sites), or:

Δ

, for the 3D systems (disc-type active sites).

Initial and boundary conditions:

𝑡 0: within all solution 𝐴 𝐶 𝑡 0: 𝑧 0, on electroactive site surface: # 𝐴 0 𝑧 0, on insulator 0 ∀𝑧, symmetry planes

In simulations the condition prevailing at "𝑧 → ∞" was replaced by imposing that 𝐴 𝐶 at the vertical distance 𝑧 6 𝐷𝑡 and the discretized diffusion equation was solved within the volumes (disc array) or the planes (band array) featuring the diffusional volume or area accordingly.

In addition, at any time and any 𝑧 value, a zero-flux condition (𝜕 𝐴 𝜕𝑛 ⃗ ⁄ 0) was applied at each point of the symmetry planes defining the elementary cells (global approach) or at each plane featuring the boundaries of their unit cells (Voronoi approach).

# For Nernstian ET. For non-Nernstian ETs, as in Table 1 of the main text, this condition was replaced by a classical Butler-Volmer condition relating the A surface concentration to its local flux 𝜕 𝐴 /𝜕𝑧.

B. Conformal mapping

Most of the simulations presented in the text were performed with Comsol Multiphysics 5.5 [START_REF]COMSOL Multiphysics Reference Manual, version 5[END_REF] to solve the system of partial derivative equation and boundary conditions described in Section A.

However, in the case of the band array, analytical derivations based on conformal mapping were performed in order to check the accuracy of Comsol results and to identify more accurately the trajectories of flux lines (see Figure 3). This was performed by defining the conformal mapping transformation T relating the complex coordinate 𝜚 𝑥 𝑖𝑧 describing the real space (𝑖 √ 1) to that, 𝜔 𝜉 𝑖𝜂, which defines the conformal space (See Figure SM1):

The other parameters defined as:

where P/2 is the half-pseudo period of the array (Figure 1a), i.e., the total width of the elementary cell (compare Figure 1b in the text and Figure SM1 below). The transformation [T] allows treating a wide variety of experimental situations, in particular those in which the widths of the two band electrodes of the elementary cell may have extremely different widths, it must be emphasized that it is restricted to the cases when 𝜔 𝜔 1 as seen in Figure SM1. In graphical terms this means that it encompasses all the situations leading to diffusional domains with identical widths at infinite distances from the basal array plane (compare Figure SM2).

Indeed, the flux lines separating the elementary cell in the real space are then the image of the vertical straight semi-infinite segment originating at the 𝜔 0; 1/ 𝛾 1 and going to the infinity along the 𝜂-axis as evidenced in Figure SM2. The other main interest of the transformation [T] is that, by construction, it allows a direct tracking of the flux lines in the solution, and to identify immediately those evidencing how the diffusional domains associated to each electrode are reshaping while they advance (Figure SM2a). This property was useful for checking the validity of the results obtained through the more heavy simulations performed with Comsol as proven in Figure SM2b.


V Breier

Brudanin

R Arnold

C Augier

A S Barabash

A Basharina-Freshville

S Blondel

S Blot

M Bongrand

D Boursette

R Breier

V Brudanin

J Busto

A J Caffrey

S Calvez

M Cascella

C Cerna

J P Cesar

A Chapon

E Chauveau

A Chopra

L Dawson

D Duchesneau

D Durand

V Egorov

G Eurin

J J Evans

L Fajt

D Filosofov

R Flack

X Garrido

C Girard-Carillo

H Gómez

B Guillon

P Guzowski

R Hodák

A Huber

P Hubert

C Hugon

S Jullian

A Klimenko

O Kochetov

S I Konovalov

V Kovalenko

D Lalanne

K Lang

Y Lemière

T Le Noblet

Z Liptak

X R Liu

P Loaiza

G Lutter

M Macko

C Macolino

F Mamedov

C Marquet

F Mauger

A Minotti

B Morgan

J Mott

I Nemchenok

M Nomachi

F Nova

F Nowacki

H Ohsumi

G Oliviéro

R B Pahlka

V Palusova

C Patrick

F Perrot
email: fperrot@cenbg.in2p3.fr

A Pin

F Piquemal

P Povinec

P Přidal

Y A Ramachers

A Remoto

J L Reyss

B Richards

C L Riddle

E Rukhadze

R Saakyan

R Salazar

X Sarazin

J Sedgbeer

Yu Shitov

L Simard

F Šimkovic

A Smetana

K Smolek

A Smolnikov

S Söldner-Rembold

B Soulé
email: soule@cenbg.in2p3.fr

I Štekl

J Suhonen

C S Sutton

G Szklarz

H Tedjditi

J Thomas

V Timkin

S Torre

Vl I Tretyak

V I Umatov

I Vanushin

C Vilela

V Vorobel Aa

D Waters

F Xie

A Žukauskas Aa

The Nemo-3 Collaboration

Search for the double-beta decay of 82 Se to the excited states of 82 Kr with NEMO-3

Keywords: Double beta decay, Neutrino, 82 Se, Excited State

Search for the double-beta decay of 82 Se to the excited states of 82 Kr with NEMO-3

Introduction

The search for the neutrinoless double-beta decay (0νββ) is of major importance in neutrino and particle physics. Its observation would prove the Majorana nature of the neutrino and would be the first evidence for lepton number violation. Up to now, no evidence of such a process has been found and the best half-life limits are in the 10 24 -10 26 y range [START_REF] Arnold | NEMO-3 Collaboration) (NEMO-3 Collaboration)[END_REF][2][3][4]. 82 Se is one of the best isotopes to investigate 0νββ decay. In particular, its high Q ββ -value of 2997.9±0.3 keV [5] lies above the main backgrounds coming from natural radioactivity. There exist also well-known methods of Se isotopic enrichment through centrifugal separation. This is why 82 Se is the baseline isotope for past, current or future experiments such as LUCIFER [6], CUPID-0 [7] and SuperNEMO [8]. Several studies have been performed in the past to search for 0νββ decay of 82 Se to the ground state of 82 Kr and recently new limits on the half-life have been obtained with the NEMO-3 (2.5×10 23 y [9]) and CUPID-0 experiments (3.5×10 24 y [10]).

The double-beta decay with emission of two neutrinos (2νββ) is a second order electroweak process in the Standard Model. It allows the experimental determination of the Nuclear Matrix Elements (NME) for such processes and provides a robust test for the different nuclear models. It could constrain the quenching factor of the axial-vector coupling constant g A and give the possibility to improve the quality of NME calculations for 0νββ decay [11][12][13][14]. This process has been observed for 11 double-beta isotopes with a range of measured half-lives between 10 18 -10 24 y [15,16]. For 82 Se, several experiments have measured the 2νββ decay to the ground state with the most precise half-life value to date of 9.39 ± 0.17(stat) ± 0.58(syst) × 10 19 y measured with the NEMO-3 experiment [9].

The search for ββ decay to excited states is also an interesting way to study such processes. Indeed, these decays have a very clear-cut signature using the 2e1γ channel (to the 2 + 1 state) or the 2e2γ channel (to the 0 + 1 or 2 + 2 state) which can dramatically reduce the number of background events. The disadvantages are a lower Q ββ available energy for electrons which suppresses the probability of the decay and a lower detection efficiency for electrons and γ-rays. Nevertheless, the decay to excited states is of importance to test the nuclear matrix elements. A detailed analysis for 2νββ decay of 100 Mo and 150 Nd to the excited 0 + 1 state of 100 Ru and 150 Sm, respectively, showed that corresponding NME are only suppressed by ∼ 30% when compared with the NME to ground state transition [17][18][19][20][21][22][23].

Up to now, the 2νββ decay to excited states has only been observed for two isotopes: 100 Mo and 150 Nd with typical half-lives of 10 20 -10 21 y [START_REF] Barabash | AIP Conf. Proc[END_REF]. It is important to note that this decay has been observed only to the 0 + 1 excited state (with the emission of two γ-rays) which is favoured compared to the decay to the 2 + 1 or 2 + 2 excited states. These measurements have been performed using both High Purity Germanium (HPGe) detectors by measuring only the γ-rays in the cascade [18,21,22,[START_REF] Barabash | [END_REF][26][27][28][29] and "trackercalorimeter" detectors such as NEMO-3 able to measure the energies of both electrons and γ-rays [19,[START_REF] Blondel | Optimisation du blindage contre les neutrons pour le démonstrateur de SuperNEMO et analyse de la double désintégration β du néodyme-150 vers les états excités du samarium-150 avec le détecteur NEMO 3[END_REF]. For 82 Se, there is up to now no evidence for such a decay. Stringent limits have been obtained by the LUCIFER collaboration for the (2ν+0ν)ββ decay to various excited states of 82 Kr using a HPGe detector [START_REF] Beeman | [END_REF]. Nevertheless, this technique using only γ-rays does not allow to distinguish between 2νββ and 0νββ. More recently, more stringent limits have been set by the CUPID-0 collaboration for the 0νββ decay to various excited states of 82 Kr using ZnSe scintillating bolometers [32].

In this work, we will present a detailed study of the 82 Se 2νββ and 0νββ decays to the 0 + 1 excited state of 82 Kr, expected to be the most favoured [33,34], with the full exposure of the NEMO-3 experiment. In this analysis, we have access to the full topology of the decay. It consists of the emission of two electrons sharing 1510.2 keV of energy and accompanied by two γ-rays with energies of 711.2 keV and 776.5 keV respectively, as illustrated in Figure 1. After a presentation of the NEMO-3 detector, the 82 Se source foils and the associated backgrounds, we will present a dedicated analysis tool called gamma tracking (GT) developed to reconstruct efficiently the γ-rays in such a decay. Finally, we will present the results of the 2νββ and 0νββ decays of 82 Se to the 0 + 1 excited state of 82 Kr with the full NEMO-3 exposure of 4.42 kg•y. 82 Se source foils and associated backgrounds 2.1. NEMO-3 detector NEMO-3 was a detector installed in the Modane Underground Laboratory (LSM) under 4800 m water-equivalent in order to be protected against cosmic muons. It took data from February 2003 to January 2011. It consisted of a hollow cylinder divided into 20 sectors hosting thin source foils from 7 different enriched isotopes with a typical thickness of approximately 50 mg/cm 2 (as shown in Figure 2). The main isotope to search for 0νββ decay was 100 Mo with a total mass of 6.914 kg. The second isotope of interest was 82 Se with a mass of 0.932 kg shared in 3 sectors. The five other isotopes studied were by decreasing order of mass: 130 Te (0.454 kg), 116 Cd (0.405 kg), 150 Nd (36.55 g), 96 Zr (9.43 g), and 48 Ca (6.99 g) (see [START_REF] Arnold | NEMO-3 Collaboration) (NEMO-3 Collaboration)[END_REF]35] for more details).

NEMO-3 detector,

The source foils were hung at the center of a wire chamber composed of 6180 cells operating in Geiger mode. The gas was a mixture composed of 94.85% helium, 4% ethanol, 1% argon and 0.15% of water vapour. These cells were placed inside a 25 G magnetic field produced by a solenoid surrounding the detector. Charged particles thus had a curved trajectory when crossing the tracking chamber, which allowed the identification of a negative curvature for 95% of electrons at 1 MeV. The minimal distance traveled by a particle crossing the tracker is ∼1.1 m, which corresponds to a typical minimal time of flight of 3 ns. The resolution of the tracker was 0.5 mm transverse to the wires and 8 mm in the vertical direction for 1 MeV electrons.

A calorimeter enclosed the wire chamber. It was made from 1940 plastic scintillator blocks, each one with a typical size of 20 x 20 x 10 cm 3 and coupled to a low background photomultiplier (PMT) through a light guide. The calorimeter measured the kinetic energy of the particles and the time difference between two distant hits could be recorded. The blocks had an energy resolution of 6 -7%/ E(MeV) and a time resolution of 250 ps (σ at 1 MeV).

NEMO-3 was a unique detector as it combined tracking and calorimetry techniques. A charged particle (e -, e + ...) was identified when going across and ionizing the wire chamber gas. Its track was associated to an energy deposit in a calorimeter block neighbouring the last fired Geiger cell. γ-rays were identified when energy was deposited in a calorimeter block but no track was associated. Alpha particles were identified as straight, short tracks as they could not travel more than ∼ 40 cm in the tracker due to their high ionisation energy loss.

In order to run in low background conditions, the NEMO-3 detector had to be protected from natural radioactivity. To do so, a passive shielding of 19 cm iron was surrounding the detector in order to stop external γ-rays. In addition, borated water, paraffin and wood were also used to moderate and absorb the environmental neutrons. For a more detailed description of the NEMO-3 detector, see [35].

82 Se source

Two different batches of 82 Se source were used (referred to as 82 Se(I) and 82 Se(II)). Those batches had an enrichment factor of 97.02 ± 0.05% and 96.82 ± 0.05% respectively. To produce source foils, the enriched 82 Se powder was mixed with polyvinyl alcohol (PVA) glue and deposited between ∼23 µm thick Mylar foils, producing composite source foils. The total mass of the 82 Se isotope in NEMO-3 was 932.4 ± 5.0 g. An analysis of these 82 Se foils was conducted in order to search for 2νββ and 0νββ decays to ground state and is detailed in [9].

Backgrounds

With its powerful topology reconstruction ability, the NEMO-3 detector was able to identify 2e2γ events that were selected for ββ decay to excited states. However, some background isotopes could also produce this type of event. Among them, 214 Bi and 208 Tl decays were the main sources of background as the produced particles could carry similar energies as the β and γ particles from double beta decays to excited states. These two isotopes are β -emitters from the 238 U and 232 Th radioactive decay chains, respectively, with Q-values of 3.27 and 4.99 MeV.

The main background contribution came from contamination in the source foils introduced during isotope production and residual contamination after isotope purification or during the foil production. This is described as internal contamination. In this case, those β emitting isotopes could produce two electrons coming from the same vertex via β-decay with internal conversion, β-decay followed by Møller scattering or β-decay to an excited state with a Compton scattering of the emitted photon. From these mechanisms, additional γ-rays could be produced by bremsstrahlung or from a decay to an excited state as presented in Figure 3. Prior to their installation, the activity of 82 Se foils in 214 Bi and 208 Tl had been measured by low background gamma spectrometry using HPGe detectors. These small contaminations had been also measured and cross-checked by the NEMO-3 detector itself thanks to its capability to measure own background. In NEMO-3, the 214 Bi contamination of the foils could be studied by looking for the so called BiPo effect using 214 Bi and 214 Po sequential decay events. The β-decay of 214 Bi is followed by the α-decay of 214 Po with a half-life of 164.3 µs. The analysis channel used to study such events was the 1e1α(n)γ channel. 208 Tl decays exclusively to excited states and emits mostly 2 or 3 γ-rays (99.9%). Its contamination was thus measured through the 1e2γ channel with a high γ-rays efficiency (about 50% at 1 MeV). Results and comparison of the 214 Bi and 208 Tl activities for 82 Se source foils using HPGe and NEMO-3 data are presented in Table 1 (including some other minor background isotopes [9]).

Finally, 2νββ decay to the ground state was also considered as a background source for 2νββ decay to excited states. When two electrons were produced, two extra γ-rays could be emitted via bremsstrahlung. A 2e2γ event was thus detected while excited states were not involved. Table 1: Results of the contamination measured in the 82 Se source foils by using independently NEMO-3 and HPGe data. All uncertainties are of statistical origin only, given at the 1σ level. The limit shown is at the 2σ level. The activities of 214 Bi and 208 Tl are derived from this independent analysis and are consistent with the ones already published in [9].

In addition to the internal contamination of the source foils, radioactivity from other components of the detector can produce background events, leading to γ-rays. These γ-rays then interact with the source foil and two electrons coming from the same vertex can then be reconstructed if there is either pair production with misreconstruction of the positron track, double Compton scattering or simple Compton scattering followed by Møller scattering of the produced electron. In the case of pair production, there can be annihilation of the positron which produces two photons. Considering that γ-ray interactions are involved in all those mechanisms, they have to be taken into account in the search of the 2νββ and 0νββ decay to the 0 + 1 excited state, with 2 γ-rays emitted in cascade. The different processes responsible for background production are described in Figure 4. The radioactivity of these external elements was first screened by low background γ-spectrometry. Also, when background isotopes produce a γ-ray, it can interact close to the surface of a calorimeter block and produce an electron. The latter crosses the whole wire chamber including the source foil. The initial γ-ray can also deposit energy in the calorimeter before interacting with the source and producing an electron. The contamination of external elements can thus be measured through two channels : crossing electron or (γ, e) external, i.e. Compton scattering in a scintillator block, producing a γ-ray energy deposit, followed by a Compton scattering in the source foil, emitting an electron detected in another scintillator block. An external background model was produced and can be found in [38]. A specific external background is the radon background. It comes from 222 Rn, a gaseous isotope in the 238 U chain. 222 Rn can be introduced via several mechanisms including emanation from detector materials, diffusion from laboratory air through detector seals or contamination of the wire chamber gas. This is only possible because of its long half-life of 3.82 days. Once inside the detector, mainly positive ions are produced from the radon decays. Because of their charge, they can drift and be deposited on the source foils or tracker wires. There, they decay into 214 Bi near the source material. This contamination can then be observed through the 1e1α(n)γ channel.

For the first 18 months of data-taking, there was a relatively high level of 222 Rn inside the detector. To reduce it, an anti-radon tent was built around the detector reducing the radon level inside the wire chamber volume by a factor ∼ 6 [START_REF] Arnold | NEMO-3 Collaboration) (NEMO-3 Collaboration)[END_REF]. The higher radon activity data-taking period is referred to as Phase 1 and the lower activity period that came after as Phase 2.

Both data and Monte Carlo simulations (MC) of signal and background are processed by the same reconstruction algorithm. The DECAY0 event generator [39] is used for generation of initial kinematics and particles are tracked through a detailed GEANT3 based detector simulation [40].

Gamma tracking technique

In most double-beta-decay experiments, a crucial aspect is to precisely measure the energy of the particles. Using the unique combination of tracking and calorimetry, NEMO-3 extracts other observables (angle between two electrons, track curvature, vertex position...) allowing a good discrimination of background and signal events. In addition, one of the most important features is the measurement of the time of flight of the particles inside the detector.

When looking for double-beta decays, selecting events with two electrons from the same vertex is not a strong enough criterion as seen in section 2.3. The time of flight measurement thus allows to reject external events by testing two hypothesis : the event has an internal or an external origin. This test is made possible in NEMO-3 by the knowledge of the particle track length, energy, time of flight and the energy and time resolution (σ t ) of the calorimeter. It can be conducted for charged particles for which tracks are reconstructed but also for γ-rays coming from the same vertex.

Time of flight for electrons is thus an essential parameter when looking for doublebeta decay. The next section will describe its measurement in NEMO-3 before a new method for measuring γ-ray time of flight is presented. The latter is crucial since a more accurate description of events containing γ-rays and a higher sensitivity to these events will improve the efficiency and precision in the search for decays to excited states.

Time of flight calculation

In order to construct an hypothesis on the time ordering of an event, some energy must be deposited in at least two calorimeter blocks and one particle track or more must be reconstructed inside the wire chamber. This track also has to be associated to one of the calorimeter hits. The other calorimeter hit with no associated track is identified as a γ-ray. Figure 5 illustrates an event sketch in NEMO-3 with an electron (one reconstructed track with one calorimeter hit) and a γ-ray (only a calorimeter hit) coming from the same vertex.

Before making any time of flight calculation for an event, two hypotheses must be considered : internal or external origin. Theoretical times of flight between the vertex and the calorimeter hit (t th ) that should be measured by the calorimeter (for each block hit) are then calculated for both hypotheses. The sum (external origin) or difference (internal origin) of these theoretical times is compared to the difference between times actually measured by the calorimeter (t exp ). If the given hypothesis is favoured, then the difference noted t hyp must be close to zero, taking into account the time resolution of NEMO-3. Considering the example presented in Figure 5(b), these differences for both hypotheses are expressed by the following equations :

t int = (t th e -t th γ ) -(t exp e -t exp γ ) (1)
t ext = (t th e + t th γ ) -(t exp e -t exp γ ) (2)
Nevertheless, the calculation of t hyp is only a preliminary analysis. A more advanced study is based on the probability of time of flight and needs to take into account the uncertainties on theoretical and measured times. Thus the χ 2 method is used as described by :

χ 2 hyp = 2 t hyp σ 2 tot , (3)
where σ 2 tot is the quadratic sum of all uncertainties affecting time measurements or calculations. These are the uncertainties on track lengths (for charged particles), path lenghts (for γ-ray), measured energies (due to calorimeter energy resolution) and times (due to calorimeter time resolution).

When considering external or internal events, as in section 4, the selections will be based on the chi-squared probabilities for the respective hypotheses.

Gamma tracking

Another type of time of flight calculation is possible considering only the trajectory of photons. Because of the thickness of NEMO-3 scintillators, γ-rays do not always deposit all their energy inside a single block. One photon can deposit part of its energy in a calorimeter block after Compton scattering, then hit another one and potentially more. Gamma tracking is an original and powerful analysis tool developed recently [START_REF] Hugon | Analyse des données de l'expérience NEMO3 pour la recherche de la désintégration double bêta sans émission de neutrinos. Étude des biais systématiques du calorimètre et développement d'outils d'analyse[END_REF] in order to take this effect into account and reconstruct the complete trajectory of γ-rays inside the detector, with each step from one scintillator to the next.

When a single γ-ray is produced inside a source foil with one or more charged particles and hits several scintillators, a few PMTs measure energies without associating them to reconstructed tracks. Figure 6(a) describes the approach presented in the previous section, where every unassociated hit is considered as having a different origin. Here, the second unassociated block is neither internal nor external and the event is excluded when selecting events for the 2e1γ channel. Using gamma tracking, the same event can be properly reconstructed as shown in Figure 6 In this example, the complete reconstruction of the photon can be done with only one time of flight probability calculation : between the two scintillators not associated to any track. However, when events include several unassociated calorimeter blocks, every combination has to be taken into account and evaluated. In that case, before making a complete calculation, the probability of time of flight is determined for each pair of blocks in the event, with once again the χ 2 method. All the pairs are then combined to extract all possible topologies, each associated to a combined time of flight probability, using the equation :

χ 2 GTtot = n=m n=1 χ 2 GT (Block n-1 Block n ), (4)
with m the total number of blocks involved in the chain of hit scintillator by a single γ-ray, Block n-1 Block n a pair of calorimeter blocks and χ 2 GT (Block n-1 Block n ) the χ 2 value calculated for each pair. The main drawback of this method is the computation time. To limit this effect, two additional conditions are applied : requiring an energy threshold of 150 keV for the energy deposit in each calorimeter block and only taking into account the probabilities greater than 0.1% for any combination.

Once all calculations have been performed, the topology with the highest probability is considered the most likely. This combination defines the number of photons in the event and their trajectories. The gamma tracking technique is thus key to the study of double-beta decays to excited states.

Validation of gamma tracking using calibration sources

During the data taking phase of NEMO-3, several calibration runs using three pointlike 232 U radioactive sources, labelled 1, 2 and 3, were conducted. Their activities were measured through γ-spectrometry (HPGe detectors) and are given in Table 2, column 2. These sources are especially well suited for gamma tracking studies since they decay to the 228 Th nucleus which belongs to the natural 232 Th radioactive decay chain. At the end of the chain, it produces a 208 Tl nucleus which is a β -emitter producing at least two γ-rays : e.g. 2.615 and 0.583 MeV.

We measured the activities of the 232 U sources using NEMO-3 analysis with and without gamma tracking. We can then compare the results with the activities measured by HPGe detectors. The main objectives are to confirm that the use of the gamma tracking method improves the signal efficiency and reduces the systematics. Several criteria are defined to only select events involving one electron and two γ-rays (1e2γ) since 99.8% of 208 Tl decays produce these three particles. Using Monte-Carlo simulations, the efficiency with gamma tracking is determined to be 1.16% for this topology (compared to 0.92% without gamma tracking) while 53082 data events are selected for source 3 with an acquisition time of 107.6 hours. The same analysis without gamma tracking, conducted on the same data sample, selected only 38956 events. About 27% more events were thus selected using gamma tracking, proving that part of the events involving Compton scattering are recovered, thus improving the efficiency. Figure 7 illustrates the number of scintillator blocks hit by a single γ-ray according to the path reconstruction calculated with the gamma tracking method. A reasonable agreement is obtained between data from the 232 U radioactive sources and Monte-Carlo simulations.

Futhermore, the 232 U sources activities obtained with gamma tracking are presented in Table 2 where they are compared to activities obtained without the use of gamma tracking and to the γ-spectrometry measurements. The interaction of γ-rays may induce low energy deposits in the bulk of the scintillator block. The energy response of the calorimeter does not take into account the interaction point in the scintillator block so the effect of the energy threshold may be difficult to simulate. However, the main observation is that activities measured using gamma tracking are more consistent with the γ-spectrometry results. However, activities measured through the analysis with gamma tracking are consistently lower than γ-spectrometry values. This difference is used as a way to estimate the systematic uncertainty induced by the use of the gamma tracking technique. The difference for sources 1, 2 and 3 are respectively 10.4%, 6.5% and 2.0% as reported in Table 2. As a conservative approach, the systematic uncertainty is considered to be 10%. In order to select 2e2γ events, several criteria are applied to distinguish them from background events. The candidate events must contain two electron tracks, originating from the 82 Se source foil, each with an energy deposit greater than 150 keV. The distance between the tracks' intersections with the foil should fulfill ∆ XY less than 4 cm (perpendicular to the wires) and ∆ Z less than 8 cm (parallel to the wires) so they can be considered to have a common vertex. Two γ-rays must be reconstructed using the gamma tracking technique, each with a total energy greater than 150 keV. The timing of the calorimeter hits for electrons and γ-rays must be consistent with an internal event defined as those particles simultaneously emitted from their common vertex in the 82 Se foil. There should be no α-particle tracks and no extra reconstructed γ-rays in the event. 77 data events were selected from a total of 897,409,450 in the selenium sectors for the selected runs. Figure 8 shows that this number is compatible with the number of background events expected when using these criteria, as well as the energy distribution of both electrons for data events and background. Using these criteria, the efficiency for the expected signal is 0.078%. These preselection criteria can be applied when looking for events including two internal electrons and γ-rays. In order to be more specific to the 2νββ(0 + gs → 0 + 1 ) decay, 13 an optimisation is made considering the energies of the four particles for this decay. The first energies to be optimized are the individual energies of both electrons labelled E e min and E e max . Both energies for each event are displayed using two-dimensional histograms for signal and total background, obtained from MC simulations as shown in Figures 9(a) and 9(b). For each bin of this two-dimensional histogram, the local statistical significance (noted N l σ ) is calculated and displayed in Figure 9(c). This value is defined by the following equation :

N l σ = S l √ S l + B l , ( 5
)
where S l is the signal and B l the background in each bin. The signal is given by the 2νββ to 0 + 1 state simulation with a half-life of 3 × 10 20 years which is three times higher than the 2νββ decay to the ground state half-life.

The result of the optimization procedure was tested for several Monte Carlo samples, including 2νββ to the 0 + 1 state with various half-lives. If the half-life in the sample is different from 3 × 10 20 y, the selection would not be optimal, thus the sensitivity to the 2νββ to the 0 + 1 excited state would be decreased. For samples with half-lives larger than 3 × 10 20 y, the optimization procedure gives a too loose selection w.r.t optimum, increasing the background contribution. For samples with half-lives smaller than 3 × 10 20 y, the optimization procedure results in a too strict selection, reducing the signal efficiency. Even if not optimal, the selection would not bias the half-life of the sample.

A selection criterion is defined on N l σ for the maximised total statistical significance N σ as presented in Figure 9(d). N σ is calculated over the total number of simulated signal events and expected background.

In Figure 9(a), simulations show that the signal is stronger when the energies of the two electrons E e max and E e min are in the range of [300-400] and [200-300] keV, respectively. This is due to their primary kinetic energies (with a total energy shared equal or lower than 1510.2 keV) slightly affected by the loss of energy in the source foil and in the tracking chamber. Concerning the background, the simulations in Figure 9(b) show that the energies of the two detected electrons can be much higher, up to 2.7 MeV, than those for the signal. This is due to the presence of 208 Tl isotope (Q-value of 4.99 MeV) which is one of the main backgrounds. Nevertheless, the optimization is able to remove all the events with high energy electrons, typically greater than 1.1-1.2 MeV as illustrated in Figure 9(c).

Other selections are then made on the total electron energy and total γ-rays energy and finally on the two γ-rays' individual energies as seen respectively in Figures 10 and11. Figure 10 shows some of the features of Figure 9, whereby the total energy of γrays for background can be greater than 1500 keV due to higher energy γ-rays emitted in 208 Tl decays (usually 2.61 and 0.58 MeV). The signal simulation fits the 2νββ(0 + gs → 0 + 1 ) transition with the two electrons sharing 1512.2 keV and two γ-rays with a total energy of 1487.7 keV. The optimisation process then only selects events with γ-rays sharing less than 1600 keV, taking into account the energy resolution of the detector. Figure 11 represents the third step of optimization and concerns individual γ-rays energies. By this stage, most of the 208 Tl induced events have been removed. Simulations indicate that most of the remaining background events contain two γ-rays of [300-400] and keV. These can be related to 214 Bi, since its decay can produce a 609.3 keV γ-ray and a lower energy one through bremsstrahlung, shown in Figure 3(a). Finally, most signal events are expected to have two γ-rays of [400-500] and [500-600] keV, corresponding to the 2νββ(0 + gs → 0 + 1 ) γ-rays of 711.2 and 776.5 keV. After the complete optimization process described here, the selection efficiency for the 2νββ(0 + gs → 0 + 1 ) signal calculated from MC is 0.069% with a total of 19 selected data events. The total electron energy distributions for Phase 1 and Phase 2 can be seen in Figure 12 while the total γ-rays energy distributions are shown in Figure 13. These figures also show the different background contributions that are detailed in Table 3. The largest contribution (52% in Phase 2) comes from internal contamination of the source foils and especially from 214 Bi. Radon is also responsible for 68% of background events during Phase 1 and was reduced to 28% in Phase 2. The external backgrounds account for 21% of the total expected background despite the strong criteria used to ensure that only internal events are selected.

It is also shown that there is a good compatibility with background and data events. 3: Numbers of expected background events from the main background sources in both Phases and their contribution to the total number of expected background events for the 2e2γ channel after optimisation for the study of 2νββ(0 + gs → 0 + 1 ) decay. 0.93 year of data taking are considered for Phase 1 and 3.82 years for Phase 2. The quoted uncertainties represent the statistical and systematic uncertainties, respectively. The number of selected data events for each phase is also presented.

In the absence of a significant excess of data versus background, a limit has been set. This can be performed using the following equation :

T 1/2 > × N nuc × ln(2) × (t acq -t d ) × 1 N ex , (6)
where is the detection efficiency, N nuc the number of 82 Se nuclei, t acq and t d the acquisition and dead time respectively and N ex the number of signal events that can be excluded. The method used here to obtain this last number is the CLs method [START_REF] Junk | [END_REF], that takes into account the shape of the expected signal and backgrounds as well as the number of data events and several statistical and systematical uncertainties. The systematics are detailed in Table 4. Considering then the 4.42 kg.y exposure, the 0.069% efficiency, the 21.4 expected background events and 19 data events, the limit on the 2νββ (0 + gs → 0 + 1 ) decay half-life for 82 Se is, at 90% CL :

T 2ν 1/2 ( 82 Se, 0 + gs → 0 + 1 ) > 1.3 × 10 21 y. ( 7
)
This result is compatible with limit of 3 × 10 21 y from Ref. [37] and lower than the value published by the LUCIFER collaboration, who determined a limit of 3.4 × 10 22 y [START_REF] Beeman | [END_REF] for the (2ν+0ν)ββ processes. However, the NEMO-3 technique precisely identifies the event topology and could thus independently study 2νββ and 0νββ processes. electrons do not share energy with neutrinos contrary to the 2νββ decay. The signal efficiency using these criteria increases by a factor 10 compared to the 2νββ process. It reaches 0.71%, as higher energy electrons are expected. The selection has then been optimized with these energies, taking into account a simulated signal with a half-life of 3 × 10 21 years and using the same method as the one described in Section 4.1. Applying those criteria, the final selection efficiency for this signal is 0.69% and 14 data events are selected. The total electron energy distributions for Phase 1 and Phase 2 are shown in Figure 14. The background composition is similar to what was presented in Table 3 with a high radon contribution. The details are presented in Table 5. The total γ-rays energy distributions are shown in Figure 15.

Neutrinoless double beta decay to

As for the 2νββ(0 + gs → 0 + 1 ) transition, data is consistent with background-only predictions so a limit has to be set on the half-life of the 0νββ(0 + gs → 0 + 1 ) process. The method used to calculate such a limit remains the CLs method. With the 20.1 background events and the statistical and systematic uncertainties, the limit on the 0νββ (0 + gs → 0 + 1 ) decay half-life for 82 Se (at 90% CL) is :

T 0ν 1/2 ( 82 Se, 0 + gs → 0 + 1 ) > 2.3 × 10 22 y. ( 8
)
This result is given for the 0νββ(0 + gs → 0 + 1 ) transition for 82 Se, separately from 2νββ(0 + gs → 0 + 1 ). It is compatible with limit of 3.4×10 22 y and 8.1×10 22 y obtained in the LUCIFER experiment [START_REF] Beeman | [END_REF] for the (2ν+0ν)ββ processes and CUPID-0 [32] experiment. According to the mass mechanism, a Majorana neutrino is exchanged during such a process and therefore a limit can also be set on the effective mass of the neutrino using 5: Numbers of expected background events from the main background sources in both Phases and their contribution to the total number of expected background events for the 2e2γ channel after optimisation for the study of 0νββ(0 + 1 → 0 + 2 ) decay. 0.93 years of data taking are considered for Phase 1 and 3.82 years for Phase 2. The quoted uncertainties represent the statistical and systematic uncertainties, respectively. The number of selected data events for each phase is also presented.

the following equation :

1

T 0ν 1/2 M M = G 0ν (Q ββ , Z)g 4 A M 0ν 2 m ββ m e 2 , (9)
where G 0ν (Q ββ , Z) is the phase space factor given in [43] for the transition, g A = 1.27 and M 0ν the nuclear matrix element [12,13,44,45]. The limit that can be set on the effective neutrino mass is m ββ < [42 -239] eV.

Summary and Conclusions

4 Using an innovative gamma tracking technique, the NEMO-3 data set was analysed to search for ββ decays of 82 Se to the excited states of 82 Kr with a 4.42 kg.y exposure. No evidence for the 2νββ process was found and thus an upper limit on the decay half-life was set at 90% CL : T 2ν 1/2 ( 82 Se, 0 + gs → 0 + 1 ) > 1.3 × 10 21 y. This result can nevertheless help to constrain theoretical QRPA models presented in [33,34,46].

The analysis of the 0νββ decay to excited states was conducted in a similar fashion and, as once again no extra events were observed over the expected background, an upper limit was set at 90% CL : T 0ν 1/2 ( 82 Se, 0 + gs → 0 + 1 ) > 2.3 × 10 22 y. These results are obtained for the first time with a detector which reconstructs each particle individually in the final state.

This analysis performed with 82 Se in NEMO-3 will also provide useful information for the next-generation SuperNEMO experiment which will host 100 kg of 82 Se, such as optimisation of the selected events and identification of the main background contributions.

In parallel with its search for 0νββ decay to the ground state, SuperNEMO will also look for the 2νββ and 0νββ decays to excited states with major improvements. Using thicker scintillators, the sensitivity to γ-rays and efficiency to 2νββ and 0νββ(0 + gs → 0 + 1 ) transitions will be enhanced. Backgrounds will also be reduced : more than a factor 30 for radon and a factor 100 for 214 Bi and 208 Tl. The expected sensitivities for SuperNEMO are respectively ∼ 10 23 y and ∼ 10 24 y for the 2νββ and 0νββ(0 + gs → 0 + 1 ) half-lives. A first module, called Demonstrator, with 7 kg of 82 Se is undergoing commissioning and will start taking data in 2019. Its goal is to reach a sensitivity on the 0νββ halflife of 5 × 10 24 y in 17.5 kg . y exposure with the demonstration of a "zero"-background experiment [47].

Figure 1 :

 1 Figure 1: Decay scheme of the 82 Se ββ decay to the 0 + 1 excited state with the emission of two electrons sharing 1510.2 keV and two prompt γ-rays with energies of 711.2 and 776.5 keV [36].

Figure 2 :

 2 Figure 2: Cross-sectional view of the NEMO-3 detector. The detector consists of source foils (1), scintillators (2), photomultipliers (3) and a wire chamber (4).

Figure 3 :

 3 Figure 3: Mechanisms producing 2e2γ events from internal contamination of β emitters inside the source foils. β-decay to excited state followed by Møller scattering and bremsstrahlung (3(a)), β-decay to excited state with internal conversion and double bremsstrahlung (3(b)), β-decay to excited state followed by Compton scattering and bremsstrahlung (3(c)).

Figure 4 :

 4 Figure 4: Mechanisms producing 2e2γ events from external contamination of the NEMO-3 detector emitting a γ-ray interacting inside the source foil. Double Compton scattering of the external γ-ray with bremsstrahlung in 4(a), pair production from the external γ-ray with double bremsstrahlung effects and poor reconstruction of the positron in 4(b), Compton scattering of the external γ-ray followed by a Møller scattering of the electron and a bremsstrahlung in 4(c).

Figure 5 :

 5 Figure 5: Event sketch with track reconstruction (5(a)) and scintillator association (5(b)). The reconstruction defines a γ-ray as energy deposit in a scintillator without any associated track. It can link it to the vertex in Figure 5(b).

  (b) : the second unassociated hit is paired with the first one under the assumption of Compton scattering and the event satisfies the 2e1γ channel conditions.

Figure 6 :

 6 Figure 6: Example of an event reconstruction without using gamma tracking (6(a)). Only one of the two scintillators not associated to a track is consistent with the internal hypothesis, the other is neither internal nor external. The same event is reconstructed with gamma tracking (6(b)), the second scintillator can be associated to the first one under the assumption of Compton scattering.

Figure 7 :

 7 Figure 7: Number of scintillator blocks hit by a single γ-ray according to the path reconstruction calculated with the gamma tracking method, logarithmic scale. Data were acquired using γ-rays from the 232 U radioactive sources and are compared to Monte-Carlo simulations.

4 .

 4 Double beta decay to the excited states 4.1. Two neutrino double-beta decay to 0 + 1 excited state As mentioned in Section 1, ββ decays to the 0 + 1 excited state consist in the simultaneous emission (compared to the NEMO-3 time resolution) of two β and two γ particles.

Figure 8 :

 8 Figure 8: Sum of the electron energies distributions in the 2e2γ channel after the preselection criteria described in the text are applied, for Phase 1 in Figure 8(a) and Phase 2 in Figure 8(b). Data are compared to the MC prediction for the different backgrounds. The background coming from the 2νββ 82 Se decay to g.s. is completely negligible and thus not visible in the two plots.

Figure 9 :

 9 Figure 9: The distributions of signal 2νββ(0 + gs → 0 + 1 ) and background events from MC simulation are represented in Figures 9(a) and 9(b) respectively, as a function of both electrons' individual energy. The local statistical significance N l σ distribution for the 2νββ(0 + gs → 0 + 1 ) transition as a function of both electrons' individual energy is calculated for each bin of the 2-D histogram and represented in Figure 9(c). The total statistical significance Nσ as a function of a cut on the local significance in Figure 9(d) allows the optimization of this cut (dotted red line). Selected bins with high local statistical significance in Figure 9(c) are separated from the removed ones by the dotted red line.

Figure 10 :

 10 Figure 10: Total electron energy vs total γ-rays energy distributions for 2νββ(0 + gs → 0 + 1 ) signal simulation (10(a)) and background (10(b)). Local statistical significance distributions for each bin of this histogram (10(c)) with optimisation cut (dotted red line) on total statistical significance Nσ presented in Figure 10(d).

Figure 11 :

 11 Figure 11: γ-ray 1 energy vs γ-ray 2 energy distributions for 2νββ (0 + gs → 0 + 1 ) signal simulation (11(a)) and background (11(b)). Local statistical significance distributions for each bin of this histogram (11(c)) with optimisation cut (dotted red line) on total statistical significance Nσ presented in Figure 11(d).

Figure 12 : 2 Figure 13 :

 12213 Figure 12: Total electron energy distributions after selection for the 2νββ(0 + gs → 0 + 1 ) transition, for Phase 1 in Figure 12(a) and Phase 2 in Figure 12(b). Experimental data events are compared to the MC simulation for the different backgrounds. The dotted red line represents the simulated signal for a half-life of 3 × 10 20 years.

Figure 14 :

 14 Figure 14: Total electron energy distributions after selection for the 0νββ (0 + gs → 0 + 1 ) transition, for Phase 1 in Figure 14(a) and Phase 2 in Figure 14(b). Experimental data events are compared to the MC simulation for the different backgrounds. The dotted red line represents the simulated signal with a half-life of 3 × 10 21 years.

Figure 15 :

 15 Figure 15: Total γ-rays energy distributions after selection for the 0νββ (0 + gs → 0 + 1 ) transition, for Phase 1 in Figure 15(a) and Phase 2 in Figure 15(b). Experimental data events are compared to the MC simulation for the different backgrounds. The dotted red line represents the simulated signal with a half-life of 3 × 10 21 years.

Table 2 :

 2 Comparison of the 232 U sources activities measured by respectively γ-spectrometry (HPGe detector) and NEMO-3 analysis without and with the gamma tracking technique using the 1e2γ topology. Uncertainties in column 2 are respectively statistics and systematics. Columns 4 and 6 present the relative differences between HPGe and analysis activities (without and with gamma tracking).

	232 U	HPGe (Bq)	No GT act. (Bq) ∆ noGT (%) GT act. (Bq) ∆ GT (%)
	1	7.79 ± 0.04 ± 0.21	6.56 ± 0.08	15.8	6.98 ± 0.07	10.4
	2	15.91 ± 0.09 ± 0.43	13.92 ± 0.13	12.5	14.88 ± 0.11	6.5
	3	32.76 ± 0.17 ± 0.89	30.00 ± 0.17	8.4	32.11 ± 0.14	2.0

Table

			Expected		Contribution to
			events		total background (%)
			Phase 1	Phase 2	Phase 1	Phase 2
		214 Bi	1.14 ± 0.05 ± 0.12 4.28 ± 0.09 ± 0.43	12.1	36.1
	Internal	208 Tl Others 0.06 ± 0.03 ± 0.01 0.29 ± 0.14 ± 0.02 0.43 ± 0.02 ± 0.07 1.58 ± 0.04 ± 0.23	4.6 0.6	13.3 2.4
		Total	1.64 ± 0.07 ± 0.20 6.15 ± 0.49 ± 0.68	17.3	51.8
	Radon	6.44 ± 0.63 ± 0.65 3.26 ± 0.31 ± 0.33	68.0	27.5
		214 Bi	0.38 ± 0.19 ± 0.04 1.49 ± 0.75 ± 0.15	4.0	12.5
	External	208 Tl Others 0.74 ± 0.37 ± 0.08 0.78 ± 0.39 ± 0.08 0.29 ± 0.10 ± 0.03 0.18 ± 0.06 ± 0.02	2.9 7.8	1.6 6.6
		Total	1.39 ± 0.43 ± 0.15 2.46 ± 0.05 ± 0.85	14.7	20.7
	Total background 9.47 ± 0.77 ± 1.00 11.87 ± 1.17 ± 1.26	100.0	100.0
	Data events	7	12	-	-

Table 4 :

 4 Values of the 1σ systematic uncertainties included in the calculation of the limits on 2νββ decay to excited states and their methods of estimate. The estimated uncertainties come from the comparison of the activity measurements of calibration sources between NEMO-3 and HPGe ( 232 U, 207 Bi, 90 Y), the uncertainties on background measurements and uncertainties specific to the detector or 82 Se sources.

Table

			Expected		Contribution to
			events		total background (%)
			Phase 1	Phase 2	Phase 1	Phase 2
		214 Bi	1.25 ± 0.05 ± 0.13 4.57 ± 0.10 ± 0.46	13.5	41.9
	82 Se foils	208 Tl Others 0.02 ± 0.01 ± 0.01 0.07 ± 0.03 ± 0.01 0.24 ± 0.01 ± 0.03 0.82 ± 0.02 ± 0.12	2.6 0.2	7.6 0.6
		Total	1.50 ± 0.06 ± 0.17 5.46 ± 0.11 ± 0.59	16.3	50.1
	Radon	6.50 ± 0.64 ± 0.65 2.97 ± 0.27 ± 0.30	70.4	27.3
		214 Bi	0.45 ± 0.22 ± 0.05 1.33 ± 0.62 ± 0.14	4.9	12.3
	Detector	208 Tl Others	0.79 ± 0.09 ± 0.08 1.13 ± 0.14 ± 0.12 0 0	8.5 0	10.3 0
		Total	1.24 ± 0.24 ± 0.13 2.46 ± 0.62 ± 0.26	13.4	22.6
	Total background 9.24 ± 0.69 ± 0.95 10.89 ± 0.69 ± 1.15	100.0	100.0
	Data events	6	9	-	-

Acknowledgements We thank the staff of the Modane Underground Laboratory for their technical assistance in running the experiment. We acknowledge support by the grants agencies of the Czech Republic (grant number EF16 013/0001733), CNRS/IN2P3 in France, RFBR in Russia (Project No.19-52-16002 NCNIL a), APVV in Slovakia (Project No. 15-0576), STFC in the UK and NSF in the USA.