Maryna Hrytsaienko

Mathieu Gallart

Marc Ziegler

Olivier Cregut

S Tamariz

R Butte

N Grandjean

Bernd Honerlage

Pierre Gilliot

O Crégut

R Butté

B Hönerlage

Dark-level

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Beyond fundamental aspects, the ongoing interest in semiconductor quantum dots (QDs) is motivated by their potential for a wide range of applications such as sources of single and entangled photons [1][2][3], quantum cryptography [4] or spintronics [5]. For a review see for instance [6]. Key conditions to use QDs for nanophotonic applications rely on the full understanding and, if possible, the control of the different dynamics associated with electronic excitations, namely population, spin-and phase relaxation. For example, generating indistinguishable photons from single QDs requires to have a short excitonic coherence time T 2 and a short radiative lifetime T 1 , which is favored in the Purcell regime [7,8].

Among the wide range of available QD systems, GaN/Al(Ga)N QDs offer attractive features because of the large band offsets at the AlN/GaN interface, enabling three-dimensional confinement of electron-hole (eh) pairs up to room temperature [9]. Advantage was taken of this to achieve single UV-photon emission at 300 K [10,11] and even above [12]. Long-lived optical ori- * mathieu.gallart@ipcms.unistra.fr entation was also experimentally observed [13]. Such experimental accomplishments go hand in hand with strong efforts regarding QD growth in order to control their structural parameters, especially their size and surface density [14,15]. A fundamental and intrinsic feature of wurtzite GaN/AlN quantum nanostructures remains the presence of a strong built-in electric field along the [0001] growth direction. This electric field is caused by the difference in both spontaneous and piezoelectric polarizations between the GaN and AlN materials [16], which can reach values of several MV/cm [17,18]. The presence of such a static electric field gives rise to the quantum confined Stark effect (QCSE), which has a strong influence on both the recombination energy and the lifetime of the e-h pairs. By inducing a spatial separation of electron and hole wave-functions, i.e., the formation of a dipole, the electric field decreases the oscillator strength, hence leading to a strong dependence of the radiative lifetime on the dot height [17,19,20] that spans several orders of magnitude from the picosecond to the microsecond range. The consequences are dramatic for the largest QDs. First, as their radiative recombination probability is decreased, other dephasing mechanisms become dominant and the radiative lifetime of those QDs will exceed the most optimistic values of the dephasing time that can be expected for this type of system. Furthermore, their weak oscillator strength will hinder any resonant excitation and the control of coherent phenomena. Finally, the possible transient screening of the electric field is responsible for a complex recombination dynamics, which is non exponential and strongly depends on the number of injected e-h pairs [17]. These are major drawbacks if one considers for instance the possibility of generating and manipulating entangled states. This is what motivates the present work on small GaN/AlN QDs whose emission energy lies above the strained bulk GaN bandgap and for which a significant value of the oscillator strength is preserved. Because of their short radiative lifetime -in the nanosecond range-these QDs are also expected to be less sensitive to any possible screening of the built-in electric field, be it due to the carrier injection regime (low to high) or to surrounding charges.

In this work, we investigate the distribution of radiative lifetime that occurs within a single GaN/AlN QD layer because of the dot size dispersion. We evidence that, beyond the variations stemming from the changes in dot height that give rise to an increasing weight of the QCSE, the fluctuations in their lateral size [21] explain the reduced variation in the radiative lifetime, within a factor about two, for dots whose emission energy covers a spectral range of ∼600 meV.

II. SAMPLE DESCRIPTION AND EXPERIMENTAL SET-UP

The sample consists of a single plane of GaN/AlN QDs deposited on a c-plane sapphire substrate by ammoniasource molecular beam epitaxy (MBE) in the Stranski-Krastanov growth mode. The template is composed of a commercial 1 µm thick AlN buffer layer (DOWA Electronics) grown by metal-organic vapor phase epitaxy and a 100 nm thick AlN layer deposited by MBE. The GaN QDs were grown under constant low ammonia flux (NH 3 beam equivalent pressure of 3 • 10 -7 Torr) at 750 • C and 4.5 monolayer (ML) GaN coverage. After forming the QDs the Ga flux was stopped and the dots were capped with a 55 nm thick AlN barrier. The dot density under such growth conditions is in the 10 10 -10 11 cm -2 range [15]. Time-resolved photoluminescence (TRPL) measurements were performed at low temperature (T = 10 K). The excitation consists in the fourth harmonic of a pulsed ytterbium laser (Tangerine model from Amplitude) emitting at 257 nm (4.82 eV) with a pulse duration of 270 fs and a repetition rate of 50 kHz. The excitation energy is below the AlN bandgap. Moreover, the wetting layer shows a thickness between 2 and 3 MLs. As observed in a similar sample, the steady-state photoluminescence (PL) of the wetting layer occurs at energies distributed between 4.5 and 5.0 eV. We therefore assume that e-h pairs are mostly created in the wetting layer even if the excited discrete states of the QDs likely play a role in the photogeneration process.

The low repetition rate also ensured that the dots relax to their ground state between two successive excitation pulses. The energy per pulse is 0.2 nJ and the diameter of the excitation spot focused on the sample was about 0.4 mm, giving an incident photon flux of 2 • 10 11 photons • cm -2 . The TRPL signal was sent into a 25 cm-focal length spectrometer with a 150 grooves/mm grating and dispersed onto the photocathode of a Hamamatsu streak camera with a spectral resolution of 0.3 nm and a time resolution close to 10 ps.

III. RESULTS

Figures 1(a) and 1(b) present the streak camera image of the QD emission recorded at low temperature (T = 10 K) and the corresponding time-integrated PL spectrum, respectively. The PL spectrum consists of a broad band modulated by fringes due to Fabry-Perot interference occurring within the epilayer stack (see Appendix A). Because of the small size of the GaN QDs, the quantum confinement increases the recombination energy far above the GaN bandgap. A comparison with a sample where the emission of GaN QDs with different heights was investigated [22] indicates that the thickness of those QDs would range between 6 and 8 MLs. Figure 2 presents normalized PL spectra recorded at different delay times after excitation. The emission, which is centered around 3.8 eV and has a full width at half maximum (FWHM) of ∼500 meV, exhibits a continuous redshift as a function of time. Both the quantum confinement and the QCSE depend on the QD height [19,23]. Due to the size-dispersion of the dots, PL emission experiences an inhomogeneous broadening that we assumed to be Gaussian. This is justified by the quasi-linear dependence of the PL emission energy with the QD height as we will show in our calculations (cf. Sec. III C). An energy shift is observed between t = 0 and t = 3.8 ns (∼90 meV) that is much smaller than the FWHM of the PL band (∼500 meV). Two phenomena can contribute to this shift. First, because of the QCSE, the spectral contributions from dots with different heights will display different dynamics. Thus, at long time delay, the contribution from the low energy side of the PL band is enhanced with respect to the high energy emission that decays faster. The larger QDs, which emit at lower photon energies, show indeed a higher QCSE that increases the e-h spatial separation and, consequently, their radiative lifetime. Another explanation could be a possible screening of the built-in electric field: The initial e-h pair population modifies the QD band profile, inducing a blueshift of the PL emission at short time delay that will disappear with the recombination of this population [START_REF] Henini | Handbook of self assembled semiconductor nanostructures for novel devices in photonics and electronics[END_REF].

Unfortunately, both effects have the same physical ori-gin and are always present and interdependent, leading to a complex PL dynamics that is challenging to interpret. However, in the limit of low carrier injection, i.e., for less than one e-h pair per dot, the screening of the electric field becomes negligible. For these reasons, we performed our experiment at low excitation intensity. Taking into account the above-mentioned photon flux and the absorption of the sample ( 3.5 %, measured at room temperature on a commercial UV-Vis-NIR Perkin-Elmer Lambda 950 spectrophotometer), we estimate that, in the limit where each absorbed photon creates one e-h pair, their maximum surface density would be 7 • 10 9 cm -2 . Compared to the surface density of QDs given above of 10 10 -10 11 cm -2 , it is reasonable to assume that there is no more than one e-h pair per dot knowing that e-h pairs will most likely not all end up in a QD (capture efficiency <100%). Hence, screening effects are expected to be negligible in the present set of experiments. This is also supported by modeling of the emission as demonstrated hereafter.

B. Time evolution of the PL spectra

In order to study the PL dynamics as a function of QD height, we have extracted several series of time-profiles integrated on spectral slices within the emission spectrum. If, as discussed above, there is no screening effect, the PL transients should exhibit an exponential decay for sufficiently narrow spectral slices. Obviously, this is not the case as can be seen in Fig. 3 that displays a few typical decay profiles spectrally integrated over 5 nm wide regions together with the corresponding fitting curves. The fitting function is the sum of three contributions: two exponential decays, together with a constant background. The latter represents a very long recombination process that is not resolved in our experiments. It will be discussed later on. For the first two contributions, at each spectral position, the fast component does not depend on the photon energy while the slower one follows the energy dependence expected in the context of the QCSE. In the following, we denote τ L and τ S the time constants related to the slow (L: long delays) and fast (S: short delays) components, respectively.

