Capucine Tong

Amaury Delamarre

Romaric De Lépinau

Andrea Scaccabarozzi

Fabrice Oehler

Jean-Christophe Harmand

Stéphane Collin
email: stephane.collin@c2n.upsaclay.fr

Andrea Cattoni
email: andrea.cattoni@c2n.upsaclay.fr

GaAs/GaInP nanowire solar cell on Si with state-of-the-art V oc and quasi-Fermi level splitting

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

Tandem solar cells consisting of III-V semiconductor layers with optimal bandgap stacked on top of a conventional Si cell are one of the most promising ways to overcome the Shockley-Queisser efficiency limit for single-junction solar cells of ∼ 30 % (1 sun) [1]. The current power conversion efficiency record for Si-based multijunction solar cells is 35.9 % and was obtained by direct wafer bonding and substrate removal [2]. Mismatch in the lattice constant and thermal expansion coefficient between the materials renders the direct growth of high-quality III-V layers on Si difficult and leads to maximum efficiencies of 25.9 % [3].

Owing to their small dimensions, III-V nanowires (NWs) can accommodate both lattice and thermal coefficient mismatches without the introduction of misfit dislocations at the interface [4,5]. Furthermore, NWs absorb a higher amount of light than their planar counterparts due to light trapping phenomena [6][7][8]. Numerical simulations indicated a theoretical efficiency of 33.8 % for a tandem GaAsP NW solar cell (1.7 eV) on Si [9]. The simulations were carried out on a NW cell of composition GaAs 0.77 P 0.23 in a pi-n radial structure, taking advantage of the orthogonality between the direction of light absorption and carrier collection, even though the debate about the advantages and disadvantages of axial vs. radial NW junctions is still ongoing [10]. Given that ternary III-V compositions are difficult to control and less studied than binary compounds, most efforts in the photovoltaics community have been focused on fabricating GaAs-based NW solar cells [11][12][13][14][15].

To date, the best GaAs NW on Si tandem solar cell has reached an efficiency of 11.4 % [16], while the best single junction GaAs NW solar cell directly grown on Si has demonstrated an efficiency of 7.7 % [15]. So far, the performances of GaAs NW solar cells grown on Si have been primarily limited by their open-circuit voltage (V OC ¡ 0.62 V for NW array solar cells [15]). Note that in the same study, a V OC of 0.77 V was obtained for a single NW solar cell, however the direct comparison with NW array solar cells is difficult and should consider the fact that a single NW absorbs light from an area, the absorption crosssection, much larger than its geometrical cross-section [7]. The maximum V OC achievable is given by the quasi-Fermi level (qFl) splitting (∆µ) generated in the absorber material and related to the density of photogenerated charge carriers [17,18]. The qFl splitting can be determined from absolute photoluminescence (PL) measurements, and compared to the electrically measured V OC to gain information on the origin of voltage losses in the device [17][18][19][START_REF] Smestad | Solar Energy Materials and Solar Cells[END_REF][START_REF] Rumyantsev | Solar Energy Materials and Solar Cells[END_REF][START_REF] Bauer | [END_REF][23][24]. Other parameters influencing the performance of a NW solar cell include the NW crystal quality, vertical yield, doping concentration and surface passivation [25].

While a good NW crystal quality is essential to make a performant solar cell, it is common to find a high density of stacking faults in self-catalyzed GaAs NWs due to crystal phase mixing [26][27][28]. Ensuring a high yield of vertical NWs is also critical and several studies have been dedicated to its optimization [29][30][31][32][33][34][35][36][37]. Studies show that the various substrate preparations and growth parameters controlling the NW growth direction make it difficult to obtain reproducible results. However, studies have shown empirically that GaP NW arrays systematically present higher vertical yields than GaAs NW arrays [35,38], and that growing a GaP stem prior to GaAs NW growth can enhance the yield of vertical NWs [14,27,37,38]. Nevertheless, the presence of a valence band (VB) offset between the GaP and the Si may introduce an additional barrier for hole extraction, and its impact on NW solar cell performances has not yet been investigated.

Semiconductor doping is key to functional p-n junctions, however the introduction of dopants during NW growth greatly affects the thermodynamics and kinetics of the sys-tem often resulting in axial and radial doping gradients in NWs [39][40][41][42][43]. Doping compensation mechanisms can also prevent high doping concentrations to be achieved [44][45][46], rendering surface passivation essential for most III-V NWs. This is important for GaAs NWs as they are particularly prone to surface depletion due to high surface recombination velocity [47][48][49][50][51][52].

Finally, all NWs need to be ohmically contacted to efficiently extract charge carriers from the solar cell. Indium tin oxide (ITO) is often used as a transparent conductive oxide (TCO) to collect the electrons while minimizing shadowing effects. But its interface with GaAs or GaInP has been shown to form an important Schottky barrier [53][54][55][56][START_REF] Lepinau | 47th IEEE Photovoltaic Specialists Conference (PVSC)[END_REF], hindering the overall device performance.

In this work, we present a GaAs/GaInP NW solar cell directly grown on Si by molecular beam epitaxy (MBE). Current-voltage measurements of the cell indicate a V OC of 0.65 V, which is at the state-of-the-art for a NW array solar cell on Si. Hyperspectral photoluminescence measurements evidence a qFl splitting of 0.84 eV, significantly higher than qV OC , indicating good optoelectronic quality NWs. The various reasons behind the lower qV OC are discussed. Optical and electronic simulations are used to study the main absorption and electrical losses. In this study, one-sixth of a two-inch p-doped (B) Cz Si(111) wafer supplied by Siltronix is used (nominal resistivity: 0.010.02 Ω cm). The surface is cleaned by three successive cycles of oxidation in oxygen plasma reactive ion etching (RIE) and etching in hydrofluoric acid (HF) diluted in water to 5 % concentration. A 37 1 nm thick silica layer is deposited via plasma-enhanced chemical vapor deposition (PECVD). The substrates are coated with an electron-sensitive resist (495PMMA, 2 %) and patterned in a Vistec EBPG 5000+ electron beam lithography machine, operating at 100 kV, resulting in a hexagonal array of holes with a 500 nm pitch. This hole pattern is transferred into the SiOx layer by capacitively coupled plasma-reactive ion etching (CCP-RIE) at 7 mTorr, with a mixture of SF 6 and CHF 3 gases (flow rate ratio of 8:20) and a bias voltage of 190 V. To avoid plasma-induced damage and amorphization of the first Si layers, about 4 nm of SiOx is left inside the holes. The PMMA resist is removed by successive immersion in trichloroethylene and acetone with sonication, then rinsed with isopropanol. To ensure complete removal of organic residues, the sample is then subjected to an oxygen plasma cleaning by CCP-RIE. The remaining SiOx inside the holes is removed by a dip in a 1 % HF solution (diluted in water), during 30 s. After rinsing with de-ionized water and drying with N 2 , the sample is stored overnight in a box filled with N 2 before being transported the next day to a Riber32 MBE chamber. The substrate is outgassed in the preparation chamber at 500-550 C until a pressure of 3 x 10 -9 Torr is reached. Prior to growth, the substrate is held at 690 10 C for 10 min in the MBE growth chamber on a rotating holder (typically 10 rpm), after which the growth temperature is stabilized at 620 ± 5 C as monitored with a 0.9 µm IR-CON pyrometer. The growth procedure is initiated by a 1 min Ga pre-deposition with an equivalent thickness of 42 monolayer (ML) GaAs(100), and a Be flux of 6 x 10 11 cm -2 s -1 . A Be-doped GaP stem (p-GaP) is grown during 80 s at a P 2 /Ga flux ratio of 1, before linearly closing the P valve while opening As over 30 s. This results in a linearly graded GaP/GaAs transition which contributes to flatten the VB offset between GaP and GaAs [14]. The main Bedoped GaAs core segment (p-GaAs) is then grown under an As 4 /Ga flux ratio of 1.2, and the same Be flux during 12 min. The Ga flux during the GaP stem and GaAs core growths is set to correspond to a 2 Å s -1 growth rate on GaAs(100). The liquid Ga catalyst droplet is then consumed in situ during 20 min at the same temperature by shutting down the Ga flux and using an As 4 flux at 16 % of its value during core growth. A nominally intrinsic GaAs shell (i-GaAs) is grown during 30 min at a substrate temperature of 590 ± 10 C under a Ga flux of 0.74 Å s -1 and an excess of As 4 (As 4 /Ga flux ratio of 7.1). The substrate temperature is finally set to 470 ± 10 C for the 25 min growth of a Si-doped GaInP shell (n-GaInP), aiming at a planar growth rate of 1.1 Å s -1 , a planar composition of In 0.58 Ga 0.42 P and a planar doping level of 6 x 10 18 cm -3 . The NW structure is terminated by an n-type delta-doped GaInP capping layer (n + -GaInP). GaInP is grown during 32 s before stopping the Ga and In flux and providing an increased Si flux under excess P 2 for another 32 s, resulting in a Si surface concentration of 1.4 x 10 12 cm -2 in one atomic plane. This cycle is repeated 10 times, resulting in an equivalent GaInP layer grown during 5.35 min with an average Si concentration of 1.7 x 10 19 cm -3 .

Characterization of as-grown NW arrays.

From the scanning electron microscope (SEM) image of the as-grown NW array (Figure 1b), the NWs height and diameter could be determined as h NW = 1800 130 nm and D NW = 220 20 nm, respectively. The p-GaAs NW core diameter reported in Figure 1a was obtained from the measurement of NW cores previously grown in nearly identical conditions. The i-GaAs shell thickness was deduced by subtracting the p-GaAs NW core diameter from the total width of the As-containing region, as determined by cross section Energy-Dispersive X-ray spectroscopy (EDX), see Figure 6 in the appendix. Similarly, the GaInP shell thickness was determined from the lateral extension of the (In,P)containing regions of the cross-section EDX maps. Note that these EDX-determined shell thicknesses are rough estimates and significant errors are expected due to several factors. First, the diffusion dynamics of group III adatoms vary depending on the species and the growth temperature. Some regions in the shell may collect more or less material due to shadowing of surroundings NWs, resulting in non-uniform thickness. Furthermore, the NW core diameter is also found to vary between growths and after droplet consumption (typically within 20 nm). Finally, the geometrical surface orientation relative to the angle of the MBE cell leads to severe growth rate differences between the NW head and sidewall, as much as a factor of six, as measured by Transmission Electron Microscopy (TEM), see Figure 6 in the appendix. On average the GaInP NW sidewalls are found to be richer in Ga, which is why we correct the nominal target Ga fraction to Ga 0.42 In 0.58 P, instead of using the theoretical GaAs-matched composition of Ga 0.51 In 0.49 P. Quantitative EDX analyses of the Ga content variation in the GaInP shell can be found in Figure 7 in the appendix. Extended defects such as dislocations and twins are also visible on the scanning TEM image in Figure 6 in the appendix.