PL dynamics at short and very long delays

The fact that, before reaching a plateau, the dynamics is bi-exponential indicates that we are dealing with two independent relaxation channels. Concerning τ S , its value is close to 0.3 ns. Being independent of the photon energy (and therefore of the dot size), it does not contribute to the spectral shift described in the previous paragraph. Since τ S is not affected by the dot height, it leads us to assume that it is related to a non-radiative process. Unfortunately, our experimental data do not enable us to identify in an unambiguous manner the nature of the electronic states that constitute the ending step of this transfer. However, as this short decay is not usually seen in experiments performed on non-unintentionally doped high quality GaN/AlGaN quantum wells (QWs) (see for exemple [START_REF] Stokker-Cheregi | [END_REF]), we postulate that these features are specific to QDs. Processes involving the trapping of a single carrier outside the dots, like in colloidal nanocrystals (see for example [START_REF] Gómez-Campos | Hole surface trapping in cdse nanocrystals: Dynamics, rate fluctuations, and implications for blinking[END_REF]), are very unlikely, especially when taking into account the large difference in bandgap energy between GaN and AlN. Indeed, it would imply bringing an energy of several hundreds of meV to ionize the QDs, which can be excluded at 10 K. As less than one e-h pair is created within each dot, scattering events between two such pairs giving rise to an ionization are also unlikely. Moreover, such processes where a single carrier remains within the dot are known to favor nonradiative recombination like Auger processes that are not observed here. Instead, we could have either an extrinsic process, that involves for example the trapping of a carrier at a point defect within the QD or at the interface with the barrier material, or an intrinsic process, like the relaxation toward the fundamental state of the QD that would be a dark level with a different spin configuration [11].

We can now focus on the very long component that shows up as a constant background in our experiments. It could be tentatively explained by the reverse process where the e-h pairs that have experienced a relaxation toward a dark state will slowly flip back to their bright radiative state.

PL dynamics at intermediate delays

Unlike the short decay characterized by τ S , τ L shows an evolution with the photon energy that is characteristic of a height-induced variation of the radiative lifetime of a QD that undergoes the QCSE. Indeed, the experiments being performed at 10 K, the decay that we measure is purely radiative. In Fig. 4, we have plotted the corresponding recombination rate γ (E) = (τ L (E))

-1 as a function of emission energy. Its variation can be written as a Taylor expansion to second order around E 0 , the initial position of the PL maximum:

γ(E) = γ 0 + γ 1 (E -E 0 ) + 1 2 γ 2 (E -E 0 ) 2 , (1)
with

γ 1 = dγ dE E0 and γ 2 = d 2 γ dE 2 E0 .
FIG. 4: Evolution of the experimental radiative rate γ(E) deduced from the PL decays recorded at T = 10 K as a function of QD emission energy (black squares) and its fit using a second order polynomial (red continuous line).

The best agreement is obtained for

E 0 = 3.869 eV, (2)
γ 0 = 0.76 ± 0.02 ns -1 , (3)
γ 1 = 1.2 ± 0.1 ns -1 eV -1 , (4)
γ 2 = 3.1 ± 1.3 ns -1 eV -2 . ( 5
)
These values show large uncertainties, mainly because fitting the dynamics at different energies can lead to approximate values of γ in the spectral region where the signal to noise ratio of the PL is poor. We thus used another method that allows a better precision especially on the second derivative: We perform an analysis by fitting the entire spectrum recorded at different delay times (Fig. 2) in order to provide a better averaged value of γ(E).

Unsurprisingly, the spectra displayed in Fig. 2 fit nicely with a Gaussian lineshape after removing the oscillations due to the Fabry-Perot effect (see Appendix A). Indeed, the homogeneous linewidth is negligible with regard to the inhomogeneous broadening that gives rise to the observed linewidth due to the size distribution of the QDs. If screening effects are negligible, the relation between the dot height and the emission energy is expected to be a bijection. The radiative recombination probability is independent of time and the dynamics for one given height of quantum dots is exponential. Then, considering an ensemble of QDs with heights following a normal distribution, the temporal evolution of the corresponding PL intensity spectrum I P L (E, t) can be written:

I P L (E, t) ∝ e - (E-E 0 ) 2 Γ 2 inh • e -γ(E)t , (6)
where Γ inh is the inhomogeneous broadening, E 0 is the energy of the PL maximum at zero delay, and, as above, γ (E) = 1/τ L (E) is an implicit function describing the energy dependence of the radiative lifetime. As long as γ (E) contains only linear and quadratic terms in E, the spectral shape will remain Gaussian at any time, with simply a width and a central value that change together with the maximum intensity.

The first derivative of Eq. ( 6) with respect to energy cancels at the maximum of the Gaussian curve, as detailed in Appendix B:

E max (t) = E 0 - Γ 2 inh γ 1 • t 2 + Γ 2 inh γ 2 • t . ( 7
)
To first order, the time evolution of the PL maximum position is thus a linear function of time. The FWHM is given by the second derivative of Eq. ( 6):

∆E(t) = 2 Γ inh 2 + Γ 2 inh γ 2 • t . (8)
The last Gaussian parameter is the intensity:

I P L (E max , t) ∝ exp -γ 0 • t + 1 2 Γ 2 inh γ 2 1 2 + Γ 2 inh γ 2 t • t 2 . ( 9
)
Using this analytic approach, the temporal behavior of the PL lineshape provides information about a possible screening of the electric field.

The time evolution of the PL lineshape parameters as well as the corresponding fits are plotted in Fig. 5, while the intensity time evolution is shown in Fig. 6. For the latter, we have to take into account the fast component decay τ S = 0.3 ns. Therefore, we performed the fit (red curve in Fig. 6) for times longer than 0.6 ns. The experimental data are nicely described by the above-mentioned expression demonstrating that the temporal evolution of the PL line is mainly due to the variation in the radiative lifetime within the QD distribution. The fitting procedure gives:

Γ inh = 0.242 ± 0.001 eV, ( 10
)
γ 0 = 0.77 ± 0.02 ns -1 , (11)
γ 1 = 1.05 ± 0.02 ns -1 eV -1 , (12)
γ 2 = 3.08 ± 0.36 ns -1 eV -2 . ( 13
)
As expected, this method leads to an improved precision regarding the numerical coefficients considered in the Taylor expansion of γ(E).

C. Dependence of the radiative lifetime on the dot height

The coefficient γ 1 in the Taylor expansion of γ(E) is positive, which reflects the increase in the recombination probability when decreasing the dot height due to the smaller weight of the QCSE. Hence, ∆E(t) is a monotonous decreasing function of time. The longer the delay time, the narrower the distribution of populated QDs because the larger radiative recombination rate of small QDs is such that only the larger dots emitting at lower photon energies will contribute to the PL spectra, thus leading to the observed decrease in the FWHM.

The time constants extracted by means of the fitting procedure are summarized in Fig. 7 together with literature data related to samples containing different heights FIG. 6: Time evolution of the PL intensity at E max (t). FIG. 7: Radiative lifetime τ L of GaN/AlN QDs as a function of emission energy: black diamonds (present work), red squares (Ref. [17]), blue circles and green hexagons (Ref. [20]). Black down triangles correspond to the fast time constant τ S (present work) ascribed to a nonradiative relaxation process. The continuous lines are results of calculations of τ rad based on a 2D envelope-function modeling for different values of the built-in electric field.

of GaN/AlN QDs. Our experimental data complete nicely the set of data from references [17,20]. In order to estimate the emission energy and the radiative lifetime of e-h pairs in GaN/AlN QDs, the envelopefunction equation was solved numerically for both electrons and holes, using a finite difference method. The confinement within the layer plane, i.e., perpendicular to the growth direction is, in a first step, neglected as the built-in electric field has the most important impact on the optical properties of III-nitride based nanostructures along the c-axis. Thus, the QDs are first mod- eled as two-dimensional (2D) QWs subjected to a static electric field along the growth direction, the carriers being free to propagate in the well plane. We took into account the blueshift of the GaN bandgap to 3.71 eV at low temperature (T = 10 K) due to the strain originating from pseudomorphic growth on the AlN layer.

The radiative lifetime τ rad is then given by the expression: 27] where ω is the photon energy, J eh is the overlap integral of electron and hole envelope-functions and A is an adjustable parameter. The resulting function τ rad ( ω) is plotted as a continuous line in Fig. 7 for different values of the electric field F and A = 20 ns•eV -1 . Such a simple calculation provides curves that qualitatively reproduce the experimental variation in the radiative lifetime over six orders of magnitude for values of the electric field ranging between F = 5 MV/cm and F = 9 MV/cm, the best agreement being obtained for a mean built-in electric field value F = 7 MV/cm. The mean value of our experimental data agrees well with compiled data from the literature showing the variation of the radiative lifetime for a large range of dot heights.