The doping levels are estimated in the core and shell of previous GaAs NW samples grown in similar conditions using cathodoluminescence (CL) spectroscopy by fitting the spectrum to a generalized Planck law, as described in our earlier work [46,[START_REF] Chen | [END_REF]. The core p-type doping with Be is thus expected here to be close to p = 5 10 18 cm -3 , whereas the shell n-type doping with Si is limited to n ¡ 1 10 18 cm -3 [46]. The size and carrier concentration of the p-GaP stem are difficult to measure. Estimations from growth parameters suggest a roughly 80 nm long and 40 nm wide p-GaP stem with a doping concentration similar to the one in the p-GaAs core. Potential hole extraction barriers in the processed NW device arising from the introduction of this ad-ditional heterojunction at the foot of the NW are discussed in section 5.

B. Nanowire solar cell fabrication

The as-grown NWs are embedded in a BCB matrix (Cyclotene 3022-46) for planarization. Plasma CCP-RIE etching and a dip in 1 % HF follow to reveal the NW heads. Photolithography is then used to deposit 150 nm of ITO and Ti/Au contacts on specific areas of the sample to electrically connect NWs of the same cell and ensure good mechanical and electrical contact with measurement probes. The back contact is formed by depositing a 300 nm Al layer on the p-Si(111) wafer. The final device is made up of several circular cells with diameters ranging from 140 m to 1 mm. The champion cell, on which this article will be focused, is shown in Figure 2a. It has an outer diameter of 600 m and an inner diameter of 440 µm. After characterization, a cross section view of the same cell is obtained by Focused Ion Beam (FIB) milling and tilted SEM imaging (Figure 2b).

III. Optoelectronic properties

A. Current-voltage (JV) characteristics

The electrical performances of 13 NW solar cells are measured using a Keithley 2601B source measure unit in a two-probe configuration at room temperature, in the dark and under a 1-sun illumination. The latter is provided by a Newport 94011A-ES solar simulator equipped with an AM1.5G filter. The recorded open-circuit voltage (V OC ), short-circuit current (J SC ), fill factor (FF) and efficiency values are summarized in the appendix (Figure 8). The JV characteristics of the champion cell in the dark and under a 1-sun illumination are shown in Figure 3a. A high V OC of 0.65 V is obtained, along with a J SC of 15.5 mA/cm 2 and FF of 40 %. External quantum efficiency (EQE) integration across the relevant wavelength range (Figure 3c) resulted however in a J SC value of 14.2 mA/cm 2 . The latter will be the official J SC value reported as J SC determination from EQE has the advantage of being independent of the spectral shape of the light source, unlike J SC determination from J-V measurement. The resulting power conversion efficiency of 3.68 % is at the state-of-the-art of solar cells with GaAs-based radial junction NWs grown on Si [14,27]. However, the performances of this NW solar cell are still far behind the ones of the best GaAs thin film solar cells. The cross-over between the dark and the illuminated curves at 0.7 V suggests an electron barrier in the conduction band which is reduced upon photodoping. The slope of the illuminated JV curve near the short-circuit current point indicates either an incomplete collection of charge carriers, or a low shunt resistance which reduces the amount of current flowing through the solar cell junction. The semilogarithmic plot of the JV curve in the dark is shown in Figure 3b along with a one-diode model fit. A high and unfittable reverse-current can be observed, while the forward current is accurately fitted with a high ideality factor of 4 and a recombination parameter of J 0 = 4.1 10 -6 A cm -2 . Satisfying shunt and series resistance of 7.2 10 3 Ω cm 2 and 0.17 Ω cm 2 , respectively, were also fitted from the forward current.

B. External quantum efficiency

A custom-built setup is used to measure the EQE of the NW solar cell. Light emitted by a Xenon lamp is filtered using a grating monochromator and guided to the setup optical axis with an optical fibre. The photon flux over the spectral range of interest is first calibrated with a reference Si diode with known responsivity. Then, the sample is subject to the same illumination conditions, while the current at its terminals is collected at zero bias with a Keithley 2600 source measure unit. The resulting EQE of the NW cell is shown in Figure 3c. Between 700 and 800 nm the EQE remains relatively constant around 60 % and is higher than the EQE values measured in other radial junction GaAs-based NW solar cells [14,27]. Below 700 nm the EQE decreases steadily, although experiments showed that a high-temperature annealing could increase the short wavelength part of the EQE by improving contact conductivity [38]. Due to light scattering, specular reflectivity measurements could not be used to explicitly differentiate absorption issues from carrier collection issues, so simulations of the absorption in the different regions of the NW solar cell are carried out.

C. NW absorption modelling.

A frequency-domain modal method known as the Rigorous Coupled-Wave Analysis (RCWA) [START_REF] Hugonin | RETICOLO CODE 1D for the diffraction by stacks of lamellar 1D gratings (classical diffraction)[END_REF][START_REF] Lalanne | [END_REF][61][62] is used to model the absorption in the GaAs NWs, the GaInP contact and the ITO layer. A top view of the unit cell (hexagonal lattice) used to generate the NW array solar cell in the simulations is shown in the appendix in Figure 9 (left). A cross-section view of the simulated NW solar cell along a single NW is also shown in Figure 9 (right).

The simplified design incorporates the main elements of the NW solar cell previously presented, with the GaInP segment only present at the top of the NW as it is six times thicker than the GaInP shell and is expected to be the main source of photon absorption in GaInP. A more refined design is not justified also due to the morphology dispersion between the different NWs of the NW array composing the cell.

The absorption in each material (Figure 3c, dashed curves) is calculated using 20 Fourier orders and the optical constants used for the RCWA simulations are plotted in Figure 10 in the appendix. The simulated absorption in the GaAs NWs (brown) is slightly higher than the experimentally measured EQE of the solar cell (red), although the general trend of the EQE is well reproduced by the simulations. The calculated short-circuit current, assuming that all the photo-generated carriers are collected, is 18.77 mA/cm 2 . The difference of about 4.5 mA/cm 2 may be due to a non-perfect carrier collection or to an underestimation of the parasitic absorption inside the GaInP contact and/or the ITO. The parasitic absorption inside the GaInP contact is relatively high (2.84 mA/cm 2 ) even without considering the 34 nm-thick GaInP shell in the simulations.

The performances of the NW solar cell presented could be hindered by poor crystalline quality NWs, defective interfaces and contacts both at the foot (interface between the GaP stem and the Si substrate) and at the top of the NWs (contact between the GaInP shell and the ITO). As high optoelectronic quality NWs are essential for a performant solar cell, hyperspectral PL is used to assess the quasi-Fermi level splitting developed in the NWs under illumination. The latter has been shown to represent the maximum V OC that a solar cell can achieve [17,18], and is discussed in the following section.

IV. Absolute photoluminescence (PL) and quasi-Fermi level (qFl) splitting

To better assess the optical properties of the NWs, a spatially and spectrally resolved hyperspectral imager set-up calibrated in absolute units is used [17]. LED excitation at 450 nm is carried out through a microscope objective with numerical aperture (NA) 0.6. Measurements are performed at illumination intensities varying from below 1 sun to 81 suns. The measurement at 12 suns is carried out with a 532 nm laser due to a set-up change.

According to the generalized Planck law, the absolute photon flux (φ PL ) emitted from a semiconductor under illumination into a hemisphere above its flat surface can be related to the carrier temperature (T) and the qFl splitting (∆µ), assuming the latter is constant (equation 1). The qFl splitting of a material has previously been shown to correspond to the maximum achievable V OC of the device [17,18]. Here, A(hω) is the absorption probability of an incident photon of energy hω on a surface element, h the reduced Planck constant, c the speed of light in vacuum and k B the Boltzmann constant. The absolute photon flux φ PL is expressed as:

φ PL = A(hω) 4π 2 h3 c 2 (hω) 2 exp( (hω)-∆µ k B T ) -1 . (1)
This equation assumes an entirely collected Lambertian emission profile which has been shown to be an accurate description of the angular emission profile of nonoptimized NW arrays [63]. In non-degenerate semiconductors where the Fermi levels lie a few kT away from the band edges, the Fermi distribution can be approximated by a Boltzmann distribution:

φ PL = EQE(hω) (hω) 2 4π 2 h3 c 2 exp( -hω k B T )exp( ∆µ k B T ). (2)
In this equation, the absorption probability A(hω) has been approximated by the EQE(hω), since only the latter can be experimentally measured. We consider this approximation to be reasonable considering the good match between the EQE and the modeled A(hω) in Figure 3c. The semi-logarithmic plot in Figure 4a shows the absolute photoluminescence spectra of the cell under a 12 suns and 81 suns illumination (dashed curves) after averaging the hyperspectral PL images over an area of roughly 50 x 50 µm 2 . Figures 11a and11b in the appendix present the hyperspectral PL maps of the champion cell emission at 870 nm, under an illumination of 12 suns and 81 suns, respectively. Figures 11c and11d depict the qFl splitting maps of the same areas obtained by fitting the PL spectra using the generalized Planck law (equation 2) between 1.43 and 1.52 eV, with a fixed carrier temperature of 330 K. Quasi-Fermi level splitting variations of less than ± 30 meV can be observed across the maps, demonstrating the relatively high qFl splitting homogeneity over the characterized device surface. The solid curves in Figure 4a correspond to the PL spectra divided by the cell's EQE (from Figure 3c), while the red lines are their generalized Planck law fit. The ∆µ obtained under a 12 suns and 81 suns illumination are 0.95 eV and 1.01 eV, respectively. These two data points are indicated in light blue and cyan, respectively, in Figure 4b.

Due to lower signal-to-noise ratio at weaker illumination intensities (¡ 12 suns), instead of directly measuring spectrally resolved images we measured simple luminescence images, that integrate all wavelengths above the 850 nm long-pass filter, therefore presenting more intense signals. This way we obtain the relative variation of the luminescence from the 81 suns illumination down to 0.7 sun, which is proportional to exp( ∆µ k B T ) according to equation 2. Since ∆µ was determined to be 1.01 eV at 81 suns, ∆µ for lower illuminations can be deduced, as displayed in Figure 4b (dark blue squares). After interpolation, a qFl splitting of 0.84 eV is determined for a 1 sun illumination. This qFl splitting value is considerably higher than the measured qV OC at 1 sun (black squares in Figure 4b), and to the best of our knowledge constitutes a first-time quantification of the qFl splitting in NW array solar cells. Furthermore, as here the qFl splittings are determined by using a generalized Planck law fit in which the absorption probability A(hω) is approximated by the EQE(hω) (equation 2), the resulting qFl splitting values are lower limits. Upper limits of the qFl splittings can be deduced by using a 30 % higher absorption, as obtained from the optical simulations in section 3.3. The latter would translate into an 8 meV increase in the qFl splitting values (kT ln(1.3), where T = 330 K). The ideality factors n ∆µ and n V OC indicated in Figure 4b are obtained by fitting the slope of the variation of ∆µ and V OC , respectively, with illumination intensity and they are discussed at the end of this section. The high qFl splitting of 0.84 eV at 1 sun indicates a better optoelectronic quality than what could be deduced from electrical measurements and TEM observations alone. It also evidences significant room for device improvement, which can only be achieved if the mismatch between the qFl splitting and the qV OC is clarified.

A first reason could be inhomogeneous NW properties. Let us consider a situation where the NWs can be grouped into two distinct families: family 1 and family 2. In family 1, the NWs sustain a qFl splitting ∆µ 1 , and emit a luminescence signal Aφ bb exp( ∆µ 1 k B T ). In family 2, the NWs present a negligible qFl splitting ∆µ 2 , and emit a negligible luminescence. Such a situation relies on two assumptions: (i) The effective carrier lifetimes of both NW families are significantly different. (ii) Resistances or current barriers must be present in between NWs, which could originate from processing issues, e.g. local incomplete BCB removal or inhomogeneous TCO deposition. Otherwise, the qFl splitting of each individual NW would equal its voltage, and the voltage of both NW families would equilibrate. It is possible that points (i) and (ii) are true for our solar cells, as they have been observed in similar NWs (VLS grown GaAs NWs on Si) by atomic force microscopy (AFM) [64] and electron beam induced current (EBIC) measurements [65], and in GaN/InGaN core-shell NWs by scanning photocurrent microscopy (SPCM) [66]. Therefore, ∆µ 1 and ∆µ 2 can differ, as well as their luminescence emissions. Nevertheless, they cannot be distinguished by luminescence imaging, since the NW-to-NW distance is in the order of the diffraction limit (500 nm). The luminescence density recorded on an area comprising an equal number of NWs from families 1 and 2 will be an average of their respective emission, i.e. Aφ bb exp(

∆µ 1 k B T ) × 0.5 ≈ Aφ bb exp( ∆µ 1 -20 meV k B T
). This means that the determined qFl splitting (here ∆µ 1 -20 meV) is close to its highest value. Simultaneously, the voltage can take any value in between ∆µ 1 q and ∆µ 2 q , depending on the local resistances. In the extreme situation where NWs 1 are totally disconnected and NWs 2 are connected without resistance, we have V = ∆µ 2 q . In the opposite case (NWs 1 connected, NWs 2 disconnected), V = ∆µ 1 q . The fact that the qFl splitting obtained (0.84 eV) is higher than the measured qV OC of the cell (V OC = 0.65 V, qV OC = 0.65 eV) could also be due to additional recombinations at heterointerfaces. This has been observed in other types of solar cells, e.g. CIS [67], CIGS [68,69] or perovskites [70]. In our case, the efficient extraction of the majority charge carriers may be prevented by the p-Si/p-GaP/p-GaAs heterointerface at the foot and by the i-GaAs/n-GaInP/ITO heterointerface at the top, and they may recombine at the dislocations near the i-GaAs/n-GaInP interface introduced due to GaInP shell composition variation along the NW axis. More details on charge generation, separation, recombination and collection could be obtained by using for example two-photon induced photocurrent and photocurrent-voltage responses [71].

The difference in the ideality factors obtained from the variation of ∆µ and V OC , respectively, with illumination intensity (for more information on the ideality factors' determination please refer to page 51 of the book by D. Abou-Ras et al. [START_REF] Abou-Ras | Advanced Characterization Techniques for Thin Film Solar Cells: Second Edition[END_REF]) could also be explained by inhomogeneous NW properties and/or the presence of potential barriers at interfaces and recombination currents. As the ideality factor provides insight on the location of the dominating recombination mechanism in a solar cell [START_REF] Scheer | Chalcogenide Photovoltaics[END_REF], the higher ideality factor n V OC = 2.6 compared to the optical ideality factor n ∆µ = 1.4 suggests the presence of an additional recombination channel at the GaAs absorber interfaces. However, these ideality factors are still lower than the one obtained from a one-diode model fit of the dark currentvoltage characteristics of the cell (n = 4, Figure 3b), and further studies are required to elucidate this matter.

V. Modelling hole transport at the p-type contact

The GaP stem introduced at the p-type contact of the NW structure (brown in Figure 1a) to enhance the vertical yield of the NWs is an additional heterostructure which could be a barrier for hole extraction [START_REF] Gudovskikh | [END_REF]. To study the impact The light blue and cyan squares correspond to qFl splittings obtained from spectrally resolved luminescence images, while the dark blue squares are obtained from simple images that integrate all wavelengths. The record V OC for planar GaAs cell (red star) corresponds to the one obtained by Alta Devices [START_REF] Kayes | IEEE Photovoltaic Specialists Conference[END_REF] of this barrier, we simulated a simplified one-dimensional structure of the cell previously presented, using the solar cell device simulator SCAPS [76]. Since charges in the NWs are separated radially, the modelled structure illustrated at the top of Figure 5a first follows the NW axis and then the direction perpendicular to the NW axis. Figure 5a presents the simulated band diagrams of the modelled structure under a 1 sun illumination and an applied bias corresponding to the maximum power point (V MPP ), for different acceptor doping concentrations in the GaP stem, ranging from 1 × 10 18 cm -3 to 8 × 10 18 cm -3 . All band diagrams have been shifted such that the conduction bands (Ec) are aligned on the n-GaInP side. The parameters used in the SCAPS simulations are detailed in Table I in the appendix.

The VB offset (∆E V ) at the GaP/Si interface directly impacts the simulation results, however its experimentally determined value varies significantly in the literature (from 0.6 eV [77] to 1.1 eV [78]). An average ∆E V of 0.9 eV was used in the simulations presented in this section (corresponding to a GaP electron affinity of 3.8 eV [START_REF]Characteristics and Properties[END_REF]). The simulation results obtained using the lower and upper bound VB offset values found in the literature are shown in Figure 12 in the appendix. As observed in Figure 5a, the increase of acceptor doping concentration in GaP brings about a narrowing of the VB potential barrier there. This is also evidenced by the current-voltage curves in Figure 5b through the improved V OC , fill factor, and hence efficiency, with increased dop-

p-Si(111) p-GaAs n-GaInP i-GaAs p-GaP n + -GaInP p-Si(111) p-GaP a) b) c)
// to NW axis ⊥ to NW axis 0.09 0.10 0.11 -

Hole Current Density (mA/cm 2 ) Distance (µm)

Tunneling (p=1e18cm -3 ) Total (p=1e18cm -3 ) Tunneling (p=3e18cm -3 ) Total (p=3e18cm -3 ) Tunneling (p=5e18cm -3 ) Total(p=5e18cm -3 ) Tunneling (p=8e18cm Increasing GaP doping p = 1e18 cm -3 p = 3e18 cm -3 p = 5e18 cm -3 p = 8e18 cm -3 p = 1e18 cm -3 (0.96 V, 15.0 mA/cm 2 , 77.0 %, 11.1 %) p = 3e18 cm -3 (1.03 V, 15.0 mA/cm 2 , 78.0 %, 12.1 %) p = 5e18 cm -3 (1.08 V, 15.0 mA/cm 2 , 78.7 %, 12.8 %) p = 8e18 cm -3 (1.12 V, 15.0 mA/cm 2 , 80.6 %, 13.6 %) ing concentration. In fact, as the potential barrier width in GaP is decreased with increasing doping concentration, hole tunnelling becomes more efficient. This is illustrated in Figure 5c where the current density of intra-band tunnelling holes at V MPP (grey/black curves) is almost doubled at the GaP/Si interface when the doping concentration changes from 1 × 10 18 cm -3 (light grey) to 8 × 10 18 cm -3 (black). For a doping concentration above 5 × 10 18 cm -3 , the hole tunnelling current densities (two darkest curves) become comparable to the total hole current density (red lines), indicating that a high doping concentration can help overcome the VB potential barrier in GaP thanks to intraband tunnelling. The more efficient transport of holes via intra-band tunnelling can also be observed in the band diagrams in Figure 5a, through a lower potential drop across the GaP/Si interface. The doping concentration required for holes to predominantly tunnel through the potential barrier depends on the GaP/Si VB offset (Figure 12 in the appendix), making it difficult to determine a target threshold doping concentration. If the GaP/Si VB offset is closer to the lower bound value of 0.6 eV from the literature [77], then simulations show that the solar cell's performances will not vary much when the GaP doping concentration is increased from 1 × 10 18 cm -3 to 8 × 10 18 cm -3 . However, the greater the GaP/Si VB offset the more the cell's performances will be impacted by the GaP doping concentration, as seen in Figures 12c and12d in the appendix.

It is difficult to precisely determine the size of the actual GaP stem in the NW structure and its doping concentration. However, previous CL studies [46] have shown that the Be doping concentration in our GaAs NW cores, which are grown via the same VLS mechanism, can be superior to 5 × 10 18 cm -3 . It is thus reasonable to assume that the doping level in the GaP stem is similar and hence that with a GaP/Si VB offset between 0.6 -1.1 eV, the potential barrier in the VB does not significantly impact hole extraction at the p-type contact in our devices.

VI. Conclusions

A core-shell GaAs/GaInP NW solar cell directly grown on Si by MBE was successfully fabricated. The cell demonstrated a PCE of almost 3.7 % and a state-of-the-art V OC of 0.65 V. A quasi-Fermi level splitting of 0.84 eV at 1 sun (up to 1.01 eV at 81 suns) was determined from hyperspectral PL measurements, which is considerably higher than qV OC . This indicates NWs of good optoelectronic quality and significant room for the improvement of our devices. Numerical 1D simulations showed that the potential barrier at the p-type contact can be overcome if the GaP stem is sufficiently doped. As we expect the hole concentration in the GaP stem to reach 5 × 10 18 cm -3 , we are confident that our NW solar cell's performances can be further enhanced by improvements at the n-type contact (GaInP/ITO).

  area growth of III-V nanowire arrays on Si 1. Patterned substrate preparations.

Figure 1 Figure 2

 12 Figure 1 (a) Schematic of the dimensions of a core-shell GaAs/GaInP heterojunction NW with a delta-doped GaInP capping layer. (b) Tilted (35 • ) SEM overview of an as-grownGaAs/GaInP NW array with a 500 nm pitch.

Figure 3 Figure 4

 34 Figure3(a) JV curves in the dark (blue) and under a 1 sun illumination (red), (b) semi-logarithmic plot of the JV curve in the dark (blue), along with a one-diode fit (black). (c) EQE of the GaAs/GaInP NW solar cell (red) and simulated absorptions for GaAs, GaInP and ITO (dashed) using the simplified design shown in Figure9in the appendix.

Figure 5 (

 5 Figure 5 (a) Band diagrams obtained from 1D SCAPS simulations under 1 sun illumination and V MPP bias for different GaP acceptor doping concentrations. (b) Current-voltage characteristics (values in brackets correspond to V OC , J SC , FF, and η) and (c) hole current densities obtained from the same SCAPS simulations.

Figure 6 (Figure 7 Figure 8 Figure 9

 6789 Figure 6 (left) Scanning transmission electron microscopy (STEM) image taken with a high-angle annular dark field (HAADF) detector. The white arrow indicates the GaAs [111]B growth direction. (right) Corresponding In, Ga, As, and P elemental maps of the final NW section, indicated by the red square in the STEM-HAADF image, obtained by energy dispersive X-ray analysis (EDX). The core-shell heterostructure can be seen from the 112 zone-axis.

Figure 10 Figure 11

 1011 Figure 10 Refractive indices (n) and extinction coefficients (k)used in absorption simulations carried out using the RCWA method.

χ

  Si = 4.05 eV χ GaP = 3.80 eV χ Si = 4.05 eV χ GaP = 3.51 eV χ Si = 4.05 eV χ GaP = 4

Figure 12

 12 Figure 12 Impact of the VB offset (∆E V ) at the GaP/Si interface on the SCAPS simulation results. (a) Comparison between the lower and upper bound GaP/Si VB offset values found in the literature (∆E V =0.6 eV[77] and ∆E V =1.1 eV[78], respectively) and the average value of 0.9 eV used in the simulations presented in the main article. The corresponding electron affinities for Si and GaP (χ Si and χ GaP ) are also indicated. (b) Band diagrams obtained from 1D SCAPS simulations under 1 sun illumination and V MPP bias for different GaP acceptor doping concentrations. (c) Corresponding current-voltage characteristics under 1 sun illumination. The figures of merit V OC , J SC , FF, and η are summarized in table (d).

Table I

 I Parameters used in the 1D SCAPS simulations presented in Figure5.

				n + -GaInP	6	1.85	4.01	11.8	6.24 x 10 17
				n-GaInP	28	1.85	4.01	11.8	6.24 x 10 17
				i-GaAs	24	1.424	4.07	12.9	4.7 x 10 17
				p-GaAs	50	1.424	4.07	12.9	4.7 x 10 17
	The material parameters were extracted from references [79-82]. The acceptor	doping concentration in p-GaP (in bold in the table) was varied between 1 x 10 18	and 8 x 10 18 cm -3 in the SCAPS simulations.	p-GaP Graded GaP/GaAs p-Si	58 22 100 Thickness (nm)	2.26 linearly graded 1.12 Bandgap (eV)	3.8 linearly graded 4.05 Electron affinity (eV)	11.1 linearly graded 11.9 Dielectric permittivity	1.8 x 10 19 linearly graded 2.8 x 10 19 CB effective density of states (cm -3 )	VB effective density of states	(cm -3

Acknowledgements

The authors thank G. Patriarche for STEM-EDX images, L. Largeau for cross-sectional sample preparation by FIB, and Ombeline Lafont, Baptiste Brenguier and Hao Xu for their contribution to early photoluminescence experiments. The authors also thank C. Cornet for insights on the GaP/Si valence band offset values.

The authors acknowledge financial support from the French RENATECH network and the French government in the framework of the Programme dInvestissement dAvenir (ANR-IEED-002-01) and the ANR projects Nanocell (ANR-15-CE05-0026) and PIANIST (ANR-21-CE09-0020-03). This work was also partly funded in the frame of the EMPIR 19ENG05 NanoWires project. The EMPIR (European Metrology Programme for Innovation and Research) initiative is co-funded by the European Union's Horizon 2020 research and innovation programme and the EMPIR Participating States.

Conflicts of interest

The authors declare no conflicts of interests.

D. Hyperspectral PL and qFl splitting maps

From Figures 11c and11d, quasi-Fermi level splitting variations of ± 20 meV and ± 25 meV can be observed respectively for the 12 suns and 81 suns PL measurements. This demonstrates the relatively high qFl splitting homogeneity over the characterized device surface.


Eugenio Clerico
email: clerico@stats.ox.ac.uk

George Deligiannidis

Benjamin Guedj

Arnaud Doucet

A PAC-Bayes bound for deterministic classifiers

We establish a disintegrated PAC-Bayesian bound, for classifiers that are trained via continuous-time (non-stochastic) gradient descent. Contrarily to what is standard in the PAC-Bayesian setting, our result applies to a training algorithm that is deterministic, conditioned on a random initialisation, without requiring any de-randomisation step. We provide a broad discussion of the main features of the bound that we propose, and we study analytically and empirically its behaviour on linear models, finding promising results.

Introduction

Effectively upper-bounding the generalisation error of modern learning algorithms is one of the main open challenges for the statistical learning theory community [START_REF] Zhang | Understanding deep learning requires rethinking generalization[END_REF]. Originally, properties of the hypothesis space, such as VC dimension and Rademacher complexity [START_REF] Bousquet | Introduction to Statistical Learning Theory[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] Shalev-Shwartz | Understanding Machine Learning -From Theory to Algorithms[END_REF], were used to establish worst-case generalisation bounds, holding uniformly over all the possible algorithms and datasets. However, these results are often vacuous when applied to the current over-parameterised models, the modern perspective focuses on algorithmand data-dependent bounds [START_REF] Mcallester | Some PAC-Bayesian theorems[END_REF][START_REF] Bousquet | Stability and generalization[END_REF][START_REF] Dwork | The algorithmic foundations of differential privacy[END_REF][START_REF] Hardt | Train faster, generalize better: Stability of stochastic gradient descent[END_REF][START_REF] Xu | Information-theoretic analysis of generalization capability of learning algorithms[END_REF]Clerico et al., 2022b;[START_REF] Lugosi | Generalization bounds via convex analysis[END_REF]. Among the various approaches proposed in the literature, the PAC-Bayesian framework [START_REF] Guedj | A primer on PAC-Bayesian learning[END_REF][START_REF] Alquier | User-friendly introduction to PAC-Bayes bounds[END_REF] has obtained promising empirical results for deep networks [START_REF] Dziugaite | Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data[END_REF][START_REF] Zhou | Non-vacuous generalization bounds at the ImageNet scale: a PAC-Bayesian compression approach[END_REF][START_REF] Pérez-Ortiz | Tighter risk certificates for neural networks[END_REF]Clerico et al., 2022a).

Typically, the PAC-Bayesian bounds are upper bounds on the expected population loss of stochastic classifiers, holding with high probability on the random draw of the training dataset. This setting, where the model's parameters are required to be intrinsically random, does not match with the more standard scenario of a non-randomised predictor. Although most available empirical PAC-Bayesian studies on neural networks focus on specifically designed stochastic architectures [START_REF] Dziugaite | Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data[END_REF][START_REF] Zhou | Non-vacuous generalization bounds at the ImageNet scale: a PAC-Bayesian compression approach[END_REF][START_REF] Pérez-Ortiz | Tighter risk certificates for neural networks[END_REF]Clerico et al., 2022a), there are also de-randomisation methods that allow to apply PAC-Bayesian ideas to deterministic classifiers. For instance, one can approximate a deterministic model with a noisy version, whose parameters are stochastically perturbed. Under suitable stability assumptions, a PAC-Bayesian analysis of the stochastic model can bring a generalisation bound for the deterministic one [START_REF] Neyshabur | A PAC-Bayesian approach to spectrally-normalized margin bounds for neural networks[END_REF][START_REF] Nagarajan | Deterministic PAC-Bayesian generalization bounds for deep networks via generalizing noise-resilience[END_REF][START_REF] Miyaguchi | PAC-Bayesian transportation bound[END_REF]. Alternatively, one can average a stochastic algorithm to get a deterministic one. Then, PAC-Bayesian bounds on the loss of this averaged predictor can be obtained by Jensen's inequality [START_REF] Germain | PAC-Bayesian learning of linear classifiers[END_REF][START_REF] Letarte | Dichotomize and generalize: PAC-Bayesian binary activated deep neural networks[END_REF][START_REF] Biggs | Differentiable PAC-Bayes objectives with partially aggregated neural networks[END_REF] (if the loss function is convex), or through margin-based arguments [START_REF] Langford | PAC-Bayes & margins[END_REF][START_REF] Banerjee | De-randomized PAC-Bayes margin bounds: Applications to nonconvex and non-smooth predictors[END_REF][START_REF] Biggs | On margins and derandomisation in PAC-Bayes[END_REF]. Finally, besides the PAC-Bayesian bounds in expectation, there are disintegrated bounds that hold with high probability with respect to a random realisation of the stochastic algorithm. For instance, in the case of a stochastic network, this means that we have a generalisation bound for the deterministic classifier obtained by drawing all the parameters once [START_REF] Catoni | PAC-Bayesian supervised classification: The thermodynamics of statistical learning[END_REF][START_REF] Blanchard | Occam's hammer[END_REF][START_REF] Rivasplata | PAC-Bayes analysis beyond the usual bounds[END_REF][START_REF] Viallard | A general framework for the disintegration of PAC-Bayesian bounds[END_REF].

When trying to apply the PAC-Bayesian framework to study a given non-stochastic model, the first of the methods proposed above is the most natural. However, in the present work, we show that for this purpose it is also possible to leverage the disintegrated PAC-Bayesian bounds. Our starting point is noticing that even training a deterministic classifier does usually involve some randomness: the standard procedure when training a neural network almost always starts with the random draw of the initial weights from some known simple distribution [START_REF] Goodfellow | Deep Learning[END_REF]. For a model trained via continuous-time (nonstochastic) gradient descent, we show that it is possible to exploit this source of stochasticity and obtain a disintegrated PAC-Bayes bound that holds with high probability on the random training dataset and initialisation (Proposition 2). Evaluating this bound only requires to be able to evaluate the density for the initial and final configurations of the network, and to compute the integral along the trajectory of the training objective's Laplacian. Although the bound necessitates the ideal situation of continuous training dynamics, to our knowledge this is the first PAC-Bayesian result that directly apply to a standard non-stochastic neural networks with a fully deterministic training, without strong assumptions on the architecture and the need of the introduction of any auxiliary stochastic model. Besides, we find our bound remarkably easy-to-prove, clean, and simple, with no strong hypotheses made on the training objective, which is only required to be twice differentiable.

The paper is organised as follows. First, we set up the learning framework and introduce the notations in Section 2. Section 3 gives a quick overview on the disintegrated PAC-Bayesian bounds, while in Section 4 we establish our main result (Proposition 2). Section 5 discusses the main features of our bound, and finally Section 6 presents some analytical and empirical quantitative results.

Notation and setting

In a standard supervised learning framework, we consider pairs (x, y) ∈ X × Y, where x denotes an example and y its label. For simplicity, we can imagine that there is a mapping f : X → Y that associates to each x a unique correct label y = f (x). A learning algorithm takes a labelled training dataset s = {x 1 . . . x m }, made of m inputs whose labels {f (x 1 ) . . . f (x m )} are known, and outputs a function f : X → Y, which ideally is a good approximation of f . More generally, the algorithm will return an hypothesis h ∈ H (e.g., values of a neural network's parameters), which is understood to parameterise a function f h : X → Y. For simplicity, we will always assume that H ⊆ R N , for some dimension N > 0. We call the algorithm stochastic when its output h is a random variable on H, whose law can depend on s.

The performance of each hypothesis can be assessed via a loss function , which takes a pair of labels (y, y ) and returns a real value, representing how "far apart" y and y are. In practice, for any hypothesis h, one can always compute the empirical loss on the training dataset

L s (h) = 1 m x∈s (f h (x), f (x)) .
However, we are essentially interested in how well h can predict the label of an x which does not belong to s. Assuming that the population of examples follows a distribution µ, the actual relevant quantity is the population loss

L X (h) = X (f h (x), f (x))dµ(x) .
Upper-bounding L X , having access to L s only, is usually referred to as the generalisation problem. Typically, the training dataset is made of m i.i.d. samples from µ (i.e., s ∼ µ m = µ ⊗m ). We are interested in upper bounds on L X (h) (where h is the hypothesis picked by the algorithm) that holds with high probability on the random draw of s (or of (h, s) in the stochastic setting).

We will consider an algorithm whose output is obtained by optimising an objective C s : H → R via continuous-time gradient descent. Clearly, C s depends on the training dataset s and, in practice, it can coincide with the empirical loss, although this is not necessarily the case, as one might use a surrogate loss for the training or add some regularisation term. We require that the objective is twice differentiable, and we define the gradient flow (h 0 , t) → Φ t (h 0 ) as the solution of the dynamics

∂ t Φ t (h 0 ) = -∇C s (Φ t (h 0 )) ; Φ 0 (h 0 ) = h 0 .
As h 0 is fixed before starting the training, we will often omit the explicit dependence on it, and simply write the solution of the above dynamics as h t . Fixed a time horizon T > 0, our algorithm's output is given by h T .

We will focus on the case of a random initialisation, with law ρ 0 . We define the distribution ρ t as the push-forward of ρ 0 under the gradient flow, that is

ρ t = Φ # t ρ 0 .
In this way, if we draw h 0 from ρ 0 , then we automatically get that h t ∼ ρ t . Note that h t is non-random when conditioned on the initialisation h 0 . In particular, although the algorithm is in fact stochastic, we can see it as a deterministic model that has only a random initialisation. This is actually the common setting for modern neural networks, whose initial configuration is usually drawn from some simple distribution (e.g., [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]).

Disintegrated PAC-Bayes bounds

The main focus in the development of the PAC-Bayesian theory has been on bounds in expectation, and only recently [START_REF] Rivasplata | PAC-Bayes analysis beyond the usual bounds[END_REF] and [START_REF] Viallard | A general framework for the disintegration of PAC-Bayesian bounds[END_REF] have brought back interest in the disintegrated bounds, which actually date back to the works of [START_REF] Catoni | PAC-Bayesian supervised classification: The thermodynamics of statistical learning[END_REF] and [START_REF] Blanchard | Occam's hammer[END_REF]. We refer to [START_REF] Alquier | User-friendly introduction to PAC-Bayes bounds[END_REF] for a detailed introductory exposition of the PAC-Bayesian framework that also discusses some disintegrated results.

We consider a stochastic model, which takes a sample s ∼ µ m and returns a random hypothesis h ∼ ρ s , where the superscript s that we use in this section stresses explicitly that ρ depends on s. 1 We denote the joint law of (s, h) as

µ m * ρ s , that is d(µ m * ρ s )(s, h) = dµ m (s)dρ s (h).
As it is standard in the PAC-Bayesian literature, we will call ρ s the posterior distribution. We denote as π a prior distribution on H, whose only requirement is to be data-agnostic, in the sense that it cannot depend on the specific dataset s used for the training. We write µ m ⊗ π for the law of a pair (s, h), where s ∼ µ m and h ∼ π are independent.

A disintegrated PAC-Bayesian bound is an upper bound on L X (h), where h is the output of a stochastic algorithm, that holds with high probability over (s, h) ∼ µ m * ρ s . The next result, from [START_REF] Rivasplata | PAC-Bayes analysis beyond the usual bounds[END_REF], allows to derive several of these bounds. We provide a short proof for completeness.

Proposition 1. Let Ψ : R 2 → R be an arbitrary measurable function, and define ξ as the expectation of e Ψ(Ls(h),L X (h)) under (s, h) ∼ µ m ⊗ π. For any δ ∈ (0, 1), for any posterior distribution ρ s absolutely continuous with respect to the prior π, with probability at least 1 -δ on (s, h) ∼ µ m * ρ s we have that

Ψ(L s (h), L X (h)) ≤ log dρ s dπ (h) + log ξ δ .
Proof. By Markov's inequality, for (s, h) ∼ µ m * ρ s we have that with probability at least 1 -

δ e Ψ(Ls(h),L X (h))-log dρ s dπ (h) ≤ 1 δ X m ×H e Ψ(Ls(h),L X (h))-log dρ s dπ (h) dµ m (s)dρ s (h) .
By noticing that for all s ∈ X m

H e Ψ(Ls(h),L X (h))-log dρ s dπ (h) dρ s (h) = H e Ψ(Ls(h),L X (h)) dρ s dπ (h) χ dρ s dπ (h)>0 dρ s (h) = H e Ψ(Ls(h),L X (h)) dπ(h) ,
we easily get the desired result.

Choosing Ψ(u, v) = 2m(u -v) 2 one finds the disintegrated version of a standard PAC-Bayes bound from [START_REF] Mcallester | PAC-Bayesian model averaging[END_REF]. Alternatively, defining kl(u v) = u log u v + (1 -u) log 1-u 1-v (the relative entropy between two binary Bernoulli distributions) and setting Ψ(u, v) = mkl(u v), we get a tighter bound, which is the disintegrated counterpart of a bound due to [START_REF] Langford | Bounds for averaging classifiers[END_REF] and [START_REF] Maurer | A note on the PAC Bayesian theorem[END_REF]. These results follow from the fact that in both cases ξ can be upper-bounded by 2 √ m, if the loss is bounded in [0, 1] [START_REF] Bégin | PAC-Bayesian bounds based on the Rényi divergence[END_REF].

Corollary 1. Assume that the loss is bounded in [0, 1]. Then, for any δ ∈ (0, 1), for any posterior distribution ρ s absolutely continuous with respect to the prior π, each of the following bounds hold with probability at least

1 -δ on (s, h) ∼ µ m * ρ s L X (h) ≤ L s (h) + log dρ s dπ (h) + log 2 √ m δ 2m ; L X (h) ≤ kl -1 L s (h) log dρ s dπ (h) + log 2 √ m δ m ,
where kl -1 (u|c) = sup{v ∈ [0, 1] : kl(u v) ≤ c}.

Disintegrated bounds under gradient flow dynamics

We state here our main result, a disintegrated bound for an algorithm trained by continuous gradientdescent optimisation. We recall some notation from Section 2. We denote as C s a twice differentiable training objective, and we consider a random initial state h 0 ∼ ρ 0 . The dynamics are described by

∂ t h t = -∇C s (h t ) ,
where we omit the explicit dependence of h t on s and h 0 . The algorithm's output is given by h T , with T > 0 some fixed time horizon. We define ρ t as the push-forward of ρ 0 under the gradient flow, so that h t ∼ ρ t . For simplicity, for all t we assume that ρ t admits a density with respect to the Lebesgue measure, which we will again denote as ρ t , with a slight abuse of notation.

It is natural to suppose that ρ 0 is chosen without relying on the specific dataset s used for the training, so that, in particular, we can select ρ 0 as our PAC-Bayesian prior. On the other hand, ρ T depends on s through the objective C s , and plays the role of the posterior. Sampling the initialisation h 0 of the model from ρ 0 , and then following the deterministic training dynamics up to T , we get a final hypothesis h T , which can be seen as a sample from ρ T . In particular, saying that an event holds with probability at least 1 -δ under (s, h T ) ∼ µ m * ρ T is equivalent to state that it has probability at least 1 -δ under (s, h 0 ) ∼ µ m ⊗ ρ 0 . We can thus use Proposition 1 to obtain a PAC-Bayesian generalisation bound for an algorithm that, once drawn s and h 0 , is deterministic.

So far, we have a bound that involves the term ρ T (h T )/ρ 0 (h T ). Typically, the density ρ 0 is known, and evaluating ρ 0 (h T ) does not pose any problem. On the other hand, computing ρ T (h T ) is less straightforward. However, thanks to the continuity of the gradient descent dynamics, it is possible to keep track of the evolution of ρ t (h t ) along the training trajectory via a simple trick, usually referred to as "instantaneous change of variables" in the normalising flow literature [START_REF] Chen | Neural ordinary differential equations[END_REF].

Proposition 2. Consider the dynamics ∂ t h t = -∇C s (h t ), where C s : H → R is twice differentiable. Let Ψ : R2 → R be a measurable function. Fixed T > 0 2 and δ ∈ (0, 1), with probability at least 1 -δ on the random draw (s, h 0 ) ∼ µ m ⊗ ρ 0 , we have

Ψ(L s (h T ), L X (h T )) ≤ log ρ 0 (h 0 ) ρ 0 (h T ) + T 0 ∆C s (h t )dt + log ξ δ ,
where ∆ denotes the Laplacian with respect to h and ξ = X m ×H e Ψ(Ls(h),L X (h)) dµ m (s)dρ 0 (h).

Proof. Since sampling (s, h T ) from µ m * ρ T is equivalent to drawing (s, h 0 ) ∼ µ m * ρ 0 and following the dynamics up to T , Proposition 1 yields directly that, with probability at least 1 -δ on (s,

h 0 ) ∼ µ m ⊗ ρ 0 , Ψ(L s (h T ), L X (h T )) ≤ log ρ T (h T ) ρ 0 (h T ) + log ξ δ .
The gradient flow's continuity equation states that, for all h ∈ H,

∂ t ρ t (h) = ∇ • (ρ t (h)∇C s (h)) = ∇ρ t (h) • ∇C s (h) + ρ t (h)∆C s (h) .
It follows that

∂ t (ρ t (h t )) = ∂ t ρ t (h t ) + ∇ρ t (h t ) • ∂ t h t = ρ t (h t )∆C s (h t ) , which yields log ρ T (h T ) ρ 0 (h 0 ) = T 0 ∆C s (h t )dt . Hence log ρ T (h T ) ρ 0 (h T ) = log ρ 0 (h 0 ) ρ 0 (h T ) + T 0 ∆C s (h t )dt ,
and we conclude.

Note that the time horizon T must be chosen a priori, as it cannot depend on the specific s and h 0 used for the training. However, through a classical union bound argument, one can slightly loosen the bound while allowing to select the optimal T from a finite pool of possible values. Concretely, it is possible to pick the best T from a set containing K potential choices at the cost of replacing the term log ξ δ with log Kξ δ , in the bounds of Propositions 1 and 2. Clearly, for ⊆ [0, 1], the choices of Ψ that lead to Corollary 1 will again bring explicit bounds. Setting

Ψ(u, v) = 2m(u -v) 2 yields L X (h T ) ≤ L s (h T ) + log ρ0(h0) ρ0(h T ) + log 2 √ m δ + T 0 ∆C s (h t )dt 2m , (1)
while from Ψ(u, v) = mkl(u|v) we get

L X (h T ) ≤ kl -1   L s (h T ) log ρ0(h0) ρ0(h T ) + log 2 √ m δ + T 0 ∆C s (h t )dt m   , (2)
where each of the above two bounds holds with probability higher than 1 -δ on the random initialisation and training dataset.

We remark that Proposition 2 holds also for time-dependent objectives (as long as they do not depend on ρ t ), as the gradient flow's continuity equation is still valid. Leveraging this, we can specialise the result for training with (potentially randomly selected) batches. We refer to Section 5.5 for further details.

Finally, an interesting alternative scenario is when the continuous-time dynamics are stochastic. For instance, one might want to look at the Langevin's diffusion dh t = -∇C s (h t ) + σdB t , where σ > 0 is a fixed coefficient and B t a Brownian motion. However, in this case we do not end up anymore with an easily tractable exact expression for log ρ T (h T ) ρ0(h0) (see Appendix A).

Analysis and discussion

A few comments on the Laplacian term

An interesting feature of the bound in Proposition 2 is the integral of the Laplacian of the optimisation objective along the training path:

T 0 ∆C s (h t )dt . (3)
Under the gradient flow dynamics, ∆C s (h t ) is keeping track of how the local probability density is locally varying. Indeed, from a PAC-Bayesian point of view, for the bound to be small we need to end up in some point where the posterior density is not too high compared to the initial prior. If we are following a trajectory characterised by high values of the Laplacian, we will see a sharp increase of the density. Intuitively, we can picture the situation as if we were attracting the nearby paths and bringing further probability mass around us. In this case we are likely to end with a final ρ T that is much larger than the initial ρ 0 , and hence, potentially, a loose bound.

It is worth noticing that, in the one-dimensional case, (3) takes a particularly simple form, as it only depends on the initial and final "velocities". To see this, note that

∂ t C s (h t ) = C s (h t )∂ t h t = -C s (h t )C s (h t ),
with C s and C s denoting the first and second derivatives of

C s . It follows that C s (h t ) = -∂ t log |C s (h t )|, which yields T 0 C s (h t )dt = log C s (h 0 ) C s (h T ) .
On the other hand, things are not that simple in the multidimensional case. With a few algebraic manipulations, we can write

T 0 ∆C s (h t )dt = log ∇C s (h 0 ) ∇C s (h T ) - h [0:T ] ∇ • τ (h) δh , (4)
where τ (h) stands for the unit tangent vector to the trajectory in h, and h is attracting the nearby trajectories. We refer to Appendix B for the derivation of (4).

Although we are not aware of any generalisation bound that involves the exact term (3), the Laplacian's integral might find a natural interpretation through the connection between flatness in the loss landscape and good generalisation. It dates back to [START_REF] Hochreiter | Flat minima[END_REF] the belief that an algorithm, able to find wide local optima of the loss, performs well on unseen inputs. Re-popularised by [START_REF] Keskar | On large-batch training for deep learning: Generalization gap and sharp minima[END_REF], this idea has been the object of an intense debate, attracting some criticisms as the definitions of flatness are often susceptible to re-parameterisations of the hypothesis space [START_REF] Clerico | Chained generalisation bounds. COLT[END_REF][START_REF] Neyshabur | Exploring generalization in deep learning[END_REF]. However, numerous recent works have provided empirical and theoretical support, as well as new insights, to the flat minima's argument [START_REF] Izmailov | Averaging weights leads to wider optima and better generalization[END_REF][START_REF] He | Asymmetric valleys: Beyond sharp and flat local minima[END_REF][START_REF] Dziugaite | Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data[END_REF][START_REF] Neu | Information-theoretic generalization bounds for stochastic gradient descent[END_REF]. Going back to our setting, if we assume that the objective C s coincides with the empirical loss, we clearly see that sharper loss's local optima correspond to higher values of the Laplacian term. As a consequence, converging towards flatter minima will typically yield a tighter bound.

Role of the time horizon

Most of the standard PAC-Bayesian bounds (in expectation) contain the Kullback-Leibler divergence between the prior and the posterior. This quantity is infinity if the posterior is a Dirac delta, and this is why these results cannot apply directly to the deterministic setting. Even the disintegrated framework is not exempt from this pathological behaviour, which is reflected by the fact that in Proposition 2 we cannot let the time horizon T to be infinite even if h t converges, as the posterior would be degenerate and typically consist in a sum of Dirac deltas centred on the local optima. Concretely, we can easily see that if an algorithm is converging to a local minimum (where C s is approximately quadratic), then (1) will asymptotically go as

L X (h) -L s (h) O( √ T )
for large T (see Appendix C for more details). However, this is not the only possible asymptotic behaviour.

For instance, we can consider situations where h t is not converging but diverging. As we will see in Section 6, this is what happens for linear models trained to optimise a linear objective or the cross-entropy.

Although characterising the typical spectrum of the Hessian of the loss during the training evolution is a challenging open problem, recent empirical studies suggest that after an initial phase where the Laplacian can potentially be negative, later stages of the dynamics are characterised by a few dominant positive Hessian's eigenvalues [START_REF] Sagun | Empirical analysis of the Hessian of overparametrized neural networks[END_REF][START_REF] Ghorbani | An investigation into neural net optimization via Hessian eigenvalue density[END_REF]. As a consequence, in most situations the Laplacian's integral will soon become monotonically increasing with T . Moreover, we can expect that log ρ0(h0) ρ0(h T ) will grow with T as well, since h 0 is likely to lie in a region where the prior density is high, and the training will typically push h t farther and farther from it. Overall, in order to achieve non-vacuous bounds, we need to stop the training at a finite time horizon T , which will play the role of a trade-off between how well we can fit the training dataset, and how much we can keep the bound small. A final remark, already mentioned earlier, is that T must be chosen independently of the dataset s and initial state h 0 used for the training. However, adding a penalty log K to the bounds of Propositions 1 and 2 allows us to pick the best time horizon among a set {T 1 . . . T K } of K possible choices. This is the consequence of an elementary union argument, as for each T k we can consider a bound holding with probability at least 1 -δ/K.

Scaling with the dimension of the hypothesis space

When dealing with generalisation bounds, a natural question is how they scale with the dimension N of the hypothesis space. In our case, the answer is not easy, as it is hard to say how the different components of bound depend on N . For instance, we can expect the Laplacian to scale linearly with N , as it is the sum of N terms. If T = O(1) and the Laplacian's integral dominates the bound, (1) becomes O( N/m), which is a standard behaviour. On the other hand, if the log-density term log ρ0(h0) ρ0(h T ) is the main contribution to the bound, the scaling with N might be more complex, as it will depend on the particular choice of the prior density.

Interestingly, the optimal time horizon T can scale with N as well, and empirical results on MNIST show that for a linear model we can have T ∼ O(1/N ), and hence obtain non-vacuous bounds even in a highly over-parameterised setting (see Figure 1). Intuitively, when N is large, the algorithm can quickly find an hypothesis that fit well the training data and is "close" to h 0 , hence requiring a low training time to be reached. On the other hand, fewer learnable parameters might need a longer training time before reaching low values of the empirical loss. Despite the simple setting of our experiment, we believe that this is a promising result as the bound seems unaffected by the curse of dimensionality. However, a more detailed empirical and theoretical analysis is needed in order to gauge how well our bound performs with more complex models, such as deep networks. (2) from exploding, even in highly over-parameterised settings. Note that the results reported here were obtained via discretised gradient dynamics (cf. Section 5.4 and Figure 2). We refer to Appendix E.3 for further details on the experiment. (2) obtained with the time-steps δt,1 = 2 × 10 -4 and δt,2 = 2 × 10 -6 respectively. The small relative error suggests that the corrections to the bound due to the time-discretised dynamics are likely to be negligible. We refer to Appendix E.3 for further details on the experiment.

Discretised dynamics

One of the main limitations of the present work is that the bound of Proposition 2 requires a continuoustime gradient descent. Apart from a few simple settings, where analytic solutions are available, the standard gradient-based optimisation algorithms consider discretised dynamics, whose trajectories are made of jumps and only approximate the continuous flow. As a consequence, in order to leverage our results and get rigorous generalisation bounds for the standard gradient-descent optimisers, one would need to control the errors introduced by the time discretisation via a careful analysis, which is beyond the scope of the present work. However, we can at least conjecture that, for small enough learning rates, the results obtained via the discretised dynamics are a good approximation of the exact solution of the gradient flow. In Figure 2, the relative error between the values of the bound (2), obtained from two different time discretisations, stays small throughout the training, suggesting that the discretised results can be good approximations of the actual continuous dynamics.

Interestingly, [START_REF] Barrett | Implicit gradient regularization[END_REF] remarked that the solutions of the discretised dynamics h k+1h k = -ε∇C s (h k ) follow closely the continuous paths yielded by ∂ t h t = -∇(C s (h t ) + ε ∇C s (h t ) 2 ), with an error after a single step of order O(ε3 ) 3 . Consequently, if we aim at studying the discrete dynamics, we could obtain a more precise approximation of an actual generalisation bound by replacing the term ∆C s with ∆(C s (h t ) + ε ∇C s (h t ) 2 ). If higher accuracy is required, one can add further correction terms to C s , as it is standard in the backward error analysis of Runge-Kutta methods [START_REF] Hairer | Geometric numerical integration[END_REF].

Training via batches

So far, we have considered the case of a training objective C s which is fixed throughout the training. However, the bound of Proposition 2 holds also for time-dependent objectives. In particular, we can apply our result to the case in which the training is performed by optimising the empirical loss over different batches of the training dataset. Let us denote as s 1 . . . s M the (temporarily ordered) sequence of batches used for the training, where a batch needs to appears multiple times if it is used more than once. It is worth noticing that the batches and their order are allowed to be selected randomly, as it is usually the standard stochastic gradient descent optimisation. Then, we can suppose that from t = 0 to t = t 1 we use C s1 as optimisation objective, from t = t 1 to t = t 2 , we use C s2 , and so on until t M = T . We can get a bound as the one of Proposition 2, with the only difference that T 0 ∆C s (h t )dt will now be replaced by

t1 0 ∆C s1 (h t )dt + t2 t1 ∆C s2 (h t )dt + • • • + t M t M -1 ∆C s M (h t )dt .
(5)

Data-dependent priors

For most of PAC-Bayesian results, the prior is required to be chosen independently of the dataset used to evaluate the bound. However, we can use an additional dataset s (independent of s) to learn the prior. This strategy was first proposed by [START_REF] Seeger | PAC-Bayesian generalisation error bounds for Gaussian process classification[END_REF] and has lead to the tightest currently available empirical PAC-Bayesian bounds on MNIST and CIFAR10 (Pérez-Ortiz et al., 2021;Clerico et al., 2022a).

In our setting we might proceed in a similar way. Till a fixed time t 0 , we follow the dynamics ∂ t h t = -∇C s (h t ), and we select ρ t0 as the PAC-Bayesian prior. Then, we take C s as training objective and we go on as usual, starting from h t0 and following the gradient flow until the time horizon T . In order to evaluate the bound, we need to compute log ρ T (h T ) ρt 0 (h T ) . Proceeding as in the proof of Proposition 2, we can get

log ρ T (h T ) ρ t0 (h T ) = log ρ t0 (h t0 ) ρ t0 (h T ) + T t0 ∆C s (h t )dt = log ρ 0 (h t0 ) ρ t0 (h T ) + t0 0 ∆C s (h t )dt + T t0 ∆C s (h t )dt .
Now, the only problem is to evaluate log ρ t0 (t T ), as we do not have a closed form for ρ t0 . However, it is possible to exactly compute this term as

log ρ t0 (h T ) = log ρ 0 ( ĥT -t0 ) + t0 0 ∆C s ( ĥT -t0+t )dt , ( 6
)
where ĥt is the solution of ∂ t ĥt = -∇C s ( ĥt ) satisfying ĥT = h T . We hence get

log ρ T (h T ) ρ t0 (h T ) = log ρ 0 (h t0 ) ρ 0 ( ĥT -t0 ) + t0 0 (∆C s (h t ) -∆C s ( ĥT -t0+t ))dt + T t0 ∆C s (h t )dt .
Plugging this expression in Proposition 1 we readily obtain a bound whose prior ρ t0 is data-dependent. We refer to Appendix D for an explicit statement of the bound and a derivation of (6).

Linear models

To give an idea of the behaviour of the bounds that we propose, we report here some analytic and empirical results. We consider models linear in the parameters, in the form

F i h (x) = j h ij Φ j (x) ,
where the Φ j are some fixed mappings R → R. For binary classification tasks, where y ∈ {-1, +1}, we let F h be a real function and the model's prediction f h corresponds to the sign of F h . Conversely, for general multi-class problems the predicted label f h is the index of the largest component of F h .

As it is standard in the classification setting, we can assume that is the 01-loss4 , so that L s and L X represent the empirical and population error. Unless otherwise stated, we will always consider an optimisation objective in the form

C s (h) = Ls (h) = 1 m x∈s ˆ (F h (x), f (x)
), with ˆ denoting some twice differentiable surrogate loss function. Applying the chain rule to ˆ , we find that

∆ ˆ = i ∂ ˆ ∂F i ∆F i + ii ∂ 2 ˆ ∂F i ∂F i ∇F i • ∇F i ,
where ∆ and ∇ always refer to derivatives with respect to h. In the linear models that we are considering ∆F = 0 and

∇F i • ∇F i = δ ii Φ 2 . Thus ∆ ˆ = Φ 2 ∆ F ˆ , (7)
with ∆ F ˆ = i ∂ 2 F i ˆ , the Laplacian of ˆ with respect to F .

Linear loss

A first simple example that we can solve analytically is a binary classification problem with ˆ (F, y) = -F y .

The surrogate loss is linear in F , and so (7) yields ∆ ˆ = 0. In particular, the Laplacian's integral in our bounds vanishes. We can easily solve the evolution dynamics and find that

h T = h 0 + γT ,
where, recalling that f (x) is the correct label of the example x, we have defined γ = 1 m x∈s f (x)Φ(x). Assuming that under ρ 0 the components of h are all independent centred random variables, with variance 1/N (N being the dimension of h), we find that log ρ0(h0) ρ0(h

T ) = N 2 ( h T 2 -h 0 2 ). So, (1) becomes L X (h T ) -L s (h T ) ≤ N 4m (2T h 0 • γ + T 2 γ 2 ) + 1 2m log 2 √ m δ .
Tipycally, h 0 ∼ O(1), and we expect h T ∼ T γ to dominate the dynamics. Notice that the bound scales linearly with T . This is different from the O( √ T ) behaviour that we had encountered in section 5.2. Indeed, here h is not approaching a local optimum, but diverging at constant speed γ.

We tested empirically the behaviour of the bounds (1) and ( 2) with a linear surrogate loss, by training a model where the learnable vector h had 10000 components. The dataset was made of 500 datapoints, coming from Gaussian clusters in R5 and labelled into two categories. The map Φ was in the form x → φ(W x), with φ the ReLU function and W a (5 × 10000)-matrix whose entries were independent draws from a centred Gaussian with variance 1/5. We refer to Appendix E.1 for further experimental details and results. Figure 3 compares the behaviour of the two bounds with the test errors obtained from a held-out dataset, for 50 different possible values of the time horizon 5 and for δ = 5 × 10 -3 . Note that for large T (1) is linear, as expected. 2) for a linear model with 10000 parameters, trained with a linear surrogate loss, for different time horizons T . The data came from eight different Gaussian distributions on R 5 that had been assigned binary labels. The training dataset counted 500 datapoints and a held-out test dataset was used to evaluate the test errors reported in the plot. We refer to Appendix E.1 for further details. 2) for a linear model with 10000 parameters, trained with a quadratic surrogate loss, for different time horizons T . The learning task was the same as in Figure 3. Again, we used a training dataset made of 500 datapoints and we evaluated the test errors on a held-out test dataset. The parameters α and β in Cs were both set to 1. We refer to Appendix E.2 for further details.

Quadratic loss

Another natural situation where we can solve things explicitly is the case of a quadratic surrogate loss

ˆ (F, y) = 1 2 (F -βy) 2 .
Here, β is some fixed positive value, and we are still considering a binary classification problem, where the true labels are +1 or -1 and the network's prediction is the sign of F . Now, ∆ F ˆ = 1 and hence (7) yields ∆ ˆ = Φ 2 . We consider a training objective in the form

C s (h) = α h 2 2N + 1 m x∈s ˆ (F h (x), f (x)) ,
with α ≥ 0 and N the dimension of h. The regularising term α h 2 /(2N ) is here introduced to computationally stabilise the algorithm. 6 We easily find that

T 0 ∆C s (h t )dt = ΓT , with Γ = α + 1 m x∈s Φ(x) 2 .
The training dynamics are now governed by ∂ t h t = βγ -Θh t , where again we have defined γ

= 1 m x∈s f (x)Φ(x), while Θ = α N Id + 1 m x∈s Φ(x)Φ(x) is a N × N matrix. Hence, assuming that Θ is invertible, we have h T = h 0 + (Id -e -T Θ )(Θ -1 βγ -h 0 ) .

If again we consider a Gaussian initialisation, (1) becomes

L X (h T ) -L s (h T ) ≤ N 2 ( h 0 2 -h T 2 ) + ΓT 2m .
As the objective is quadratic, for large T this bound will behave as O( √ T ), as it happens more generally every time the algorithm is converging to a local minimum of the objective (see Section 5.2).

Figure 4 shows the bounds obtained by optimising a quadratic objective with α = β = 1. The dataset and the setting were those that we used for the linear loss's training, which we described in the Section 6.1. Again, we used δ = 5 × 10 -3 . We refer to Appendix E.2 for further experimental details and results.

Cross-entropy loss

For a general classification task, which can possibly have more than two labels, when f h = argmax i F i h a common optimisation objective is the cross-entropy loss

ˆ (F, y) = log i e F i -F y .
In this case we do not have simple analytic dynamics, and we will need to rely on some numerical solution.

We have

∂ t h ij t = - 1 m x∈s e (htΦ(x)) i k e (htΦ(x)) k -δ if (x) Φ j (x) ,
where δ if (x) denotes the Kronecker delta between the labels i and f (x), and

∆ ˆ = Φ 2 1 - k e 2F k ( k e F k ) 2 .
Contrarily to what happened with the quadratic and linear losses, where the objective's Laplacian stayed constant during the training, here we are in the interesting situation where the Laplacian tends to decrease during the training. Indeed, for ˆ (F, y) to be small, we need F y to be much larger than the other components of F , so that k e F k e F y . But this also means that the Laplacian is small, since

( k e 2F k )/( k e F k ) 2
1. This is a different behaviour from what happens when the algorithm is converging to a local minimum of the objective, as as now we can expect that the Lapacian's integral will have a sub-linear behaviour in T . Indeed, in the present situation, we are not approaching a local optimum. Conversely, h t is slowly diverging, reaching flatter and flatter regions of H.

We trained linear models of different sizes on MNIST, using the cross-entropy loss as objective. Since the dynamics does not admit a simple analytic solution, we performed the classical gradient-descent algorithm on batches of the training dataset, with finite time-steps. As discussed in Section 5.4, what we obtain in this way are not rigorous generalisation bounds, as we are not following the continuous dynamics. However, we believe that the result that we get are a good approximation of the ones that one would obtain with the exact solution of the gradient flow (see Figure 2 and the discussion in Section 5.4). The feature mapping Φ that we used was again the composition of a ReLU function and a linear transformation x → W x, where the components of W were independently sampled from a normal distribution, with variance 1/784 (784 being the dimension of the input). We refer to Appendix E.3 for further experimental details.

Table 1 compares the bound (2) with the test error evaluated on a held-out dataset, for different dimensions of the hypothesis space H, and always with δ = 5 × 10 -3 . We chose the best bound among those obtained with 25 different values of T , at the price of adding a penalty factor log 25 (see Section 5.2).

As we had already stressed in Figure 1 (which refers to the same experiment), note that the bounds are not negatively affected by the high over-parameterisation, essentially because the optimal time horizon T decreases roughly as 1/N . The tightest bound was achieved for a network of width 10000 (i.e., N = 10 5 ). For this model, Figure 5 shows the evolution of the bounds (1) and ( 2) during the training. The values of the empirical error and of the penalty factor (log ρ0(h0) ρ0(h T ) + log 25 √ m δ + T 0 ∆C s (h t )dt)/m are also reported. Figures 6 compares the contributions of the Laplacian's integral and the log-density term, showing that in this case it is the latter that dominates.

As a final remark, we notice that the bounds of Table 1 are comparable to those reported by [START_REF] Clerico | Wide stochastic networks: Gaussian limit and PAC-Bayesian training[END_REF] for a stochastic single-hidden-layer network. Notably, the results therein were obtained by using a PAC-Bayesian bound as optimisation objective, while here we are directly optimising the empirical cross-entropy loss, without any regularisation term.

Conclusion

We established disintegrated PAC-Bayesian bounds for deterministic algorithms, holding with high probability on the random training dataset and initialisation. We tested our findings empirically on linear models, obtaining promising results that are non-vacuous even for highly over-parameterised settings. It is worth noticing that most of empirical PAC-Bayesian results in the literature are derived by using a PAC-Bayesian bound as the optimisation objective [START_REF] Dziugaite | Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data[END_REF][START_REF] Pérez-Ortiz | Tighter risk certificates for neural networks[END_REF][START_REF] Clerico | Wide stochastic networks: Gaussian limit and PAC-Bayesian training[END_REF]Clerico et al., , 2022a)), while our bounds directly apply to the standard (continuous-time) empirical risk optimisation. Yet, our results on MNIST closely match the bounds obtained by [START_REF] Clerico | Wide stochastic networks: Gaussian limit and PAC-Bayesian training[END_REF]. We defer to future work a more extensive empirical analysis, focusing on data-dependent priors and non-linear models.