A • ( ω) -3 |J eh | -2 [
In order to be more quantitative, let us note that the QD ensemble in our sample exhibits a large inhomogeneous broadening. To illustrate this aspect, we consider Fig. 8, which shows that the variation in the radiative lifetime does not exactly match the curves (continuous lines) that only take into account the variation of the QCSE due to the changes in dot height. In fact, the fluctuations occurring within the dot ensemble are likely to also originate from the lateral confinement with an increasing weight of the latter when the QD lateral size becomes comparable to the exciton Bohr radius. Indeed, the 2D QW model works better for the largest dots for which variations in the oscillator strength are essentially due to changes in the e-h overlap along the growth axis. For small dots, the in-plane localization of the envelopefunctions is no longer negligible. To estimate this effect on the energy and the lifetime of an e-h pair, we performed a simple calculation based on the variational method (see Appendix C): the in-plane confinement potential is supposed to be step-like with a cylindrical symmetry. The in-plane envelope function of the carriers is considered to be a Gaussian centered in the cylindrical dot and for which the spatial extent is the variational parameter. This method is known to give reasonable qualitative results in spite of its simplicity [28]. For a given QD height, a variation in the localization radius has little influence on the lateral overlap while it strongly affects the confinement energy. As a consequence, the variation in the lifetime versus energy becomes less pronounced. This is illustrated in Fig. 8 that shows as dashed lines the effect of the variation in the lateral dot size on the transition energy and the radiative lifetime for various dot heights. We thus see that the slope of the dashed lines nicely matches the relatively flat variation in the radiative lifetime as a function of photon emission energy observed experimentally. Obviously, more realistic calculations would consider the different lateral confinement experienced by electrons and holes. First, these carriers have different effective masses and the wave function of the lighter particle, the electron, should leak into the barrier. Second, if the dot shape looks like truncated pyramids, as confirmed by high-resolution transmission electron microscopy measurements [29], the confinement potential would act differently on the two particles which are pushed away to opposite sides, the electron to the top of the pyramid and the hole in the wetting layer. Hence, both effects will change the lateral overlap of the electron and hole wave functions [START_REF] Butté | Effects of polarization in optoelectronic quantum structures[END_REF]. However, our approach shows that such a detailed description of the QD wave functions is not mandatory. The lateral confinement plays a major role by blueshifting the emission energy without changing the e-h overlap and thus the radiative lifetime.

Finally, we can interpret these data in the following way: from one dot to the other any variation in the mean radiative lifetime is primarily governed by a change in the dot height that involves the QCSE as already demonstrated in many publications (see, e.g., Fig. 7 that displays literature data). However, the present analysis shows that, within a given QD sample, changes in the QD emission energy by several hundreds of meV accompanied with only slight changes in the radiative lifetime can be understood through fluctuations in the lateral dot dimensions within a given QD subset corresponding to a well-defined height. In this respect, we recall that the identification of GaN/AlN QD subsets with a well-defined height within the same sample was already reported for dots exhibiting a variation in the biexciton binding energy that was correlated to the exciton energy and the dot radius within a QD subset with a given height [22]. Let us note that for QDs emitting around 3.5-4.1 eV, as those studied here, the effect of the built-in electric field is less important than in higher dots where the QCSE has more dramatic effects on the emission energy and the oscillator strength and for which only a mean radius has to be considered. Thus, in dots emitting above the bulk GaN bandgap, which are of prime interest for single photon emission [10,11], the weaker weight of the QCSE enables to reveal otherwise hidden effects originating from changes in the e-h lateral confinement, be it what governs variations in the biexciton binding energy [22] or a relative insensitivity of the radiative lifetime while covering a wide spectral range in the QD emission energy (here, up to ∼600 meV). FIG. 9: Correlation between the calculated transition energy and the QD height for different values of the built-in electric field. The straight lines correspond to calculations based on a 2D QW model (no lateral confinement), while the dashed one correspond to calculations based on the variational method for small dots with a lateral radius of 3 nm.

The reported behavior confirms that the main mechanism behind the relaxation of the e-h pair population is the radiative recombination process and that the radiative lifetime follows the expected evolution for QDs subjected to the QCSE. Moreover, our results strengthen the fact that, if care is taken to avoid any screening of the built-in electric field, the relation between the radiative lifetime and the emission energy tends to a universal trend (Fig. 7). This can be well understood within a simple model including a mean value of the built-in electric field that is the same for all GaN/AlN QD samples.

From a practical point of view, the calculations presented in Fig. 9 allow to estimate that the QD height lies between 5 and 8 MLs (1.3 and 2.08 nm), since 1 ML of GaN along the c-axis is equal to ∼ 0.26 nm [START_REF] Vurgaftman | [END_REF].

IV. CONCLUSION

We have performed low temperature TRPL experiments on a single plane of GaN/AlN QDs, whose emission energy mostly lies above the bulk GaN bandgap. The PL data show a three-component dynamics: a first exponential decay with a short time-constant (0.3 ns), independent of photon energy, a second one that follows the expected evolution for the radiative lifetime with increasing QD height, and finally a very long component with a plateau-like behavior.

We attribute the fast dynamics to a nonradiative relaxation process of e-h pairs toward a dark state and the concomitant existence of a very long component, which shows up as a constant background in the PL decays, that corresponds to the reverse process where the dots come back to their original state.

Concerning the dynamics of the radiative recombination process, we have established that, in the absence of any screening of the built-in electric field, the PL spectrum of QDs with a normal height distribution presents a Gaussian lineshape with a characteristic time evolution that is easily described by an analytical model: after a pulsed excitation the position of the line maximum undergoes a redshift that varies linearly with increasing delay time and is accompanied by a decrease in the FWHM. This behavior originates from the combined effect of the quantum confinement and the QCSE, and constitutes a criterion to conclude on the possible occurrence of screening and/or many-body effects in the population dynamics within GaN/AlN QDs. Moreover, our analysis enables an improved precision concerning the determination of the energy-dependence of the QD radiative lifetime when compared to a method that involves a systematic study of the decay-profiles corresponding to different narrow PL-energy windows.

By completing previous results from other groups, our data provide new evidence that the dependence of the effective PL lifetime with emission energy follows a global trend that is an intrinsic feature of the GaN/AlN QD system where the built-in electric field acts as the most dramatic parameter. Numerical simulations, which show a good qualitative agreement over an energy range of 1 eV and six orders of magnitude for the effective PL lifetime when assuming a built-in electric field value of 7 MV/cm, enable to estimate that the dot height in our sample is distributed around a mean value of 1.6 nm (∼ 6 MLs). Remarkably, in small QDs emitting at high photon energies, around 3.5-4.1 eV, and hence for which the impact of QCSE is less pronounced, we can evidence the effects of lateral confinement on the carriers. The latter affects the wave function extent and the emission energy but only slightly the carrier overlap and the corresponding oscillator strength (equivalently the radiative lifetime).

of each spectral component is exponential with a decay time τ (E).

The PL intensity can then be written:

I P L (E, t) ∝ exp - (E -E 0 ) 2 Γ 2 inh exp (-γ (E) t) , (B1)
where we have introduced the radiative relaxation rate:

γ(E) = 1 τ (E) . ( B2
)
When γ only depends linearly and quadratically on E, i.e., when its Taylor expansion is stopped at the second order, the shape of I P L as a function of photon energy E remains a Gaussian at any time. The energy position of the PL maximum E max (t) and the FWHM ∆E (t) versus time are obtained by canceling the first and the second derivatives of Eq.( 6) relative to the photon energy E, respectively. The calculation of these derivatives is straightforward and gives:

∂ ∂E I P L (E, t) = - dγ dE • t - 2(E -E 0 ) Γ 2 inh • I P L (E, t) , (B3)
and

∂ 2 ∂E 2 I P L (E, t) = -t • d 2 γ dE 2 + t 2 • dγ dE 2 + 4t Γ 2 inh dγ dE (E -E 0 ) - 2 Γ 2 inh + 4 Γ 4 inh (E -E 0 ) 2 • I P L (E, t) . (B4) 1. γ is independent of E If γ = γ 0 is independent of E, the first derivative can- cels for: E max = E 0 , (B5)
and the second derivative for:

E ± = E 0 ± Γ inh √ 2 . (B6)
Hence, the FWHM is given by:

∆E = E + -E -= √ 2Γ inh . (B7)
We recover the simple case of a Gaussian lineshape whose spectral position is independent of time and whose intensity decays with a lifetime τ 0 = γ -1 0 .

2. General case: γ is a function of E If now γ is a quadratic function of E that we can expand around E 0 :

γ(E) = γ 0 + γ 1 • (E -E 0 ) + 1 2 γ 2 • (E -E 0 ) 2 . (B8)
The first derivative of I P L cancels for:

E max (t) = E 0 - Γ 2 inh 2 dγ dE • t = E 0 - Γ 2 inh 2 (γ 1 + γ 2 • (E max (t) -E 0 )) • t, E max (t) = E 0 - Γ 2 inh γ 1 • t 2 + Γ 2 inh γ 2 • t . (B9)
The time evolution of the FWHM is obtained by canceling the second derivative, which gives:

E ± (t) = E max (t) ± Γ inh √ 2 1 + Γ 2 inh 2 d 2 γ dE 2 • t = E 0 - Γ 2 inh γ 1 .t ∓ Γ inh 2 + Γ 2 inh γ 2 .t 2 + Γ 2 inh γ 2 • t . (B10) ∆E(t) = E + (t) -E -(t) = 2 Γ inh 2 + Γ 2 inh γ 2 • t . (B11)
The intensity is then given by the PL emission maximum at E max (t):

I P L (E max , t) ∝ exp -γ 0 • t + 1 2 Γ 2 inh γ 2 1 2 + Γ 2 inh γ 2 t • t 2 .

(B12)

Appendix C: Variational calculation of in-plane wave-functions

A variational calculation enables to estimate the influence of an in-plane confinement potential on the energies and radiative lifetimes of QDs. We consider as a first approximation that the potential is separable in the z and (x, y) directions. The built-in electric field is taken into account along the vertical direction and does not play any role on the lateral confinement. The Hamiltonian of a particle with effective mass m * trapped in an abrupt potential V 0 with a cylindrical symmetry is given by:

H = - 2 2m * 1 ρ ∂ ∂ρ ρ ∂ ∂ρ + V 0 Θ (ρ -ρ 0 ) , ( C1
)
where ρ is the radial coordinate, Θ is the Heaviside step function and ρ 0 is the spatial extension of the lateral potential. The confinement energy can be found using the variational method by minimizing: E (ρ 0 , σ) = φ (σ) |H|φ (σ) .

(C2) as a function of the extension parameter σ, where φ is a normalized Gaussian envelope-function:

φ (ρ) = 1 σ √ π exp - ρ σ 2 . ( C3
)
The problem is solved analytically and gives:

σ = ρ 0 √ 2 ln 2m * ρ 2 0 V 0 2 -1/2 , (C4)
and

E (ρ 0 ) = 2 2m * ρ 2 0 1 + ln 2m * ρ 2 0 2 V 0 . ( C5
)
The procedure is realized for electrons and holes with effective masses m * e = 0.2m 0 and m * h = 0.33m 0 [32], respectively, with m 0 the free electron mass, giving the confinement energies E e and E h , as well as the in-plane extensions σ e and σ h of the corresponding envelopefunctions. The overlap integral can also be evaluated analytically:

FIG. 1 :FIG. 2 :

 12 FIG. 1: (a) Streak camera image of GaN/AlN QDs recorded at T = 10 K using a 4.82 eV pulsed laser for an incident photon flux of 2 • 10 11 photons • cm -2 , (b) corresponding time-integrated PL spectrum.

FIG. 3 :

 3 FIG. 3: Low temperature (T = 10 K) PL decays of GaN/AlN QDs depending on the emission energy corresponding to different heights of dots.

FIG. 5 :

 5 FIG. 5: Time evolution of both the center E max (t) and the FWHM of the Gaussian function ∆E(t) modeling the low temperature PL spectra of GaN/AlN QDs recorded at different delay times.

FIG. 8 :

 8 FIG. 8: Radiative lifetime τ L of GaN/AlN QDs as a function of emission energy. The continuous lines are issued from 2D-QW calculations of τ rad for different values of the built-in electric field assuming a continuous variation in the dot height. The dotted lines show the effect of the QD lateral size with radii ranging between 3 and 20 nm. Each curve corresponds to a different QD height expressed in MLs. The size of the black diamonds is proportional to the PL intensity, which reflects the weight of the corresponding dot height within the QD distribution.

+∞ 0 φ

 0 e (ρ) φ h (ρ) 2πρdρ = 2 σ e σ h σ 2 e + σ 2 h . (C6)

ACKNOWLEDGMENTS This work was funded by the French National Research Agency (ANR) through the Programme d'Investissement d'Avenir under contract ANR-17-EURE-0024 within the Investissement d'Avenir program ANR-10-IDEX-0002-02.

Appendix A: Removing the oscillations from PL spectra

As shown in Fig. 10, PL spectra are modulated by fringes due to Fabry-Perot interferences occurring in the epilayer stack. We remove these oscillations by applying a Gaussian fit to time-sliced PL spectra in order to get access to the time evolution of the fitting parameters: the energy position of the PL maximum E max , its width and its intensity. In order to obtain an analytical expression of the spectro-temporal evolution of the PL spectrum emitted by an ensemble of GaN/AlN QDs, we first make the assumption that the height distribution is normal, resulting in a Gaussian lineshape. If the density of injected e-h pairs is small enough, the recombination probability is time-independent and it exists a function τ (E) describing the variation of the radiative lifetime as a function of emission energy. The consequence is that the dynamics


PhD Daniel Romero

PhD Mireia Calvo

PhD Virginie Le Rolle

Nathalie Béhar

Phillipe Mabo

PhD Alfredo Hernández
email: alfredo.hernandez@inserm.fr

MD Nathalie Behar

MD Philippe Mabo

Formal

Philippe Mabo

; Al- Fredo Hernández

Daniel Software

Romero

Philippe Mabo

Philippe Mabo

Philippe Mabo

Al- Fredo Hernández

Philippe Mabo

Multivariate ensemble classification for the prediction of symptoms in patients with Brugada syndrome

Keywords: Brugada syndrome, depolarization disorders, ensemble classifier, heart-rate recovery, sudden cardiac death

Multivariate Ensemble Classification for the Prediction of Symptoms in Patients with Brugada Syndrome

Daniel Romero 1 • Mireia Calvo 1 • Virginie Le Rolle 2 • Nathalie Behar 2 • Phillipe Mabo 2 • Alfredo Hernández 2 1 Introduction
Brugada syndrome (BS) is a genetic pathology associated with a high risk for sudden cardiac death (SCD) in patients with structurally normal hearts. It is diagnosed when a coved type-1 electrocardiogram (ECG) morphology (STsegment elevation ≥ 2mm) is observed in at least one right precordial lead (V1 and/or V2), placed in the 2nd, 3rd or 4th intercostal space. This pattern can be observed either spontaneously or after a provocative drug challenge test using intravenous Na + channel blockers [START_REF] Priori | ESC Guidelines for the management of patients with ventricular arrhythmias and the prevention of sudden cardiac death: The Task Force for the Management of Patients with Ventricular Arrhythmias and the Prevention of Sudden Cardiac Death of the European Society of Cardiology (ESC)[END_REF][START_REF] Priori | Executive summary: HRS/EHRA/APHRS expert consensus statement on the diagnosis and management of patients with inherited primary arrhythmia syndromes[END_REF].

Implantable cardioverter defibrillator (ICD) therapy is recommended (class I) in BS patients who survived an aborted cardiac arrest and/or have suffered spontaneous sustained ventricular tachycardias (VT). It should also be considered (class IIa) in those patients with a spontaneous type-1 ECG pattern and history of syncope. The decision of implanting an ICD is far more complex in those individuals who remain asymptomatic after BS diagnosis, which represent around 60% of the BS population. Indeed, although these patients show a lower risk of arrhythmic events, as reported in the FINGER study [START_REF] Probst | Long-term prognosis of patients diagnosed with brugada syndrome: Results from the finger brugada syndrome registry[END_REF], this risk is not insignificant and the consequences of not implanting an ICD in these cases is extreme [START_REF] Sacher | Outcome After Implantation of a Cardioverter-Defibrillator in Patients With Brugada Syndrome A Multicenter Study-Part 2[END_REF]. Therefore, one of the challenges regarding the treatment of this pathology is the identification of useful quantitative markers to identify asymptomatic patients that may have a higher risk of SCD and that may properly benefit from the implantation of an ICD.

Several clinical and experimental electrophysiological studies have suggested that both depolarization and repolarization disorders may significantly contribute to the underlying pathophysiology of BS [START_REF] Meregalli | Pathophysiological mechanisms of Brugada syndrome: Depolarization disorder, repolarization disorder[END_REF]. The presence of QRS late potentials (LP), fragmented QRS complexes [START_REF] Morita | Fragmented QRS as a marker of conduction abnormality and a Accepted manuscript predictor of prognosis of Brugada syndrome[END_REF], the prolongation of both the PR and QRS durations [START_REF] Postema | Local Depolarization Abnormalities Are the Dominant Pathophysiologic Mechanism for Type 1 Electrocardiogram in Brugada Syndrome. A Study of Electrocardiograms, Vectorcardiograms, and Body Surface Potential Maps During Ajmaline Provocation[END_REF], as well as the appearance of wider S waves in inferolateral leads and a rightward deviation of the QRS loop axis in its terminal quarter [START_REF] Morita | Sinus node function in patients with Brugada-type ECG[END_REF], are some depolarization disorders that have been associated with BS, mainly related to a slowing conduction within the right ventricular outflow tract (RVOT). Regarding repolarization, most alterations observed from experimental studies have been related to an imbalance of the ionic currents contributing to phase 1 of the action potential (AP) of cardiac cells.

Furthermore, the role of the autonomic nervous system (ANS) seems to have a significant impact in unmasking symptoms in this pathology, particularly during exercise testing. Tan et al. [START_REF] Tan | Exercise-induced ECG changes in brugada syndrome[END_REF] reported that exercise resulted in an increase of the J-point amplitude in BS patients. Makimoto et al. [START_REF] Makimoto | Augmented STsegment elevation during recovery from exercise predicts cardiac events in patients with brugada syndrome[END_REF] evaluated the relationship between the post-stress parasympathetic activation and ST-segment changes. They reported that parasympathetic reactivation during early recovery after exercise, assessed through the heart rate recovery (HRR) index [START_REF] Cole | Heart-rate recovery immediately after exercise as a predictor of mortality[END_REF], tends to be higher in BS patients with prior VF episodes. Furthermore, recent works from our group showed that symptomatic BS patients exhibit greater fluctuations in sinus node response to ANS in 24-h Holter recordings [START_REF] Behar | Heart rate variability and repolarization characteristics in symptomatic and asymptomatic Brugada syndrome[END_REF] and an increased parasympathetic modulation during incremental exercise and early recovery [START_REF] Calvo | Multivariate classification of Brugada syndrome patients based on autonomic response to exercise testing[END_REF].

Accepted manuscript

Multivariate approaches that involve several electrophysiological markers have demonstrated their usefulness for assessing VF risk in BS patients [START_REF] Kawazoe | Risk stratification of ventricular fibrillation in Brugada syndrome using noninvasive scoring methods[END_REF], increasing prognosis accuracy. Since both repolarization and depolarization disorders mentioned above are manifested as a consequence of this genetic disease, studies considering both scenarios simultaneously may provide more meaningful insights than if considered separately. Moreover, the assessment of cardiac electrophysiological markers from both ventricular repolarization and depolarization periods, in addition to other autonomic-related parameters during controlled autonomic maneuvers, such as physical exercise testing, could have a great impact in the prognosis of BS patients.

In this study we aimed at designing a specific machine learning (ML) pipeline to: 1) identify potentially predictive risk markers in BS patients based on their symptomatology that may offer new insight; 2) propose a multivariate classifier to discriminate between symptomatic and asymptomatic BS patients, with a good enough trade-off between performance and generalization. To achieve these goals, depolarization and repolarization ECG-derived features were obtained from 12-lead ECG recordings acquired during a standardized exercise stress test. The HRR marker was also considered as an autonomicrelated parameter. The parameters with the greatest predictive potential were identified through a resampled-based regularization approach to obtain sparse and interpretable models, with high generalization capability. Then, the final model was validated after further model assessment of an ensemble classifier proposed herein, using a wrapper-based strategy for balancing simplicity and performance.

Materials and methods

Study population

The study population comprises 110 consecutive patients (mean age = 44.6 ± 13.7 years) suffering from Brugada syndrome, enrolled in a prospective, multicentric study led by the cardiology department of the Rennes University Hospital (CHU de Rennes), in collaboration with other french hospitals. (1996 and 2000). All patients provided their written informed consent to participate in the study.

From the total population, 25 patients presented clinical symptoms including syncope or aborted SCD related to VF (symptomatic group) based on their medical history. The remaining patients (asymptomatic group) had not presented BS-related symptoms when the study was conducted. As indicated on

Accepted manuscript

active guidelines, all patients on the symptomatic group have received an ICD. On the asymptomatic group, 19 of 85 (22%) have received an ICD, based on a positive EPS (Electrophysiological Study) test. Age of participants ranged from 19 to 74 years old (mean: 44.6 ± 13.7) and 74.5% were males. Structural heart disease was discarded by physical examination, patient's clinical history, and resting and exercise electrocardiogram. Moreover, no significant left ventricular hypertrophy was observed on any patient during echocardiographic screening. More details of this clinical study, including information on the genetic background, can be found on [START_REF] Behar | Heart rate variability and repolarization characteristics in symptomatic and asymptomatic Brugada syndrome[END_REF]. Table 1 summarizes the baseline characteristics of this population.