A main limitation of our framework is the necessity of continuous-time dynamics. However, it might be possible to obtain rigorous results for the standard discretised gradient-descent optimisation via a careful numerical analysis. Alternatively, a scenario which might be worth exploring is the limit of infinite width, where the NTK dynamics [START_REF] Jacot | Neural tangent kernel: Convergence and generalization in neural networks[END_REF]) might provide interesting continuous-time exact results. Finally, we mention as possible further developments a deeper study of stochastic training dynamics, where the gradient flow's continuity relation is replaced by the Fokker-Planck equation, and the case of a training objective that depends on the bound itself.

Appendix A Langevin dynamics

We briefly discuss what happens if we try to use a Langevin dynamics instead of a deterministic gradient method, that is dh t = -∇C s (h t )dt + σdB t .

Here, σ > 0 is a diffusion coefficient and B t denotes the standard Brownian motion. Note that the introduction of noise prevents the bound from diverging at large T . Indeed, under suitable regularity conditions for C s , ρ t will converge to the Gibbs equilibrium distribution

ρ ∞ (h) = 1 Z e -2Cs(h)/σ 2 ,
with Z the normalising constant. As ρ ∞ is bounded by 1/Z, asymptotically log ρ T (h T ) ρ0(h T )

-log(Zρ 0 (h T )) does not diverge with T .