During a physical exercise test performed on a cycle ergometer (Ergoline 900 Egamed, Piestany, Slovakia), standard 12-lead ECG recordings sampled at 1000 Hz were acquired for all patients. The protocol used for this standardized test is defined as follows:

-Exercise phase (EX): 2 minutes of initial workload by pedaling at 50 W (30 W for women) followed by successive increments of 30 W (20 W for women) every 2 minutes until the patient reached at least the 80% of his/her maximal theoretical heart rate, defined by HR max = 220 -age. -Recovery phase (RE): two successive 3-minute periods, consisting of an active recovery (pedaling at a workload of 50 W), followed by a passive recovery at rest.

Preprocessing

ECG signals were preprocessed before the automatic extraction of the analyzed features. This preprocessing included automatic QRS complex detection and subsequent visual inspection to guarantee normal beat annotations, cubic spline interpolation for baseline drift attenuation, 4 th order Butterworth lowpass filtering at 45 Hz to remove muscular noise and wave delineation using an evolutionary optimization approach [START_REF] Dumont | Improving ECG beats delineation with an evolutionary optimization process[END_REF].

ECG-derived measures

Depolarization parameters

Classical and unconventional QRS parameters were evaluated at relevant stages of the exercise test for each individual lead. The classical group of parameters included the amplitudes of the R and S waves (R a and S a ), the duration of the QRS complex (QRS d ) and the maximum vector magnitude (VM) within the spatial QRS loop. The latter was evaluated from the orthogonal leads X, Y and Z, generated by applying the Dower inverse matrix on the standard leads V1-V6, I and II [START_REF] Edenbrandt | Vectorcardiogram synthesized from a 12-lead ECG: Superiority of the inverse Dower matrix[END_REF].

The unconventional group of parameters included the main slopes (R U , R D , S U ) and angles (φ R and φ S ) identified on a QRS complex. R U and R D Accepted manuscript refer to the up-and down-strokes of the R wave, respectively, while S U is the final up-stroke of the S wave. Details on how to compute these parameters can be found in [START_REF] Romero | Depolarization changes during acute myocardial ischemia by evaluation of QRS slopes: Standard lead and vectorial approach[END_REF]. The angles φ R and φ S can be calculated from the measured QRS slopes, using (1) and [START_REF] Behar | Heart rate variability and repolarization characteristics in symptomatic and asymptomatic Brugada syndrome[END_REF], where constants represent the conversion factors needed to standardize angular measures in order to match conventional ECG printouts used in clinical practice [START_REF] Romero | Detection and quantification of acute myocardial ischemia by morphologic evaluation of QRS changes by an anglebased method[END_REF].

φ R = arctan R U -R D 0.4(6.25 + R U R D ) . ( 1
)
φ S = arctan R D -S U 0.4(6.25 + R D S U )
.

(2)

Repolarization parameters

The spatial dispersion observed in the T wave morphology can be quantified by certain parameters computed from the principal components analysis (PCA) applied on the available ECG leads. The equations below define such parameters:

T-wave residuum:

T Ri = L l=4 δ i,l / L l=1 δ i,l (3)
T-wave uniformity:

T Ui = δ i,1 / L l=1 δ i,l (4)
T-wave complexity:

T Ci = L l=2 δ i,l / L j=1 δ i,l (5)
T-wave non-planarity:

T Ni = δ i,3 /δ i,1 (6)
where i is the i-th beat and δ i contains the eigenvalues of the interlead repolarization correlation matrix, Rei , sorted in descending order. The matrix Rei can be generated as Rei = N -1 n=0 e i (n)e T i (n), with e i (n) being the vector for the multilead ECG signal, corresponding to the T-wave interval of the the i-th beat. Note that for each beat, the N × L matrix E i = [e i,1 , e i,2 , . . . , e i,L ] is obtained, where each column representing the lead l can be seen as an observation of a random process e with zero empirical mean [START_REF] Laguna | Techniques for Ventricular Repolarization Instability Assessment from the ECG[END_REF].

The parameter T Ri represents the relative energy of the non-dipolar components of the ECG model described in [START_REF] Acar | Spatial, temporal and wavefront direction characteristics of 12-lead T-wave morphology[END_REF], with respect to the total energy (sum of all the eigenvalues). In that model, the sum of the first three components are considered as to represent the energy of the dipolar components in normal conditions. If local heterogeneities are present in the repolarization, larger eigenvalues associated with the non-dipolar components are expected to increase T Ri measures.

Accepted manuscript

T Ui and T Ci are designed to characterize the ST-T loop morphology. T U values closer to 1 are indicative of narrower ST-T loops that typically lie in the direction of the first eigenvector δ 1 . Likewise, T C quantifies the roundness of the loop, so that values close to 1 indicate an ST-T loop that is mainly contained in a plane. Finally, T N Pi accounts for the non-planarity of the ventricular repolarization.

Spatial dispersion of the repolarization can also be assessed through the analysis of the following ECG intervals: the corrected QT interval (QT c ) using the Bazett's formula, the T wave width, (T w ), the T-wave peak to T-wave end interval (T pe ) and the ratio between T pe and QT , (T pe /QT ). All these features were calculated in leads V5 and II, usually used for QT related measures in most clinical studies. In order to quantify T-wave symmetry, the following shape-based parameters were also assessed: the ratio of the left-(T LA ) and right-side (T RA ) areas with respect to T-wave peak, T area = T LA T RA , and the ratio of the time intervals lasted between the T-wave peak and the onset (T op ) and end of the T-wave (T pe ), respectively, T time =

Tp-To

Te-Tp = Top Tpe [START_REF] Langley | Quantification of T wave shape changes following exercise[END_REF].

Heart rate recovery

To characterize recovery dynamics, strongly related to vagal response during this period, the heart rate recovery (HRR) curve was assessed by means of the following bi-exponential decay model previously used in [START_REF] Romero | A biexponential approach for assessing parasympathetic reactivation after submaximal exercise[END_REF]:

θ HRR (t) = A exp (-τ 1 t) + B exp (-τ 2 t) (7)
where τ 1 and τ 2 are the decay rates of two exponential processes, associated with the active and passive periods of recovery in the defined protocol. Parameters A and B account for the amount of change in HR decay caused by the two exponential processes controlled by τ 1 and τ 2 .

The HRR curve was fitted after resampling the instantaneous heart rate series (expressed in beats/min) at 4 Hz to have evenly sampled values. Then, the difference between the actual value θ HRR (t) and the initial value θ HRR (t 0 ) was computed for t = {60, 120, 180} seconds, to obtain the HRR values at minutes 1, 2, and 3 after the peak effort (HRR 1 , HRR 2 , and HRR 3 ), respectively. Figure 1 shows an example of the HR evolution for a particular patient and the resulting HRR curve.

Feature extraction

Mean values ( Ȳl s ) of the ECG markers were computed using beat-to-beat measures from 15-s windows in each lead l, and at three relevant stages (s) including exercise onset (T1), the peak effort (T2) and the end of recovery (T3). For T2 and T3 stages, the 15-s windows were taken just before the maximum heart rate and minute 6 th of recovery, respectively. The changes observed between T1 and T2 (exercise period), and between T2 and T3 (recovery period) were

Accepted manuscript

then determined and defined by ∆Y l EX = Ȳl T2 -Ȳl T1 and ∆Y l RE = Ȳl T3 -Ȳl T2 , respectively, where Y represents the actual marker. Finally, the initial feature database, F , was created including the mean values of all parameters at each stage s, Ȳl s , and their dynamics during exercise (∆Y l EX ) and recovery (∆Y l RE ), together with the HRR at the first, second and third minute after the peak effort. Thus, a total of 198 features were initially included in F : 80 concerning cardiac depolarization, 115 from ventricular repolarization and 3 related to the autonomic control.

Pipeline used for designing the proposed classifier

The proposed BS classifier was designed based on two separate phases (see Fig. 2). Phase 1 is applied to obtain a subset of features with high predictive power for symptoms and is based on the elastic net regularization algorithm using repeated k-fold cross-validation. Phase 2 further exploits this reduced subset of features using a wrapper approach, where the final subset (model) is selected by repeated holdout cross-validation applied on non-regularized models. These two phases are explained in more detail in the following sections.

Phase 1: Selection of high predictive features using elastic net

Since F is composed of nearly 200 features, some of them probably highly correlated, the first step is to identify the most predictive set of features. In this work, the dimension of the original feature space is greater than the number of patients available. We have thus chosen to apply the elastic net method for feature selection.