For ( 9), the Fokker-Planck equation yields

∂ t ρ t (h) = ∇(ρ t (h)∇C s (h)) + σ 2 2 ∆ρ t (h) = ∇ρ t (h) • ∇C s (h) + ρ t (h)∆C s (h) + σ 2 2 ∆ρ t (h) . (9)
Applying Itō's lemma we get

d(ρ t (h t )) = (ρ t (h t )∆C s (h t ) + σ 2 ∆ρ t (h t ))dt + σ∇ρ t (h t ) • dB t .
Now, writing u t = log ρ t , we have

du t (h t ) = ∆C s (h t ) + σ 2 2 ∇u t (h t ) 2 + σ 2 ∆u t (h t ) dt + σ∇u t (h t ) • dB t , (10)
where we used ∆ρt ρt = ∇u t 2 + ∆u t . This expression cannot easily be integrated in t, as it happened for the deterministic dynamics, since now we have a dependence on u t (and hence on ρ t ) in the RHS. In order to obtain a computable bound, one would probably need to upper-bound the intractable quantities in the RHS, an approach that we defer to future work.

An alternative way to tackle the stochastic dynamics consists in converting the SDE (8) in an ODE inducing the same marginal distributions (an idea already exploited in [START_REF] Song | Score-based generative modeling through stochastic differential equations[END_REF]). In fact, consider the deterministic flow

∂ t ĥt = -∇ C s ( ĥt ) + σ 2 2 log ρt ( ĥt ) , (11)
where we denote as ρt the density of ĥt and we impose ρ0 = ρ 0 . The continuity equation now reads

∂ t ρt (h) = ∇ • ρt (h) ∇ C s (h) + σ 2 2 log ρt (h) = ∇(ρ t (h)∇C s (h)) + σ 2 2 ∆ρ t (h) . (12)
Comparing ( 12) with ( 9), we immediately see that ρ t and ρt undergo the same dynamics and, since ρ 0 = ρ0 , we have that ρ t = ρt for all t. In practice, this means that we have two equivalent ways to get h T ∼ ρ T : given an initial point from ρ 0 , we can follow either the stochastic dynamics (8), or the deterministic ones (11), up to time T . Moreover, in order to evaluate the factor log ρ T (h T ), needed for our PAC-Bayesian bound, instead of using (10) we can leverage the continuity equation of the deterministic flow to derive

log ρ T (h T ) = log ρ 0 ( ĥ0 ) + T 0 ∆ C s ( ĥt ) + σ 2 2 log ρ t ( ĥt ) dt ,
where ĥt denotes the solution of ( 11) satisfying ĥT = h T . Although this expression looks easier than (10), its actual computation is typically unfeasible in practice, due to the RHS's dependence on log ρ t (both explicit and implicit, through ĥt ). In the generative modelling context, score matching and neural networks techniques are used to approximate ∇ log ρ t , but this introduces approximation errors difficult to quantify [START_REF] Song | Score-based generative modeling through stochastic differential equations[END_REF].