Elastic net is a regularized regression method combining LASSO (L 1 -norm) and ridge regression (L 2 -norm) penalties for the regularization of the model coefficients [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. This approach solves the limitations of the LASSO method when dealing with highly correlated variables, which tends to select one variable and ignore the others. Elastic net makes it possible to obtain sparse models and thus to better interpret final outcomes. Two parameters are involved in the elastic net regularization: α defines the weight of L 1 versus L 2 optimization, while λ defines the regularization strength used in the optimization process. For 0 < α < 1 and a non-negative λ, elastic net solves the problem:

min β0β ( 1 2N N i=1 y i • (β 0 -x T i β) -log(1 + e β0-x T i β ) + λP α (β)), ( 8
)
where N is the number of observations in the dataset, y i is the response at observation x i , a vector containing the input features, and P α (β) represents the penalty term that interpolates between the L 1 and L 2 norms of β coefficients, defined as:

P α (β) = (1 -α) 2 ||β|| 2 2 + α||β|| 1 = p j=1 ( (1 -α) 2 β 2 j + α|β j |) (9)

Accepted manuscript

Parameter α must be defined between 0 and 1. Values of α=1 and α=0 represent the LASSO and ridge regressions, respectively. Other values of α represent elastic net optimization. We have chosen an α value closer to LASSO regression in order to select relevant variables.

Parameter λ can take different values, so that several models with distinct sets of predictive variables can be obtained. An optimal λ value generates a model that seamlessly generalizes new, previously unseen data, while achieving the right balance between simplicity and adjustment of the training data. If λ is too high, the model will be sparse, but with the risk of underfitting the data. Conversely, if λ is too low, the model will be more complex and thus there will be a risk of overfitting the data. Therefore, the best predictive model (with the optimal λ) can be found by selecting the lowest cross-validated error, estimated as the expected deviance of the model applied to new data not seen during training, among all the evaluated λ values.

To find this optimal λ value (hyperparameter optimization), stratified k-fold cross-validation (CV) procedure is often applied [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. We have implemented a repeated (N =100) K-fold cross-validation approach (K = 10) through each λ i value, with i = 1, ..., I (number of λ values), as follows:

1. For n from 1 to N -Split pseudo-randomly the dataset F in K folds -For k from 1 to K -Define L k (without the k-th fold) and T k (the k-th fold) as the training and test subsets.

-For i from 1 to I • Build a regularized, CV model M i cv (L k ) = f (L k , λ i ). • Evaluate M i cv (L k ) on
T k and store prediction errors. -Calculate the CV error for iteration n in each λ i , n (λ i ). 2. For each λ i calculate the average error over the N realizations, ˆ (λ i ). 3. Define λ opt as the optimal value with the minimal average error. 4. Fit the regularized models M i ncv = f (F, λ i ) without CV (ncv) in order to get the variables (with non-zero coefficients) selected for each λ i , including all models from λ opt to λ min (the simplest model).

The implementation of regularized logistic regression using elastic net, as presented above, requires the application of a feature scaling method, because features' coefficients with large variance are less penalized. As a result, the update speed for the coefficient estimates during the optimization approach depends on their feature scales, and when converged, some coefficients are closer to the optimum than the others. Feature scaling helps with this issue, because coefficients are at the same scale and update roughly with the same speed. Therefore, all variables entering in the regularized features selection process were scaled using the following expression:

fq = f q -min(f ) max(f ) -min(f ) (10)
Accepted manuscript

where f = (f 1 , ..., f Q ) with Q being the total number of patients and fq the normalized data for patient q.

During the above normalization step, no criteria were applied for outliers removal, since their presence in different features and for different patients can lead to a significant reduction of the already small dataset. Nevertheless, after performing the Phase 1, the preselected feature subset was carefully inspected to individually validate features' inclusion in the final evaluation step.

The final objective of this phase is to define optimal feature subsets for further processing. Once λ opt is determined, the features selected on step 4 are thus retained, and defined by the subset F 1 . Nevertheless, there may be λ values close to λ opt that produce even simpler models and with very similar CV errors. Variables selected for these other interesting values of λ can also be studied. Therefore, the models corresponding to those λ values between λ opt , denoted as M P r ncv , and λ min producing the simplest model M min ncv , represent the candidate models to be evaluated in Phase 2. Note that M min ncv can have one or more features, while P r represents the model holding all preselected features obtained with λ opt .

Phase 2: Building the classifier using a wrapper-based approach

In order to get a more robust final prediction, the final model used in the proposed BS classifier is determined in a wrapper fashion [START_REF] Kohavi | Wrappers for feature subset selection[END_REF] (see Fig. 3), through an ensemble classification system combining a set of classical, unregularized logistic regression models C, based on the preselected sparse models, M min ncv , ..., M P r ncv , C = {M min U nreg , ..., M P r U nreg }. Models C were trained using the stacked generalization (SG) approach [START_REF] Kuncheva | Combining Pattern Classifiers: Methods and Algorithms: Second Edition[END_REF]. SG is a way for combining different models or lower-level sub-models that have been trained for a particular classification task, allowing to derive higher-level models that minimize the generalization error rate, and thus the accuracy of the ensemble. In summary, the entire training set F L including only the P features of each model in C, is divided into j = 1, ..., J blocks, and each base classifier c j b is first trained on j-1 blocks of the training subset L j b , and internally evaluated on the hold out block T j b . Then, each trained fold classifier c j b (L j b ) is applied to the test set F T (it can be either the training set or an external test set) yielding the output o j b . Figure 3 illustrates this stacking procedure using 3 folds, although we set J=10 in our final design.

The final classification decision (right part of Fig. 3) is generated by combining the output of all base classifiers, by means of a fixed rule. In brief, the outputs of each individual base classifier o j b (soft outputs or posterior probabilities) were concatenated into the matrix O b , where

O b = [o 1 b , o 2 b , ..., o J b ]
, and serves as inputs for a fixed combiner. The proposed combiner is called the trimmed-weighted mean, which discards the worst and best among the lowlevel base classifiers, and calculates the weighted average with the remaining eight base classifiers [START_REF] Polikar | Ensemble based systems in decision making[END_REF]. The weight of each classifier, w j , is defined as its per-

Accepted manuscript

formance obtained on fold j not seen during training. Thus, a well-balanced average posterior probability vector is obtained as:

µ g (x) = 1 J J-2 j=1
w j o j b,g (x) [START_REF] Laguna | Techniques for Ventricular Repolarization Instability Assessment from the ECG[END_REF] where x represents the actual observation of the input database. This combiner calculates the support for a particular patient group, g= {symptomatic, asymptomatic}, by fusing all individual supports given by the retained base classifiers.

Model performance assessment

To assess the overall performance of the candidate predictive models, C, a repeated stratified holdout cross-validation strategy was used. The whole database was randomly splitted into training (L) and testing (T ) subsets using the ratio 70/30. L was used to train the ensemble classifiers as described in the previous section, resulting in E(L). Predictive performance was assessed by testing E(L) on T through the sensitivity (Se), specificity (Sp), and the area under the ROC curve (AU C) metrics. Since the population size is relatively small, the entire process was repeated N = 100 times in order to get an average estimate of the model performance and variance.

As mentioned before, the final model was determined in a wrapper fashion among those candidate subsets shortlisted during the regularization procedure. Therefore, the model producing the best metrics during the performance assessment was then used as our final BS classifier.

Results

The initial features dataset F of 198 features was first reduced by applying the elastic net procedure, thereby discarding the most irrelevant parameters and retaining a few potential candidate models.

Feature selection

Results obtained from the elastic net-based feature subsets selection and parameter tuning process are illustrated in Figure 4. The model hyperparameters were set as follows: α = 0.75 and λ varying between 0.001 and 0.1, logarithmically spaced in 40 steps. The darkest region observed in the heat map (see Fig. 4-a) corresponds to those λ values where CV errors are lower. By averaging such errors, the optimal λ value, λ opt , that produces the smallest mean error, is determined as indicated in Fig. 4-b (vertical dashed line in black). However, an alternative λ can be chosen if the error does not exceed a certain limit referring to that calculated for λ opt . A common practice is to apply the Accepted manuscript 1 standard error (SE) rule when CV is performed once [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Since this rule is no longer valid when using repeated CV because it reduces SE values as a result of averaging, we selected the most predictive model through a wrapperbased approach. In the proposed method, all λ values between λ opt and that one producing the simplest model (just one feature) were tested. The latter is represented in Fig. 4-b by the vertical dashed line in red.

In addition to CV errors, Fig. 4-b shows the number of features selected for each λ. Depending on λ, models with 1 and up to 30 features can be obtained with λ min ≥ λ ≥ λ opt (M min ncv , ..., M P r ncv ), representing a significant reduction (85%) from the initial number of variables in the dataset (P = 198).

The variables selected in each predictive model can be extracted from the coefficients' pathway generated as a function of λ during the regularization procedure. Table 2 lists the model (feature subsets) obtained for λ opt (P =30), with features ordered by importance. It shows that features related to ventricular depolarization were slightly higher in number (15 features) than those associated with repolarization (13 features). HRR-related indices deserve special attention as well, particularly when evaluated at later intervals of recovery such as at minutes 2 and 3. Of note, the importance of each feature within the model is given by its order of appearance in the list, derived from the coefficient's pathway.

On the other hand, features evaluated at specific phases (T1, T2, T3) of the test protocol were predominant (about two-thirds) as compared to those related to dynamic behaviors either during exercise or recovery. T1 (baseline phase) was the most important among the three studied phases for relevant repolarization features, while T2 and T3 for depolarization. Furthermore, recovery prevailed over exercise in terms of dynamic-related variables retained (9 features for recovery including HRR 2 and HRR 3 and 2 features for exercise). The overall performance achieved with this model, called Model 1 is included in Table 3, for both the training and testing sets.

Ensemble classifier performance

Fig. 5-a shows the average performance metrics obtained for all predictive models C = {M min U nreg , ..., M P r U nreg }, where M P r U nreg matches Model 1. As observed in Fig. 5-b, the model containing only 7 features, Model 2, was the most predictive among all wrapped candidate models, since it presented the smallest Euclidean distance (dist) between one unit and the average Se and Sp values.

The features retained for this simpler Model 2 correspond to the first seven features listed in Table 2, which are delimited by the horizontal dashed line. It is worth noting that these variables also appeared in Model 1, but in that case accompanied by other, less relevant features. Specifically, three local features related to repolarization period were preserved and only two from depolarization. HRR 2 and HRR 3 were also retained, confirming their relevance when predicting symptoms in BS. In addition, the first three features

Accepted manuscript

in Table 2 correspond to stage T1, highlighting the importance for the inclusion of basal measures. The average performance metrics of Model 2 and that of the model based on the first 3 features, Model 3, were also included in Table 3 for the training and testing subsets.

For Model 1, the results presented in Table 3 show a perfect predictive accuracy when tested with the training subset. However, when the test subset is used, Model 3 and in particular Model 2, achieved similar or even better results than Model 1 indicating that this simpler model is more generalizable and robust when applied to new, unseen patients.

Discussion

In this study, a robust predictive classifier has been implemented to discriminate between Brugada syndrome patients according to their symptomathology. To do this, ECG-derived markers that capture dynamic variations from ventricular repolarization and depolarization during physical exercise were assessed. These markers aimed at quantifying the temporal, spatial and morphological properties of the QRS complex and STT interval, since they might be affected by the electrophysiological disorders associated with this disease. Additionally, autonomic-related markers such as the HRR, estimated immediately after the maximum effort, was also considered.

First, a feature selection approach was implemented to efficiently reduce the initial feature dataset. It was achieved by a regularized optimization approach (elastic net). Then, a wrapper-based approach was used for final model evaluation. The first stage allowed to obtain sparse, predictive models while discarding irrelevant predictors [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Then, the unregularized version of these models underwent a second-round test to evaluate their predictive performance using a designed ensemble classifier. The model that presented the best balance between its overall performance and the number of predictors was selected to implement the final classifier.

The above-mentioned approach allowed to obtain a model of only 7 predictive features (Model 2) from a total of 198. From those features, the most important were related to the repolarization period, which characterized the morphological shape (T time ) and dispersion of the T wave (T pe , T pe /QT ), all of them evaluated at baseline stages. Two other features associated with the recovery process of the heart rate, HRR 2 and HRR 3 , were found to be relevant when predicting symptoms in BS, highlighting the importance of autonomic control in this pathology, specially the vagal modulation immediately after exercise [START_REF] Calvo | Multivariate classification of Brugada syndrome patients based on autonomic response to exercise testing[END_REF][START_REF] Romero | A biexponential approach for assessing parasympathetic reactivation after submaximal exercise[END_REF] or during sleep. The remaining predictors in the model were related to ventricular depolarization dynamics (∆S U V 2 RE , ∆S a V 2 RE ) evaluated during recovery. Thus, the results confirm the fact that Brugada syndrome is a complex multifactorial pathology, that needs to be addressed by considering factors related to the whole ventricular activity and the autonomic control, in an effort to obtain better predictive results [START_REF] Kawazoe | Risk stratification of ventricular fibrillation in Brugada syndrome using noninvasive scoring methods[END_REF][START_REF] Postema | Risk stratification in Brugada syndrome: Where is the finish line?[END_REF].

Accepted manuscript

The performance metrics achieved can be considered acceptable given the challenge of distinguishing these two groups, whose electrophysiological properties might be quite similar. Of note, these metrics were estimated by repeated cross-validation (hold out) 100 times, so as to obtain a realistic performance estimate when tested on new data, not seen during training. Also noteworthy is the novel configuration used in the design of the classifier, where through the fusion of several base classifiers, trained in different data subsets, a more robust prediction of the classes is obtained. Moreover, it should be noted that these results were obtained exclusively from electrocardiographic markers, without including any clinical parameters, as it is usually done in many clinical studies. Further work is directed towards the integration of such clinical data in a similar multivariate classifier.

In general, the results obtained in this study have shown the usefulness of the stress test to unmask the differences between symptomatic and asymptomatic BS patients. For example, we can highlight from the parameters included in the final model, Model 2, that the recovery of the heart rate is more informative when evaluated in later stages of this period, and not in the first minute as it is commonly reported [START_REF] Cole | Heart-rate recovery immediately after exercise as a predictor of mortality[END_REF]. In fact, in a previous study performed on the same population, we have demonstrated that the best timing for assessing significant differences among these patient groups through the HRR, was found between the minute 2.5 and 3 of recovery [START_REF] Romero | A biexponential approach for assessing parasympathetic reactivation after submaximal exercise[END_REF]. Another finding related to this period was that the two depolarization markers involved in the model, were associated with the second half of the QRS complex, and the dynamics of these markers. This suggests that both ventricular depolarization dynamics and vagal activity reactivation during that period are closely related to the symptomatology in this population. Finally and most importantly, the markers linked to ventricular repolarization appeared to be the most relevant to the model, which were tested separately in Model 3. The most remarkable in this case is that these particular markers were assessed at the beginning of the test, which gives them greater relevance because they could be measured at rest, avoiding exercise testing.

Quantifying risk level in asymptomatic Brugada patients is still a major clinical challenge. Better predictive markers are needed to improve prognosis and to optimize the therapy for a given patient. The analysis of ventricular depolarization and repolarization combined with measures related to autonomic control provides valuable information about the patients' condition and prognosis. In a study performed by Postema et al. [START_REF] Postema | Local Depolarization Abnormalities Are the Dominant Pathophysiologic Mechanism for Type 1 Electrocardiogram in Brugada Syndrome. A Study of Electrocardiograms, Vectorcardiograms, and Body Surface Potential Maps During Ajmaline Provocation[END_REF], repolarization abnormalities in Brugada syndrome were reported to be mostly induced by depolarization disorders. This was thought to be a result of the heterogeneity in the action potential duration with a ventricular endo-epi gradient. However, in Meijborg et al. [START_REF] Meijborg | Interventricular dispersion in repolarization causes bifid T waves in dogs with dofetilide-induced long QT syndrome[END_REF], repolarization abnormalities were thought to be related to an increased interventricular and LV-intraventricular dispersion in repolarization time, after dofetilide infusion in an experimental model of dofetilide-induced long QT syndrome type 2 (LQT2). Tokioka et al. [START_REF] Tokioka | Electrocardiographic parameters and fatal arrhythmic events in patients with brugada syndrome: Combination of depolarization and repolarization abnormalities[END_REF] reported that the combination of both repolarization and depolarization abnormalities enables potential identification of high-and low-risk Brugada patients. However, all these studies did Accepted manuscript not consider autonomic related markers in their analyses, which can provide an added value in the final outcomes. To our knowledge, the present work is the first in providing such quantitative comparison combining depolarization, repolarization and autonomic features for risk stratification in Brugada syndrome.

Limitations

According to results, despite the proposed classifier may seem potentially helpful to differentiate among Brugada syndrome patients, the relatively small sample population may have limited the impact of the results obtained. Moreover, the class imbalance existing between the two patient groups can also affect the prediction of the minority class in favor of the other. Therefore, all markers used in the model should be validated with larger clinical series to confirm its predictive potential. Finally, although the lack of follow-up information does not allow us to propose the use of these features as potential risk markers of future cardiac events, they remain as promising candidates that should be tested in other cohorts having such information.

Conclusions

The results obtained in this study confirm the significance of performing controlled autonomic manoeuvres, such as a physical exercise test, to elucidate the potential differences that may exist between patients with Brugada syndrome. Symptomatic and asymptomatic patients can produce different responses to this type of challenge, which can be detected by means of an exhaustive analysis of the electrophysiological properties and autonomic parameters that are likely to be affected by this pathology. The main conclusion of the study concerns the fact that both repolarization and depolarization characteristics measured at baseline and during exercise, are crucial in understanding the pathogenesis of this disease, and that the analysis of heart rate recovery at longer intervals may emerge as a sensitive risk marker. Accepted manuscript 2, respectively. Se: Sensitivity; Sp: Specificity; AUC: area under the ROC curve. Accepted manuscript

Accepted manuscript

Elastic net regularization

Accepted manuscript

Accepted manuscript

3 Fig. 1

 31 Fig. 1 An example of the HRR curve fitted on HR values during the recovery period of the exercise test. HRR 1 , HRR 2 , and HRR 3 represent the heart rate recovery measures evaluated at the first, second and third minute of recovery, respectively.

Fig. 2 FFig. 3

 23 Fig. 2 General work-flow including the regularized feature selection process (Phase 1) and the wrapper step (Phase 2) using the ensemble system. ncv: non cross-validated model; P r index of the preselected model obtained for λopt: Unreg: unregularized model; f : index of the actual wrapped model.

Fig. 4

 4 Fig. 4 Results from the elastic net regression of Phase 1, to find the most predictive features to be used in the ensemble classifier. (a) Deviance heat map obtained after 100 realizations (repeated CV: 100 times, 10-fold) using different λ values. (b) Mean ± SD of the deviance values displayed in a) (left axis), and the number of features selected for each λ value (right axis). Vertical dashed line in black indicates the optimal λ producing the lowest mean deviance.

Fig. 5 a

 5 Fig. 5 a) Average performance metrics obtained for the testing in all models C = {M min U nreg , ..., M P r U nreg }. b) Euclidean distance obtained for each model based on Se an Sp values. Model 1, Model 2 and Model 3, corresponds to all features, first 7 features, and first 3 features in Table2, respectively. Se: Sensitivity; Sp: Specificity; AUC: area under the ROC curve.

  The study protocol was approved by the respective local ethics committees: Comité d' Éthique du CHU de Rennes (ID RCB 2007-A00887-46), Comité d' Éthique du CHU Saint-Pierre, Comité d' Éthique du CHU de Nantes, Comité d' Éthique du CHU de Bordeaux, Comité d' Éthique du CHU de Brest and Comité d' Éthique du Centre Hospitalier de La Rochelle; in accordance with the ICH E6 recommendations, and the ethical principles of the Helsinki Declaration

Table 1

 1 Baseline clinical characteristics of the study population. is obtained using the Mann-Whitney test for continuous features and χ 2 test for categorical features.

	Clinical features	Symptomatic Asymptomatic	p-value
	Sex males	21 (84%)	61 (72%)	0.22
	Age (years)	46.8 ±15.5	44.0 ±13.1	0.63
	Symptoms
	Cardiac arrest	11 (44%)
	Syncope	14 (56%)
	Spontaneous Type-1	7 (28%)	27 (32%)	0.72
	ICD implanted	25 (100%)	19 (22%)	< 0.01
	SNC5A mutation (N=79)
	Positive	7 (28%)	22 (26%)
	Negative	13 (52%)	37 (44%)	0.85
	Max. HR (beats/min)	151.7 ±17.5	160.9 ±18.5	0.07
	Max. workload (Watts)	171.2 ±57.8	170.1 ±57.9	0.66
	p-value

Acknowledgments

This work was partly supported by a grant of the French Ministry of Health (Programme Hospitalier de Recherche Clinique -PHRC Regional). D. Romero acknowledges the financial support of the Fondation Lefoulon-Delalande, Institut de France, France. M. Calvo acknowledges the financial support of the social program funded by CaixaBank, Spain.

Conflict of interest

The authors declare that they have no conflict of interest.

Accepted manuscript

Table 2 Ordered subset of the 30 features selected by choosing the optimal λ value, λopt. The first seven features correspond to the smaller subset associated with the final selected model (Model 2). The numbers in the first column indicate the models in which each feature has been involved. T-wave intervals ratio in lead II during maximum effort 1

Model

Tarea II

T3

T-wave area's ratio in lead II when recovery ended

R-wave angle in lead V5 during maximum effort

Accepted manuscript