As a final remark, we mention that the dynamics (11) also appear when trying to build an algorithm that directly optimises a PAC-Bayesian bound. Indeed, for ⊆ [0, 1], the choice Ψ(u, v) = 2 σ 2 (u -v) (for an arbitrary σ > 0) in Proposition 1 leads to a bound in the form

L X (h T ) ≤ L s (h T ) + 1 4mσ 2 + σ 2 2 log ρ T (h T ) ρ 0 (h T ) + log 1 δ .
Using this bound as optimisation objective is equivalent to follow the dynamics (11), provided that we set C s (h) = L s (h) + 1 4mσ 2 + σ 2 2 (-log ρ 0 (h) + log 1 δ ). We defer to future work a deeper analysis of this connection.

B Rewriting the Laplacian's integral

We explicitly derive here (4). First, notice that

∂ t log ∇C s (h t ) = -∇Cs(ht) ∇Cs(ht) • ∇ ∇C s (h t ) . Since, for all h, ∆C s (h) = ∇ • ∇C s (h) = ∇ • ∇C s (h) ∇C s (h) ∇C s (h) + ∇C s (h) ∇C s (h) • ∇ ∇C s (h) , we get ∆C s (h t ) = ∇ • ∇C s (h t ) ∇C s (h t ) ∇C s (h t ) -∂ t log ∇C s (h t ) .
We conclude that

T 0 ∆C s (h t )dt = log ∇C s (h 0 ) ∇C s (h T ) + T 0 ∇ • ∇C s (h t ) ∇C s (h t ) ∇C s (h t ) dt .
The integral in the RHS is a line-integral along the path h [0,T ] , as ∇C s (h) is the norm of the flow's "velocity" in h. Moreover, τ (h) = -∇Cs(h) ∇Cs(h) is the unit tangent vector to the gradient flow in h. We can thus write

T 0 ∆C s (h t )dt = log ∇C s (h 0 ) ∇C s (h T ) - h [0:T ] ∇ • τ (h) δh ,
which is (4).

C Asymptotic behaviour around a local minimum

We consider the case in which the algorithm is approaching a quadratic local optimum of C s . We will hence say that for t larger than some t we have

C s (h t ) C s (h ) + 1 2 (h t -h ) H(h t -h ) ,
where H is a positive definite matrix, namely the Hessian of C s evaluated at h . Hence, for t ≥ t the objective's Laplacian is roughly constant,

∆C s (h t ) Tr[H] .
For a large enough time horizon T , we can asymptotically write

T 0 ∆C s (h t ) T Tr[H]. Moreover, since the Laplacian's integral is diverging while ρ0(h0) ρ0(h T ) → ρ0(h ) ρ0(h T )
, asymptotically we find that (1) takes the form

L X (w T ) -L s (w T ) O T Tr[H] 2m .

D Data-dependent priors

As we discussed in Section 5.6, if we want to use a prior ρ t0 , learnt on a dataset s (independent of s), we need to evaluate log ρ t0 (h T ). Define ĥt as the solution of ∂ t ĥt = -∇C s ( ĥt ) satisfying ĥT = h T . We have that

log ρ t0 (h T ) = log ρ t0 ( ĥT ) = log ρ 0 ( ĥT -t0 ) + T 0 ∂ t (log ρ t ( ĥT -t0+t ))dt .
Since, for t ∈ [0, t 0 ], ρ t is the push-forward of ρ 0 under the gradient dynamics with objective C s , we have that for all h

∂ t ρ t (h) = ρ t (h)∆C s (h) + ∇ρ t (h) • ∇C s (h)
and so in particular

∂ t (ρ t ( ĥt )) = ∂ t ρ t ( ĥt ) + ∇C s ( ĥt ) • ∂ t ĥt = ρ t ( ĥt )∆C s ( ĥt ) .
We conclude that log ρ t0 (h T ) = log ρ 0 ( ĥT -t0 ) +

) ∼ µ m ⊗ ρ 0 , we have Ψ(L s (h T ), L X (h T )) ≤ log ρ 0 (h t0 ) ρ 0 ( ĥT -t0 ) + t0 0 (∆C s (h t ) -∆C s ( tT -t0+t ))dt + T t0 ∆C s (h t )dt + log ξ δ ,
where ĥt is the solution of ∂ t ĥt = -∇C s ( ĥt ) with ĥT = h T , and ξ = X m ×H e Ψ(Ls(h),L X (h)) dµ m (s)dρ t0 (h).

E Experimental details

E.1 Linear loss

For our experiments with the linear surrogate loss (see Section 6.1), we used a binary toy dataset, which was created as follows. First, we sampled eight independent points in R 5 , from a standard normal distribution. These eight points were then used as the means of eight Gaussian distributions, symmetric and with variance 0.1. We sampled eight Gaussian clusters of 5000 points from each of this eight distributions.

We then assigned the label -1 to four clusters, and +1 to the remaining four. Finally, all the points were projected on the unit sphere. 500 of them, randomly chosen, were used as training dataset, while the remaining 39500 constituted a held-out test dataset.

The linear model that we considered was a neural network with a single-hidden layer, where the hidden nodes were kept fixed at their initial values and only the outer layer was trained. Explicitly, we had

F (x) = h • Φ(x), with Φ(x) = φ(W x).
Here, W denotes a random 5 × N matrix, whose components were independently drawn from a centred normal distribution with variance 1/5, and φ is the ReLU function φ(x) = max(0, x). We chose the distribution ρ 0 as a multivariate normal distribution, so that all the components of h 0 were independent draws from a Gaussian random variable with variance 1/N .

We tried different values for the dimension N of h: 100, 1000, and 10000. We reported in Figure 3 in the main text the results obtained for N = 10000. The parameter δ was fixed at .005, so the bounds hold with probability higher than .995. Since the values reported are already comprehensive of the factor log 50 due to the union bound, it is possible to select the horizon time T yielding the tightest result. Figure 7 compares the values of (2) for the different sizes N .

E.2 Quadratic loss

For the experiments with the quadratic surrogate loss (Section 6.2), we used the same binary dataset and models that we described in Appendix E.1. The parameters α and β that define C s were both set to 1. Figure 4 in the main text reports the results obtained for N = 10000. The value reported are already taking into account the log 50 factor (due to the union bound), and δ was again chosen to be .005.

E.3 Cross-entropy loss

We performed our experiments with the cross-entropy as surrogate loss (see Section 6.3) on the MNIST dataset, with the standard split of 60000 images used for the training and 10000 for the testing of the model. The training was performed on batches of 500 elements. The images were flattened to onedimensional arrays and renormalised. We used again linear architectures F (x) = hΦ(x), but now with h being a 10 × w matrix (w denoting the width of the network). Φ was always in the form Φ(x) = φ(W x).

Here φ acts component-wise as the ReLU function, and W is a w ×784 matrix (784 being the dimension of the input datapoints), whose components were independently drawn from a centred normal distribution with variance 1/784. The dimension of the hypothesis h was N = 10w, as the networks had 10 output nodes. We tested different values of N , as reported in Table 1. In all the bounds δ = .005.

Since closed-form solutions for the flow are not available, we used discretised gradient dynamics. For each batch s k we made an update h t+δt = h t + δ t ∇C s k (h t ), for some small δ t . Recall that, when training by batches, the Laplacian's integral takes the form (5). Here we chose that t k = kδ t , so that every time we switched to a new batch after a single time-step. In particular, we get t k t k-1 ∆C s k (h t )dt = ∆C s k (h t-1 )δ t in (5). Table 2 reports the different values of time discretisations δ t , used for the models that we trained. .6 2.00 × 10 -4 3 × 10 5

.2 6.67 × 10 -5 10 6

.06 2.00 × 10 -5

Figure 2 in the main text shows the relative error between the values of (2) obtained for the same model (N = 10 5 ) for two different time discretisations (δ t,1 = 2 × 10 -4 and δ t,2 = 2 × 10 -6 ). In both cases we trained for a time-interval of 2 × 10 -4 before switching batch. This means that for δ t,1 each batch was used for a single time-step, before switching to the next one, while for δ t,2 100 gradient-steps were performed on each batch before passing to the next one.

The values of (2) that we obtained for the different dimensions N are reported in Table 1, in the main text. The penalty log 25 is already taken into account for the values reported, as for each width the best bound has been chosen among the values obtained with 25 different time horizons. The optimal T reported in Figure 1 are the values of the best time horizon (among the 25 ones that we tried) that we obtained for each model's width.

  [0:T ] . . . δh denotes the line integral along the path h [0:T ] . The last term has a clear interpretation, as it quantifies how much h [0:T ]

Figure 1 :

 1 Figure 1: Behaviour of best time horizon T with the model dimension N , from experiments on MNIST for linear models trained to optimise a cross-entropy objective. The optimal T was found to go roughly as 1/N (cf. dotted line), preventing the bounds (1) and(2) from exploding, even in highly over-parameterised settings. Note that the results reported here were obtained via discretised gradient dynamics (cf. Section 5.4 and Figure2). We refer to Appendix E.3 for further details on the experiment.

Figure 2 :

 2 Figure 2: Relative error (B1 -B2)/B2, for a linear model with 10 5 parameters trained on MNIST via discretised gradient steps on a cross-entropy objective. Here, B1 and B2 are the values of the bound(2) obtained with the time-steps δt,1 = 2 × 10 -4 and δt,2 = 2 × 10 -6 respectively. The small relative error suggests that the corrections to the bound due to the time-discretised dynamics are likely to be negligible. We refer to Appendix E.3 for further details on the experiment.

Figure 3 :

 3 Figure 3: Behaviour of the bounds (1) and (2) for a linear model with 10000 parameters, trained with a linear surrogate loss, for different time horizons T . The data came from eight different Gaussian distributions on R 5 that had been assigned binary labels. The training dataset counted 500 datapoints and a held-out test dataset was used to evaluate the test errors reported in the plot. We refer to Appendix E.1 for further details.

Figure 4 :

 4 Figure 4: Behaviour of the bounds (1) and (2) for a linear model with 10000 parameters, trained with a quadratic surrogate loss, for different time horizons T . The learning task was the same as in Figure3. Again, we used a training dataset made of 500 datapoints and we evaluated the test errors on a held-out test dataset. The parameters α and β in Cs were both set to 1. We refer to Appendix E.2 for further details.

Figure 5 :

 5 Figure 5: Behaviour of the bounds (1) and (2) at different time horizons T , for a linear model with 10 5 parameters trained with the cross-entropy loss on MNIST. The training dataset counted m = 60000 datapoints and a held-out test dataset was used for the test errors reported. The plot also shows the values of the penalty factor (log ρ 0 (h 0 ) ρ 0 (h T ) + log 25 √ m δ +

Figure 6 :

 6 Figure 6: For the same experiment illustrated in Figure 5, comparison of the log density term log ρ 0 (h 0 ) ρ 0 (h T )

  Figure 8 compares the values of the bound (2) for different dimensions N .

Figure 7 :

 7 Figure 7: For the experiment with the linear surrogate loss, comparison of the evolution of the bound (2) for different sizes N of the the model. See Section 6.1 and Appendix E.1 for more details.

Figure 8 :

 8 Figure 8: For the experiment with the quadratic surrogate loss, comparison of the evolution of the bound (2) for different sizes N of the the model. See Section 6.2 and Appendix E.2 for more details.

  Figures 5 and 6 refer to the model with size N = 10 5 .

Table 1 :

 1 Bounds on MNIST

	Dimension N Bound (2) Test Error
	10 3	.3251	.2049
	3 × 10 3	.3034	.1967
	10 4	.2786	.1680
	3 × 10 4	.2795	.1541
	10 5	.2751	.1602
	3 × 10 5	.2812	.1534
	10 6	.2757	.1568

  Let s and s be independent training dataset, with s ∼ µ m . Let C s and C s be two twice differentiable training objectives, fix T > 0 and t 0 ∈ [0, T ). Consider an algorithm that, for an initial state h 0 , follows the dynamics ∂ t h t = -∇C s (h t ) for t ∈ [0, t 0 ], and ∂ t h t = -∇C s (h t ) for t ∈ [t 0 , T ], and finally outputs h T . Then, for any measurable function Ψ : R 2 → R, with probability at least 1 -δ on the random draw (s, h 0

	T
	∆C s ( ĥT -t0+t )dt ,
	0
	which is (6). We can hence state the following variant of Proposition 2 for data-dependent priors.
	Proposition 3.

Table 2 :

 2 Time discretisation

	Dimension N Final Time	δ t
	10 3	60	2.00 × 10 -2
	3 × 10 3	20	6.67 × 10 -3
	10 4	6	2.00 × 10 -3
	3 × 10 4	2	6.67 × 10 -4
	10 5

To be rigorous, one should actually require that s → ρ s is a Markov kernel.

To be fully rigorous, one should explicitly assume that the solution t → ht of the dynamics is well defined on [0, T ], for all initial conditions h 0 ∈ H.

For ∂tht = -∇Cs(ht) the single-step error is of order O(ε 2 ).

(y, y ) = 0 if y = y ; 1 otherwise.

The penalty factor log 50 (see end of Section 5.2) is already included in the values reported.

In the unregularised case α = 0, the matrix inversion Θ -1 , necessary to compute h T , can be source of significant computational errors. Introducing a factor α > 0 makes the algorithm more stable.

Acknowledgments: This work has been partially supported by the UK Engineering and Physical Sciences Research Council (EPSRC) through the grants EP/R513295/1 (DTP scheme) and CoSInES EP/R034710/1. The authors would like to thank Umut S ¸imşekli, Tyler Farghly, and Patrick Rebeschini for the valuable comments and suggestions.