Boubacar Traore
email: boubacar.traore@univ-rennes1.fr

Jacky Even

Laurent Pedesseau

Mikaël Képénekian

Claudine Katan
email: claudine.katan@univ-rennes1.fr

Band gap, effective masses and energy level alignment of 2D and 3D halide perovskites and heterostructures using DFT-1/2

We revisit the Slater half-occupation technique within the DFT-1/2 method to provide improved accuracy of 2D and 3D halide perovskites band gaps at a moderate computational cost. We propose an electron removal scheme from the halide states that drastically improve the predicted band gaps of 2D compounds. Concurrently, we compute the effective masses of the considered structures and show that DFT-1/2 describes them with a nice degree of accuracy when compared to available experimental data. Moreover, we assess the suitability of DFT-1/2 to compute the energy level alignments of this family of compounds. We compare the results in light of those we predict using the hybrid functional HSE06 and the many body perturbation theory within the G0W 0 approximation. We discuss the limitations of our refined DFT-1/2 scheme, which tends to overestimate the downshift of the absolute valence energy levels of the layered perovskite systems.

We anticipate that our results would be useful to initiate benchmark studies and investigate large systems such as heterostructures for which other approaches may be out of reach.

I. INTRODUCTION

Metal halide perovskites (MHP) exhibited a giant leap in their solar cell efficiencies and optoelectronic properties over the past years 1,2 . While the highest efficiencies have so far been achieved in their 3D structures, their 2D counterparts are adding interesting features related to their structural diversity, stability, photophysical and quantum confinement properties.

Understanding the latter photophysical processes is key to futher exploit these materials to their full potential. Moreover engineering complex heterostructures of 3D and 2D MHP is one of the important experimental trend nowadays 3 . In both 2D and 3D compounds, the band gap value and its nature precondition their aptness to favorably perform as solar cells or light emitting diodes. Indeed, band gap is an essential property of semiconductors and optoelectronic materials.

It is fundamentally defined as the difference between the ionization energy I and the electron affinity A : E g = I -A, where E g is the fundamental or electronic band gap. Experimentally, it is measured via spectroscopic techniques such as photoemission spectroscopy (PES) and inverse photoemission spectroscopy (IPES) 4 . From the theoretical side, predicting band gaps is quite challenging. With regards to the condensed matter community, Kohn-Sham (KS) density functional theory (DFT) 5,6 is among the most popular computational methods to investigate the fundamental properties of materials. Within DFT, the computed band gap is obtained as the difference between the eigenvalues ( ) at the conduction band minimum (CBM) and the valence band maximum (VBM) using E KS g = CBM -VBM . However, It is known that the band gap is severely underestimated within DFT at the LDA (Local Density Approximation) or GGA (Generalized Gradient Approximation) levels of theory because of the self-interaction error in the occupied states 7 and the absence of the derivative discontinuity ∆ xc 8,9 .

Various levels of theory have been developped to improve the band gap prediction of solids such as the many-body perturbation theory within the GW approximation 10,11 , hybrid functionals [12][13][14][15] and approximations using meta-GGA based potentials 16,17 , self-interaction error correction techniques 7 , DFT+U 18,19 and DFT-1/2 20 . Despite GW being formally among the most exact approximations for predicting fundamental band gaps of materials, its application to large systems (exceeding hundreds of electrons), with accurately converged results, becomes extremely cumbersome. With hybrid functionals, large systems start to be accessible but at the cost of high computational resources. Hence, efficient and computationally cost effective band gap correction schemes are highly desirable. Among those approaches is the Slater half occupation technique in the so-called DFT-1/2 method. So far, DFT-1/2 has been applied to predict the band gaps of a variety of materials ranging from classical semiconductors, metal oxides and 3D hybrid perovskites with a degree of accuracy comparable to GW [20][21][22][23][24][25][26][27][28][29] . Recently, Xue et al. proposed shLDA-1/2 (shell-LDA-1/2) to improve DFT-1/2 in the case of covalent semiconductors 30 . This resulted in the improvement of the band structure of Ge and the definition of some rules in the choice of the ratio of electron removal (e.g. 1/2e -or 1/4e -) to predict the band gap of covalent semiconductors.

The appealing feature of DFT-1/2 is its relative accuracy and efficiency coupled with a computational cost that is similar to standard LDA or GGA. It takes its root from the Slater transition state theory which was generalized to solids through Janak's theorem 31

∂E(f α ) ∂f α = (f α ) (1)
where E(f α ) is the total energy with respect to the occupation (f α ) for orbital α and (f α )

is the eigenvalue at f α . By integrating Eq. 1 between the state at which one electron is removed (f α = -1) and the neutral state (f α = 0), it is trivial to show that

E(0) -E(-1) ≈ (- 1 2 ) ( 2
)
where the equality assumes a linear variation of the Kohn-Sham eigenvalues with the occupation 32 . Eq. 2 tells that the eigenvalue at the half-ionized state is the negative of the ionization energy. This is at the heart of DFT-1/2 scheme and is best used when applied to the highest occupied molecular orbital (HOMO) or the valence band maximum (VBM) state. In general, the half-ionization is performed on the electronic orbitals that are present at the VBM or HOMO state. The latter information can be obtained from the projected density of states and/or from a symmetry analysis In other words, one should not blindly proceed with DFT-1/2 without a prior knowledge of the band edge electronic states. For instance, in the case of Ge, only half-ionized 4p states are considered because of their contribution to the VBM. 20 .

To benefit from Eq. 2 in a DFT self-consistent scheme, a self-energy (S α ; 2S α = ∂ α ∂fα ) term is introduced, which can be interpreted as the energy necessary to localize the hole 20 .

Integrating S α and combining it with Eq. 2 leads to

-S α = -1/2 0 ∂ α ∂f α df α ≈ (- 1 2 ) -(0) ⇒ E(-1) ≈ E(0) -(0) + S α (3)
where the self-energy is assumed to be independant of the occupation 33 . The self-energy potential (V S ) is chosen such that

S α = dr 3 ρ α (r)V S (r) ≈ (0) -(- 1 2 ) and V S (r) ≈ -V (- 1 2 , r) + V (0, r) (4)
ensuring that Eq. 2 is recovered. The self-energy potential is approximated by the difference between the all-electron potentials from the neutral single atoms (V (0, r)) and the corresponding half-ions (V (-1 2 , r)). In practice, the potential V S (r) is not computed on a solid but on an isolated atom where the orbital α is located. In other words, the potentials V (0, r) and V (-1 2 , r) obtained on the isolated atom of interest are used to compute V S (r) using the bottom expression in Eq. 4.

Adding V S to the atoms in an infinite crystal leads to divergence as it falls off in the form of 1/r. Therefore, it is trimmed with a step function using:

Θ(r) =    (1 -( r rcut ) 8 ) 3 , r ≤ r cut 0, r ≥ r cut V S (r) → V S (r) = V S (r) × Θ(r) (5)
The trimming parameter or the cutoff radius r cut is obtained variationally and is chosen so as to maximize the band gap. The power of 8 in Eq. 5 already provides good results but it may vary. To proceed with DFT-1/2 technique in a pseudopotential based DFT code, one has to identify a priori the electronic states contributing to VBM (or HOMO) and generate the pseudoptentials containing the self-energy potential V S correction for these states at the atomic level.

In the context of hybrid halide perovskites, DFT-1/2 has already been implemented to improve the prediction of the band gap of 3D compounds, leading to a good agreement with experiment 26,29,34 . However, the DFT-1/2 pseudopotentials, which are generated for 3D perovskites present insufficient performance on 2D layered MHP (vide infra). Moreover, to the best of our knowledge, there has not yet been any report on the optimization of DFT-1/2 for 2D hybrid halide perovskites. Here, we revisit DFT-1/2 method for halide perovskites and propose, using symmetry considerations, an electron removal scheme on the halide states that drastically improves the band gaps of 2D compounds. At the same time, we show that the effective masses are also corrected. Finally, we propose a critical assessment on the ability of DFT-1/2 to predict absolute valence energy alignments or ionization energies in light of the results that we obtained from hybrid functional (HSE06) calculations and quasi-particle correction within the GW approximation.

II. COMPUTATIONAL METHODS

First-principles calculations were performed using the SIESTA 35,36 code based on a basis set of finite-range of numerical atomic orbitals. Norm-conserving Troullier-Martins pseudopotentials were used for each atomic species to account for the core electrons. 37 1s 1 , 2s 2 2p 2 , 2s 2 2p 3 , 3s 2 3p 5 , 4s 2 4p 5 , 5s 2 5p 5 and 5d 10 6s 2 6p 2 were used as valence electrons for H, C, N, Cl, Br, I and Pb respectively. Polarized Double-ζ (DZP) basis sets with an energy shift of 200 meV and a real space mesh grid energy cutoff of 400 Rydberg were used for the calculations. We used PBE and the revisited Slater half-occupation technique in the so-called DFT-1/2 20 implementation for the band structure calculations. To proceed with DFT-1/2, we applied the -1 2 e -correction scheme to the orbitals that form the VBM state of the relevant compounds. This means that for the concerned atoms, we generated DFT-1/2 corrected pseudopotentials following the bottom expression of Eq. 4. Again, we stress that one needs to know the electronic states forming the VBM or HOMO state before using DFT-1/2. Apart from that, the main points of the theory are as described in the previous section. Spin-orbit coupling (SOC), in its on-site approximation 38 , was taken into account in all calculations unless otherwise stated. Sufficient k-point grids were used to sample the different Brillouin zones. We computed the effective masses using the parabolic approximation around the band extrema.

For the 2D/3D heterostructure relaxation, we used the van der Waals density functional with C09 exchange 39 within the van der Waals DF2 flavor to describe the exchangecorrelation term (we modified SIESTA version "master-post-4.1-251" to include DF2 flavor). This tuning was suggested to improve C09 based geometry optimization over that with the original DF1 flavor 40 . The in-plane lattice constants of the 3D part were imposed for the whole heterostructure and were not allowed to relax during the 2D/3D interface structural optimization. The c-lattice constant along z was allowed to relax along with the atomic positions. In order to mimic the structural properties of 3D bulk materials at room temperature, we used Cs as a pseudotom instead of FA (FA * PbI 3 41 , FA * = FA replaced with Cs) and consider a thickness of 9 octahedra. The 2D structure is the room-temperature BA 2 PbI 4 (BA = Butylammonium) 42 . The whole heterostructure consists of 314 atoms and 1320 valence electrons. We computed the heterostructure band offsets using the Hartree Potential alignment method as described in Ref. 43.

For the absolute valence energy level alignment of the 3D systems, we used slabs constructed from SIESTA relaxed pseudo-cubic CsPbI 3 bulk structure derived from the roomtermperature orthorhombic compound reported in Ref. [ 44]. The relaxed lattice constant a = 6.297 Å (The experimental pseudo-cubic lattice constant a = 6.201 Å).

Hybrid functional (HSE06 15 ) and G0W 0 11 calculations were performed using VASP 45,46 within the Projector Augmented Wave (PAW) formalism 47,48 . The wave functions were expanded using a plane-wave basis set with an energy cut-off of 400 eV. We used 4×4×4 and 4×4×1 Γ-centered k-grids to sample the Brillouin zones of the bulk (in HSE06 and G0W 0) and slab (in HSE06) systems, respectively. We used the Adaptively Compressed Exchange (ACE) Operator 49 to reduce the computational cost of HSE06 calculations. We used 960 empty bands for G0W 0 calculations and the convergence of the quasi-particle energies within the approximation was considered (vide infra). Spin-orbit coupling effect was included in all the calculations.

Except for the 2D/3D heterostructure and pseudo-cubic CsPbI 3 compound, the experimental crystal structures without any geometry optimization were used for all calculations.

The following structures were considered with the space groups provided in parenthesis:

3D systems: MA * PbI 3 (I4/mcm) 50 with MA * = MA (methylammonium) replaced with Cs to avoid the issue of dynamical disorder of MA at high temperature phases 51 , MAPbI 3 (Pnma) 52 , FA * PbI 3 (Pm-3m) 41 , MA * PbBr 3 (Pm-3m) 53 , MAPbBr 3 (Pnma) 54 , FA * PbBr 3 (Pm-3m) 55 , CsPbI 3 (Pnma) 56 , CsPbBr 3 (Pnma) 57 , CsPbCl 3 (Pm-3m) 58 , RbPbI 3 (Pnma) 56 , MASnI 3 (I4cm) 59 , CsSnI 3 (Pnma) 60 .

2D systems: The references of the structures used or considered are provided in Table VII. In 2D systems, the subscript n in the formula corresponds to the number of inorganic octahedra or inorganic layer thickness in the structure as, for instance, in BA 2 MA n-1 Pb n I 3n+1 (BA = butylammonium).

III. RESULTS AND DISCUSSION

A. DFT-1/2 for 3D metal halide perovskites

We start by applying the DFT-1/2 method to 3D halide perovskites with the general formula AMX 3 where A is a monovalent organic or inorganic cation (A = Cs, MA

[MA=CH 3 NH 3 ], FA [FA=CH(NH 2 ) 2 ]
), X is a halide anion (X = I, Br, Cl) and M is a divalent metal cation (M = Pb, Sn ...). In order to use the method, one has to know a priori the states contributing to the VBM of the material under study. In 3D halide perovskites, it is quite established, from first-principles calculations 61 and predicted from symmetry analysis 62 that the VBM of 3D halide perovskites is made of an anti-bonding hybridization between mainly X(p) and M(s) states (Figure 1a). Considering the room temperature structure of MA * PbI 3 in the tetragonal phase (I4mcm) 50 , Figure 1 shows the optimal r cut values for I and Pb that maximize the band gap of this compound (see Eq. 5). Comparing the DFT-1/2 band gaps (with SOC) at the optimal r cut for I (1.40 eV at r cut ∼3.76 Bohr) and Pb (1.52 eV at r cut ∼2.20 Bohr + optimized r cut from I) with the one obtained at the PBE+SOC level (0.35 eV), we observe that the major part of DFT-1/2 correction is due to the I(5p) states. This is consistent since the main contribution to the VBM of these compounds is due to the halide p states. It further supports the idea that DFT-1/2 correction on the metal cation states could generally be neglected in most semiconductors and oxides (ZnO seems to be an exception where Zn (3d) states also become important) 33 . Here, we included the correction on the metal cation (Pb and Sn) states since it brings an additional improvement as compared to experiments (vide supra with the additional 0.12 eV in MA * PbI 3 ).

Following the same recipe as for MA * PbI 3 , we inspected DFT-1/2 scheme with other 3D compounds including Sn, Br and Cl and the optimal r cut values are shown in Table I. Noteworthy, an optimized r cut for an element in a specific compound is generally transferable to the same element in a different chemical environment using a similar r cut value with only minor changes 20 .

Table II summarizes DFT-1/2 corrected bands gaps for a range of 3D MHP. Previous results obtained using ppTB-mBJ potential are also included for comparison 63 . Overall, DFT-1/2 brings a substantial improvement to the band gaps of these compounds as compared to PBE. Among the investigated compounds, DFT-1/2 predicted band gaps present a mean absolute error (MAE) of 0.17 eV as compared to the experimentally reported values.

When compared to ppTB-mBJ, both methods perform well with the deviation from the experiment being smaller for ppTB-mBJ (MAE of 0.11 eV). effect is not taken into account here since it is not necessary for the optimization of r cut .

Considering the effective masses summarized in Table III, PBE severly underestimate them by more than 50% as compared to available experimental data with a mean absolute percentage error (MAPE) of about 50.1%. Using ppTB-mBJ, the effective masses are improved as compared to PBE but tend to be overestimated with regards to the experimental values (MAPE ∼ 46.5%) as discussed thoroughly in a prior report 63 . Computed DFT-1/2 are also slightly on the high side (MAPE ∼ 20.2%) but with a better agreement with the available experimental data. This is consistent with previous reports that DFT-1/2 not only improves the prediction of band gaps but also the band curvatures i.e. the effective masses 20 . Here, we consider 2D MHP whose optical and photophysical properties are very often dominated by excitons 83 . Notewothy, the latter excitonic effects are not captured within DFT-1/2, which it is not intended for. We apply the 3D optimized DFT-1/2 pseudopotentials to compute the band gaps of a wide range of 2D halide perovskites. More specifically, this corresponds to removing 1 2 e -from X(p) and M(s) states. The results are summarized in Table VII. We notice that DFT-1/2 brings some improvements to the band gaps of 2D compounds as compared to PBE. However, the computed values are still largely underestimated when compared to experiment. For instance, DFT-1/2 predicts 2.20 eV for the 2D n=1 BA 2 MA n-1 Pb n I 3n+1 which is still 0.6 eV lower as compared to the continuum threshold from magneto-absorption experiments (2.80 eV) 84 . This becomes worse for the n=5 compound whose computed band gap (1.36 eV) is even lower than that of the room-temperature 3D MAPbI 3 structure (1.52 eV from Table II). Clearly, using the 3D optimized DFT-1/2 parameters to investigate 2D perovskite structures still underestimates band gaps.

The inferior performance of DFT-1/2 band gap prediction for 2D MHP, considering only the correction on halide p states and Pb (6s), can to a certain extent be related to the symmetry change of VBM/CBM when going from 3D to 2D. This arises from the fundamental anisotropy of the system, and is more specifically related to the presence of apical halide atoms in 2D compounds. Indeed, from the symmetry analysis implemented for 2D halide perovskites 2,85 , an additional electronic coupling appears at the VBM state: a coupling between apical X(s) and M(s) as shown in Figure 2. This additional X(s) states expected at the VBM of 2D halide perovskites, though in a minor contribution, is absent in 3D ones.

Because of the emergence of X(s) states from the symmetry analysis, we refine DFT-1/2 pseudopotentials by removing additional 1 4 e -from the s orbitals of the halide atoms. In fact, with removed 1 2 e -from X(s) orbitals, we find that the band gaps are significantly overestimated. This is not surprising since the contribution of s states to the VBM of these 2D compounds remains limited. Hereafter, DFT-1/2 with the additional 1 4 e -removed from X(s) states is referred as xs-DFT-1/2 and the corresponding r cut are also reported in Table I.

At this stage, we emphasize that this adjustment (additional 1 4 e-) is made in the spirit of the DFT-1/2 scheme: at its initial stage, the method may need to be readjusted to account for the chemical environment of atoms in a specific family of compounds. Noteworthy, once parameterized, the method becomes predictive a postoriori for similar systems where atoms exhibit similar chemical environments. For instance, we have adjusted the additional amount of electrons to be removed from X(s) within xs-DFT-1/2 on BA 2 PbI 4 system only.

The approach is then applied to all the other 2D systems with much improved computed band gaps as summarized Table VII. Table VII compiles the results obtained using the refined xs-DFT/-1/2 for the considered 2D MHP. When compared to DFT-1/2, we notice a dramatic improvement in the computed band gaps using the refined method as compared to the experimentally reported values. This is especially true for the iodine and bromine based 2D perovskites. For instance, the continuum thresholds of the Ruddlesden-Popper (RP) phases of 2D BA 2 MA n-1 Pb n I 3n+1 are nicely captured by xs-DFT-1/2. For n=4, both levels of theory lead to computed band gaps that deviate from the monotonous trend of decreasing experimental band gaps with increasing n-value. This is likely due to uncertainties in the experimental crystallographic structures.

In the case of chlorine based 2D perovskites, it appears that the overestimation of the predicted band gaps is more pronounced using xs-DFT-1/2 refinement. This is possibly related to an over-correction of Cl(3p) orbital eigenvalues as a result of removing the additional 1 4 e -from its 3s states leading to larger band gaps. We note, however, that reported experimental band gaps in Table IV for chlorine compounds are optical band gaps and a fair comparison with theory should rather rely on electronic band gaps, which are lacking. This discrepancy needs more investigations on both theoretical and experimental sides. Nevertheless, xs-DFT-1/2 brings the computed band gaps closer to experiment as compared to the normal DFT-1/2. In many of these 2D compounds, applying DFT-1/2 correction on only the inorganic states (X and M states) causes a shift of the latter states leaving behind the molecular ones. These molecular states may sometimes be spuriously located inside the band gap (underlined data in Table VII). This is an artifact of the uncorrected molecular states (Figure 3). In those cases, it is useful to implement the DFT-1/2 method to the molecular states as well. Using the optimized r cut for N and C, reported in Table I, allows to recover the correct band gaps determined by the inorganic entities (Figure 3). This procedure may be justified by the fact that these layered systems can be conceptualized as composite materials or heterostructures where the inorganic and organic frameworks exhibit specific band alignments and their band edges may be aligned differently depending on the nature of the molecular spacers 109,110 .

Incidently, even after correcting both the inorganic and organic states for (PMA) 2 PbCl 4 , some molecular states remain at the band egde implying that this compound may present a type-II band alignment as we suggested in a previous work 110 .

In Table V, we report the average of the in-plane effective masses of various 2D MHP.

Similar to what has been shown for 3D MHP, DFT-1/2 leads to a sizeable increase in computed effective masses as compared to PBE. The additional correction on X(s) states leads to a further increase and in turn a better agreement with available experimentally reported values. This shows the importance of predicting more accurate band gaps for computing effective masses that are in better agreement with experiment. C. Is DFT-1/2 more empirical than HSE06 for metal halide perovskites?

We note that within our xs-DFT-1/2 refinement, the total number of removed electrons [ 3 4 e -→ 1 2 e -(X(p)) + 1 4 e -(X(s))] somehow violates the fundamental formalism behind the DFT-1/2 theory (Eq 2). In principle, this number should be 1/2e -for a given atom. However, the refinement is in part justified by the fact that the method remains semi-empirical, with necessary adjustments to take into account different symmetries and chemical environments. In fact, the developpers of DFT-1/2 alluded to this in their original paper 20 .

They had to change the number of removed electrons for Si (3p) from 1 2 e -to 1 4 e -to reach a good match with its experimental band gap. This readjustment was motivated by the covalent nature of Si-Si bond. But for Ge(4p), the same column as Si, they kept 1 2 e -. More recently, DFT-1/2 was revisited for conventional semiconductors in order to obtain better quantitative estimates of band gaps leading to the so-called shell-LDA-1/2 (shLDA-1/2) 30 .

In shLDA-1/2 method, some electron removal schemes (e.g. 1/2e -or 1/4e -) and different cutoff radii for the relevant orbitals were adopted to include the effect of the chemical environment. This further illustrates the intrinsic limitation and the semi-empirical character of the method as previously discussed 111 . Anyhow, DFT-1/2 correction assumes a priori knowledge of the orbitals composing the VBM, by virtue of which the full ab initio like picture is clearly lost. It is worth noting that the strategy based on parameters tuning for band gap prediction is not unique to DFT-1/2. Meta-GGA based exchange potential in the TB-mBJ implementation also relies on parameters tuning 112 . Even the more popular hybrid functionals, such as HSE06, also include a parameter, namely the amount of the exact exchange (α). Interestingly, different amounts of exact exchange are also needed to accurately predict the band gaps of 3D and 2D MHP 113,114 . Here, we further inspect the impact of the amount of exact exchange on computed band gaps of 3D and 2D perovskites. The results are summarized in Figure 4 and Table VI. Our results clearly show that a unique value of exact exchange α does not allow to accurately predict the band gaps of both 3D (α ∼ 50%) and 2D (α ∼ 60%) halide perovskite compounds. Moreover, considering the 2D MHP, one can see that the same amount of exact exchange does not peform in the same manner for iodides compounds as for bromides. Thus, the issue of accurate band gap prediction based on first-principles remains an open problem in condensed matter systems. TABLE VI: Computed band gaps of 2D and 3D halide perovskites at the HSE06 level of theory using an increasing amount of exact exchange (α). The superscripts "opt" and "el" stand for optical and electronic band gaps, respectively.

D. Assessing ionization energies or absolute valence energy levels using DFT-1/2

Having critically discussed the accuracy of DFT-1/2 in predicting the band gaps of halide perovskites, we further inspect the behavior of the method in predicting the absolute valence energy levels (or the ionization energies) of these compounds. We compute the absolute valence energy levels of the different systems using the Hartree potential alignment method (Figure 5a) with the vacuum level set as the reference [117][118][119] . Indeed, by breaking the periodicity of the material along one direction (here, along z) by adding vacuum to it, forming a semi-infinite slab, allows to have the zero of the potential to be a well-defined quantity 120,121 . Here, the absolute valence energy level E abs v is obtained using:

E abs v = E bulk V BM + ∆V H + ∆E vacuum (6)
where E bulk V BM is the eigenvalue of the valence band maximum of the bulk system, ∆V

H = V slab H -V bulk H
is the shift of the Hartree potential between the slab and the bulk systems, and ∆E vacuum = 0 -E vacuum is the vacuum level shift with respect to zero energy level. When considering the quasi-particle correction within the single-shot G0W 0 approximation, we use

E abs v = E bulk V BM + ∆Σ V BM + ∆V H + ∆E vacuum (7)
where ∆Σ V BM is the shift in the eigenvalue of VBM due to the Self energy quasi-particle correction within G0W 0. Here, ∆V H is obtained from the starting mean-field theory used for the ground state calculation (PBE or HSE06). This has a minor effect on the computed E abs v since the potential shift depends on the ground state densities which are relatively well computed at the PBE or HSE06 level 121 . Indeed, considering CsPbI 3 slab, Figure 5a shows that ∆V H is comparable at the PBE and HSE06 levels, which supports the use of ∆V H from the DFT level for the absolute valence energy level alignment with the quasi-particle energies obtained from the G0W 0 approximation. The convergence of the quasi-particle energies at VBM and CBM from the G0W 0 approximation are shown in Figure 5b. For the conduction band, we added to E abs v the theoretically computed band gap of the bulk compound within the used level of theory. This then gives, ∆V H is similar in both cases which partly justifies the use of ∆V H from the DFT level for the absolute valence energy level alignment using the quasi-particle energies from G0W 0 approximation. b Convergence of the G0W 0 quasi-particle energies at VB/CB of bulk CsPbI 3 with respect to the response function cutoff and the number of empty bands. Note that, by default, VASP uses 2/3 of the plane wave cutoff for the energy cutoff of the response function.

E abs c = E abs v + E g,bulk (8)
We use CsPbI 3 slabs constructed from the relaxed pseudo-cubic CsPbI 3 bulk structure as model systems. We consider two different terminations of the CsPbI 3 slab namely CsI and PbI 2 as shown in Figure 6a. CsI-termination is reported to be the most stable surface from a previous theoretical calculation 122 . Figure 6b summarizes the predicted absolute valence energy levels for the CsPbI 3 slabs at different levels of theory. In all cases, PbI 2 -termination presents a deeper absolute valence energy level as compared to CsI-termination for these unrelaxed slab surfaces. This can be explained by the larger surface dipoles present at the PbI 2 -terminated surface, which we have thoroughly discussed in this work 123 . DFT-1/2 predicts E abs v values that are closer to those predicted at the G0W 0@HSE(50%) level. Here, G0W 0@PBE presents the least negative values followed by HSE06(50%) showing the impact of the starting point for obtaining more accurate quasi-particle excitation energies 119 .

DFT-1/2 tend to predict deeper E abs v values as compared to the other levels of theory. In comparison to G0W 0@HSE(50%), DFT-1/2 differ by 0.36 eV and 0.38 eV for CsI and PbI 2terminated surfaces, respectively. Unfortunately, a wide range of experimental values have been reported for CsPbI 3 ranging from 4.68 eV for β-CsPbI 3 124 to 5.2-5.7 eV [125][126][127][128] reaching as high as 6.22 eV 129 . This hampers a proper assessment of the level of theory based on experimental data. The disparity of experimental results stems from several origins among which are process conditions, the used substrates, surface morphologies, surface terminations, crystal purity and measurement methods 130 . Based on the relatively good agreement between DFT-1/2 and G0W0@HSE(50%) results, DFT-1/2 appears to be a promising method for computing the absolute valence energy levels. This is not really surprising since the half occupation technique at the heart of DFT-1/2 leads to the negative of the ionization energy and thus suitable for VBM or HOMO states (vide supra, Eq. 2). for the two compounds using the normal DFT-1/2 fairly agree with the reported values by Silver et al., which suggests that DFT-1/2 would be more reliable in predicting band offsets and energy levels related to the valence bands (more consistent with Eq. 2 to predict the ionization energies). The additional 1 4 e -removed from the halide s states within xs-DFT-1/2 tends to push down the VBM states by opening the band gap but at the same time deepens the absolute valence energy level impairing the agreement with experiments. We believe that going beyond the prediction of band gaps, the refined xs-DFT-1/2 for 2D halide perovskites might present some caveats when used to infer band alignments but more investigations especially on the experimental level are needed to conclude.

E. Band alignment of a 2D/3D heterostructure

We further apply the DFT-1/2 method to compute the band alignment in a sizeable halide perovskite 2D/3D heterostructure. Here, we consider the 3D FA * PbI 3 interfaced with BA 2 PbI 4 2D perovskite as explained in Section II (Figure 7a). Figure 7b,c show the computed band offsets using PBE and DFT-1/2, respectively. We predict a Type-I band alignment from both levels of theory. In particular, the 2D layer presents a barrier to holes collected from the 3D region. Referring to our discussion on the assessment of DFT-1/2 in predicting the ionization energies (vide supra), we postulate that our computed valence band offsets from the heterostructure would be more reliable than PBE and closer to experiment.

Unfortunately, for such an extended system, implementation of G0W 0 calculations are computationally prohibitive. Through this example, we highlight the computational efficiency of DFT-1/2 including spin-orbit coupling effect, which is of similar computational cost as a standard PBE calculation. In this section, we briefly discuss on the computational merits of DFT-1/2 as compared to more sophisticated levels of theory. This includes the hybrid functional HSE06 and the many-body perturbation theory within G0W 0 approximation. We consider the Coulombhole and screened-exchange (COHSEX) approximation (Phys. Rev. B 74, 045102, 2006), which is the static limit of the GW approximation and appears to be more computationally efficient. We use the relaxed CsPbI 3 pseudo-cubic (a=6.297 Å) model to structure to compare the computational efficiency of the different levels of theory. Table VII summarizes the performance of G0W 0@PBE, COHSEX@PBE, HSE06(α=50%) and DFT-1/2 to compute the band gap of the pseudo-cubic CsPbI3 structure. Because of the static approximation to the self-energy, COHSEX presents a better computational efficiency as compared to the standard G0W 0 calculation with a computing time that is about 2.5 times lower. However, DFT-1/2 is much more efficient as compared to either type of GW approximation or HSE06.

For instance, DFT-1/2 presents a computational time that is orders of magnitude lower than that of COHSEX implementation within VASP. Besides, this COHSEX calculation requires much more memory (309.3 Mbytes), which might present an additional bottleneck for GW type of calculations as compared to DFT-1/2 (55.7 Mbytes using VASP). Indeed, DFT-1/2 presents a similar computational efficiency as for a standard PBE. Noteworthy, DFT-1/2 computed band gaps reach a similar accuracy as those obtained from GW calculations. d The computation of the dielectric matrix usually requires the use of a large number of empty states. Within VASP, the of bands at the dielectric matrix step needs to be the same as in the GW step for computing the screened exchange self Σ. Hence, our COHSEX approximation also uses a large number of bands.

G. Conclusions

In summary, we present a generic inspection of the performance of DFT-1/2 for computing the band gaps, effective masses and band alignments of metal halide perovskites. We propose a refinement of DFT-1/2 method within the xs-DFT-1/2 approach that improves the prediction of band gaps and effective masses of 2D MHP. DFT-1/2 approach is appealing and brings orders of magnitudes improvement in the computational efficiency but at the same time produces band gaps that are comparable to those obtained from more computationally demanding state-of-the art HSE and GW levels of theory. The inspection of the energy level alignment further shows the capabilities of DFT-1/2 to study band alignments and predict band offsets in heterojunctions.

When compared to ppTB-mBJ, DFT-1/2 computed band gaps present a slightly larger mean absolute error as compared to experimentally reported values. However, DFT-1/2 better describe the band dispersions leading to improved effective masses as compared to ppTB-mBJ.

Regarding xs-DFT-1/2, specifically designed for layered perovskites, it improves the band gaps and effective masses as compared to normal DFT-1/2. However, beyond the computation of band gaps and effective masses, it tends to push down the VBM states by opening the band gap but at the same time deepens the predicted ionization energy values impairing the agreement with experiments. Hence, xs-DFT-1/2 might present some caveats in the context of computing band alignments but more investigations especially on the experimental level are needed to conclude.

The amount of additional electrons to be removed beyond 1 2 e -of standard DFT-1/2 is more empirical and is somehow determined through adjustments for a specific compound presenting peculiar chemical environments. Once adjusted, the method becomes predictive a posteriori when applied to other compounds within the same family.

Altogether, DFT-1/2 is an efficient computational method that can be useful in studying the electronic properties of large structure compounds with a good degree accuracy. Moreover, it can be interesting in large benchmark studies where comparison can be attained on sizeable datasets.

FIG. 1 :

 1 FIG. 1: Optimization of DFT-1/2 r cut values of I and Pb using the tetragonal MA * PbI 3 . a Band structure of MA * PbI 3 in the tetragonal phase. The wavefunction, without SOC, at the VBM occuring at the Γ point is also shown in the figure. b Optimizing r cut for (left) the pseudopotential of I by removing 1/2 e -from its 5p states (right) Pb atom with the removal of 1/2e -from its 6s states using the optimized I(5p) atom from the left. Here the optimal r cut are ∼3.76 Bohr and ∼2.20 Bohr for I and Pb, respectively. Note that spin-orbit coupling (SOC)

FIG. 2 :

 2 FIG. 2: Symmetry analysis of 2D MHP 2,85 . (left) An example of a 2D MHP structure showing the apical and equatorial halide X atoms along with the central metal M atom. (right) Electronic coupling between X and M states expected from the symmetry analysis of 2D MHP.

FIG. 3 :

 3 FIG. 3: Evidencing the necessity to correct molecular states in some 2D MHP. a Band structure of (PEA) 2 (MA) 2 Pb 3 I 10 structure using xs-DFT-1/2 (left) without DFT-1/2 correction on molecular states and (right) with DFT-1/2 correction on molecular states. b Band structure of (PMA) 2 PbCl 4 structure using xs-DFT-1/2 (left) without DFT-1/2 correction on molecular states and (right) with DFT-1/2 correction on molecular states. The molecular states are shown in the dashed rectangular boxes.

FIG. 4 :

 4 FIG. 4: Impact of the amount of HSE06 exact exchange α on band gap prediction in halide perovskites. a Variation of the band gap of the room temperature 3D MA * PbI 3 MHP with the amount of exact exchange α. b Variation of the band gap of different 2D MHP with the amount of HSE06 exact exchange α. Horizontal dashed lines correspond to the experimentally reported electronic band gap values for the compounds of the same color.

FIG. 5 :

 5 FIG. 5: Theoretical methods for the absolute energy level alignment considering CsPbI 3 bulk and slabs. a Comparison of the Hartree potential shifts computed at the PBE+SOC and HSE06+SOC (α = 50%) levels of theory using VASP. The figures are obtained using CsPbI 3 slab.

FIG. 6 :

 6 FIG. 6: Slab structures and computation of the absolute valence energy levels of MHP at different levels of theory. a CsPbI 3 slabs with CsI and PbI 2 -terminations used for the computation of the absolute valence energy levels. b Comparison of the absolute valence energy level of pseudo-cubic CsPbI 3 with different terminations at different levels of theory. Here, PBE resuts are obtained using VASP to be consistent with HSE and G0W 0@PBE c Slab structures used to compute E abs v for BA 2 PbI 4 and BA 2 PbBr 4 . d Comparison of the absolute valence energy levels of BA 2 PbX 4 (X = I, Br) with BAX-termination at DFT-1/2 and xs-DFT-1/2 corrections.BAX stands for butylammonium halide (X) termination of the surface.

2 FIG. 7 :

 27 FIG. 7: Application of DFT-1/2 to compute the band alignment of a sizeable 2D/3D heterostructure. a Constructed 2D/3D FA * PbI 3 /BA 2 PbI 4 heterostructure model. In the heterostruture, both the 2D and 3D regions were terminated with BAI. b Computed band alignment using PBE without SOC effect. c Computed band alignment using DFT-1/2 including SOC effect.

b

  SIESTA : real space mesh cutoff = 300 Ry; k-point grid = 4×4×4 point; Basis set = DZP. c VASP : kinetic energy cutoff = 400 eV; k-point grid = 4×4×4.

  IV. ACKNOWLDGEMENTThis work was performed with funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 861985 (PeroCUBE). J.E. acknowledges financial support from the Institute Universitaire de France. This work was granted access to the HPC resources of [TGCC/CINES/IDRIS] under the allocations 2020-A0100911434 and 2020-A0090907682 made by GENCI.

TABLE I :

 I Cutoff radius r cut (in atomic unit) and half-ionized atomic orbitals within DFT-1/2 and xs-DFT-1/2, respectively.

	Atom r cut (a.u.)	Half-ionized orbital
			DFT-1/2
	Pb	2.20	6s
	Sn	2.20	5s
	I	3.76	5p
	Br	3.44	4p
	Cl	3.12	3p
		xs-DFT-1/2 for halides
	I Br Cl	4.00 3.44 3.12 DFT-1/2 for C, N 5p (-1 2 e -) + 5s(-1 4 e -) 4p (-1 2 e -) + 4s(-1 4 e -) 3p (-1 2 e -) + 3s(-1 4 e -)
	C	2.40	2p
	N	2.00	2p

TABLE II :

 II Band gap (in eV) of 3D MHP obtained with PBE, DFT-1/2 and ppTB-mBJ63 .The space group is indicated in parenthesis. SOC effects were included in all calculations. MAE and MAPE stand for mean absolute error and mean absolute percentage error, respectively with respect to the experimental values.

		This work This work	Ref. 63	Exp.
	3D MHP	PBE	DFT-1/2 ppTB-mBJ	Exp.
	MA * PbI 3 (I4/mcm)	0.35	1.52	1.43	1.48-1.61 64-69a
	MAPbI 3 (Pnma)	0.64	1.69	1.64	1.65 70b
	FA * PbI 3 (Pm-3m)	0.25	1.41	-	1.48 66 -1.52 69f
	MA * PbBr 3 (Pm-3m)	0.52	1.97	2.16	2.18-2.22 67,71,72c
	MAPbBr 3 (Pnma)	0.88	2.29	2.43	∼2.34 73g
	CsPbI 3 (Pnma)	1.91	2.79	3.33	3.17 74d
	FA * PbBr 3 (Pm-3m)	0.65	2.11	-	∼2.30 73g
	CsPbBr 3 (Pnma)	0.76	2.16	2.45	2.25 57,75c
	CsPbCl 3 (Pm-3m)	0.88	2.50	2.84	2.86 75c
	RbPbI 3 (Pnma)	1.87	2.77	3.30	3.14 76d
	MASnI 3 (I4cm)	0.22	1.27	1.55	1.2-1.4 59e
	CsSnI 3 (Pnma)	0.07	1.08	1.21	1.30 60c
	MAE (eV)	1.38	0.17	0.11
	MAPE (%)	68.8	7.5	5.9

a Measured at 293 K or room temperature. b Measured at 4.2 K. c Measured at room temperature. d Measured at 90 K. e Measured at 200 K.

f High temperature tetragonal phase according to the authors (140 K -160 K) g Measured at 300 K.

TABLE III :

 III Electron (m e ) and hole (m h ) effective masses in m 0 unit of 3D MHPs calculated with PBE, DFT-1/2 and ppTB-mBJ63 . µ = m e m h / (m e + m h ) is the reduced effective mass and m h = |m v |. For MAE and MAPE, we use the average of the computed/reported reduced effective masses for each compound.

	B. DFT-1/2 for 2D metal halide perovskites
			This work	This work		Ref. 63
	3D MHP	Direction	PBE	DFT-1/2	ppTB-mBJ	Exp.
			m v m e µ	m v m e µ	m v m e µ	µ
		Γ → [001] -0.07 0.07 0.04 -0.23 0.20 0.11 -0.27 0.21 0.12	0.104 77
	MA * PbI 3 (I4/mcm)	Γ → [110] -0.10 0.09 0.05 -0.26 0.26 0.13 -0.36 0.27 0.15	(155 K -190 K)
		Γ → [111] -0.09 0.09 0.05 -0.25 0.25 0.12 -0.33 0.25 0.14
	MAPbI 3 (Pnma)	Γ → [010] -0.14 0.11 0.06 -0.31 0.25 0.14 -0.45 0.27 0.17 0.15 78 , 0.12 79 (4.2 K) Γ → [101] -0.17 0.15 0.08 -0.33 0.31 0.16 -0.50 0.34 0.20
	FA * PbI 3 (Pm-3m)	R → [001] -0.06 0.06 0.03 -0.22 0.21 0.11 -Γ → [111] -0.06 0.06 0.03 -0.22 0.21 0.11 -	--	--	0.095 69 (140 K-160 K)
	MA MAPbBr 3 (Pnma)	Γ → [010] -0.16 0.15 0.08 -0.33 0.33 0.17 -0.50 0.37 0.21 0.13 78 (4.2 K) 0.117 69 (2K) Γ → [101] -0.19 0.19 0.10 -0.36 0.39 0.19 -0.59 0.44 0.25
	FA * PbBr 3 (Pm-3m)	R → [001] -0.12 0.13 0.06 -0.27 0.32 0.15 -Γ → [111] -0.12 0.13 0.06 -0.28 0.31 0.15 -	--	--	0.13 69 (160 K-170 K)
	CsPbBr 3 (Pnma)	Γ → [010] -0.13 0.13 0.07 -0.29 0.32 0.15 -0.44 0.37 0.20 Γ → [101] -0.16 0.16 0.08 -0.30 0.36 0.16 -0.53 0.44 0.24 (m h = 0.26 81 ) 0.126 (2K) 80
	CsPbCl 3 (Pm-3m)	R → [001] -0.15 0.17 0.08 -0.29 0.40 0.17 -0.40 0.42 0.20 Γ → [111] -0.15 0.17 0.08 -0.29 0.40 0.17 -0.40 0.42 0.20	0.202 (2K) 82
	MAE (m 0 )	-	0.07	0.03		0.06
	MAPE (%)	-	50.1	46.5		20.2

* PbBr 3 (Pm-3m) R → [001] -0.10 0.11 0.05 -0.25 0.30 0.14 -0.34 0.32 0.16 Γ → [111] -0.10 0.11 0.05 -0.25 0.29 0.13 -0.35 0.32 0.17

TABLE IV :

 IV Room temperature (otherwise specified) band gap (in eV) of 2D MHP obtained with PBE, DFT-1/2 and xs-DFT-1/2. Underlined data indicate that it was necessary to apply DFT-1/2 correction to the molecular states (containing C and N) as well. SOC effects were included in all calculations. The superscripts "opt" and 'el" stand for optical and electronic band gaps, respectively. For MAE and MAPE, when the experimental electronic band gaps are not reported, we use the optical ones, which add more uncertainties to the computed statistical values.

					X(ap)	M	X(eq)
					px, p y
					s
				Layer thickness PBE DFT-1/2 xs-DFT-1/2	Exp.
				n=1 42	1.28	2.20	2.84	2.80 84 el
				n=2 86	0.87	1.82	2.45	2.44 84 el
	BA 2 MA n-1 Pb n I 3n+1	n=3 86	0.64	1.65	2.29	2.26 84 el
				n=4 86	0.64	1.67	2.32	2.15 84 el
				n=5 87	0.39	1.36	2.01	2.08 84 el
	(C m H 2m+1 NH 3 ) 2 PbI 4	m=6 42 m=10 89	1.28 1.28	2.20 2.19	2.85 2.85	2.70 88 el (5K) 2.72 90 el
				n=1 91	1.12	2.05	2.70	2.27 91 opt
				n=2 91	0.72	1.67	2.31	1.99 91 opt
				n=3 91	0.16	0.87	1.46	1.73 91 opt
	PEA 2 MA n-1 Pb n I 3n+1	n=1 92 n=3 94	1.04 0.89	1.99 1.76	2.71 2.34	2.53 93 opt (10K) 2.1 95 opt
				n=1 96	1.02	2.02	2.69	2.23 96 opt
	(3AMP)MA n-1 Pb n I 3n+1	n=2 96	0.39	1.27	1.83	2.02 96 opt
				n=3 96	0.30	1.23	1.82	1.92 96 opt
				n=1 96	1.07	1.99	2.64	2.38 96 opt
	(4AMP)MA n-1 Pb n I 3n+1	n=2 96	0.61	1.48	2.04	2.17 96 opt
				n=3 96	0.63	1.61	2.26	1.99 96 opt
	(BA) 2 (EA) 2 Pb 3 I 10	ref. 97	0.99	1.94	2.54	2.12 97 opt
	(PMA) 2 PbI 4			ref. 98	1.25	2.16	2.81	2.19 99 opt
	BA 2 MA n-1 Sn n I 3n+1	n=1 100 n=3 100	0.82 0.12	1.83 1.07	2.37 1.59	1.83 100 opt 1.50 100 opt
	BA 2 MA n-1 Pb n Br 3n+1 PEA 2 MA n-1 Pb n Br 3n+1 (BEA) 2 PbBr 4	n=1 101 n=1 103 ref. 106	1.61 1.72 1.78	2.70 2.95 3.01	3.37 3.69 3.76	3.34 102 el 2.91-2.95 104,105 opt , 3.26 105 el -
	(EA) 4 Pb 3 Br 10		ref. 107	1.29	2.64	3.42	2.75 107 opt
	(EA) 4 Pb 3 Cl 10	l		ref. 107	1.77	3.39	4.24	3.45 107 opt
	(PMA) 2 PbCl 4	m	ref. 108	2.30	3.57	4.50	3.65 108 opt
	MAE (eV)				1.42	0.42	0.29
	MAPE (%)				61.7	18.9	12.0

(C(NH 2 ) 3 )(CH 3 NH 3 ) n PbnI 3n+1

TABLE V :

 V Average in-plane electron (m e ) and hole (m h ) effective masses in m 0 unit of 2D MHP calculated with PBE, DFT-1/2 and xs-DFT-1/2. µ = m e m h / (m e + m h ) is the reduced effective mass and m h = |m v |.

		Layer thickness	PBE	DFT-1/2	xs-DFT-1/2	Exp
			m v m e µ	m v m e µ	m v m e µ	µ
		n=1	-0.24 0.17 0.10 -0.43 0.26 0.16 -0.46 0.31 0.18 0.22 83a
		n=2	-0.22 0.13 0.08 -0.42 0.22 0.14 -0.45 0.26 0.17 0.22 83a
	BA 2 MA n-1 Pb n I 3n+1	n=3	-0.14 0.12 0.06 -0.26 0.22 0.12 -0.31 0.26 0.14 0.20 83a
		n=4	-0.13 0.11 0.06 -0.26 0.21 0.12 -0.31 0.27 0.14 0.20 83a
		n=5	-0.14 0.13 0.07 -0.25 0.23 0.12 -0.29 0.27 0.14 0.19 83a
	(C m H 2m+1 NH 3 ) 2 PbI 4	m=6	-0.25 0.17 0.10 -0.45 0.26 0.16 -0.48 0.31 0.19 0.18 88b
	PEA 2 MA n-1 Pb n I 3n+1	n=1 n=3	-0.17 0.12 0.07 -0.30 0.21 0.12 -0.38 0.25 0.15 -0.20 0.18 0.10 -0.34 0.30 0.16 -0.37 0.35 0.18
		n=1	-0.23 0.14 0.09 -0.77 0.22 0.17 -0.40 0.29 0.17
	(3AMP)MA n-1 Pb n I 3n+1	n=2	-0.14 0.13 0.07 -0.26 0.22 0.12 -0.29 0.26 0.14
		n=3	-0.15 0.11 0.06 -0.31 0.21 0.12 -0.33 0.26 0.15
		n=1	-0.54 0.19 0.14 -0.90 0.28 0.21 -0.76 0.31 0.22
	(4AMP)MA n-1 Pb n I 3n+1	n=2	-0.21 0.16 0.09 -0.35 0.26 0.15 -0.39 0.31 0.17
		n=3	-0.17 0.15 0.08 -0.31 0.26 0.14 -0.35 0.30 0.16
	(BA) 2 (EA) 2 Pb 3 I 10	-	-0.28 0.16 0.10 -0.48 0.27 0.17 -0.51 0.33 0.20
	(PMA) 2 PbI 4	-	-0.29 0.17 0.11 -0.54 0.26 0.18 -0.60 0.31 0.20
	BA 2 MA n-1 Sn n I 3n+1	n=1 n=3	-0.11 0.11 0.05 -0.22 0.19 0.10 -0.25 0.22 0.12 -0.06 0.07 0.03 -0.16 0.16 0.08 -0.20 0.20 0.10
	BA 2 MA n-1 Pb n Br 3n+1	n=1	-0.25 0.23 0.12 -0.41 0.36 0.19 -0.45 0.43 0.22
	PEA 2 MA n-1 Pb n Br 3n+1	n=1	-0.31 0.26 0.14 -0.49 0.41 0.22 -0.56 0.48 0.26
	(BEA) 2 PbBr 4	-	-0.37 0.29 0.16 -0.57 0.46 0.25 -0.63 0.54 0.29
	(EA) 4 Pb 3 Br 10	-	-0.42 0.28 0.17 -0.79 0.52 0.31 -0.64 0.55 0.30
	(EA)

m=10

-0.26 0.17 0.10 -0.46 0.26 0.17 -0.49 0.31 0.19

(C(NH 2 ) 3 )(CH 3 NH 3 ) n PbnI 3n+1 n=1 -0.

35 0.14 0.10 -0.65 0.23 0.17 -0.64 0.25 0.18 n=2 -0.18 0.14 0.08 -0.32 0.24 0.14 -0.36 0.29 0.16 n=3 -0.14 0.12 0.06 -0.23 0.22 0.11 -0.28 0.26 0.13 4 Pb 3 Cl 10 --0.33 0.81 0.23 -0.50 1.43 0.37 -0.56 1.03 0.36 a Measured at 4 K. b Measured at 5 K.

  Figure 6d compares the absolute valence energy levels of the 2D BA 2 PbI 4 42 and BA 2 PbBr 4

		116
	compounds computed at both DFT-1/2 and the refined xs-DFT-1/2 levels of theory. The
	corresponding slab is shown in Figure 6c. E abs v	predicted with xs-DFT-1/2 tends to be
	much deeper (by more than 1 eV) as compared to the normal DFT-1/2. Considering
	the experimental absolute valence energy level alignments reported by Silver et al. 102 on
	BA 2 PbI 4 (-5.8 eV) and BA 2 PbBr 4 (-6.5 eV), xs-DFT-1/2 appears to overestimate E abs v .
	Interestingly, the computed E abs v

TABLE VII :

 VII Comparison of the computational performance of different levels of theory using the relaxed pseudo-cubic CsPbI 3 bulk structure. Note that for PBE and DFT-1/2 calculations, inclusion of more empty bands is not necessary and has no influence on the computed eigenvalues. Number of bands Number of cores Band gap (eV) Elapsed wall time (s) Total computing time (s) a

	PBE (SIESTA) b	160	16	0.15	52.5	840
	PBE (VASP) c	80	16	0.22	56.4	902.4
	DFT-1/2 (SIESTA)	80	16	1.32	47.1	753.6
	DFT-1/2 (VASP)	80	16	1.30	65.7	1051.2
	HSE06(α=50%)	192	192	1.23	1274.6	2.4×10 5
	COHSEX@PBE d	960	960	1.36	3728.5	3.6×10 6
	G0W 0@PBE	960	960	1.02	9376.3	9.0×10 6

a We use: Total computing time = number of cores × elapsed wall time


A new methodology to analyse and optimize territorial compensations of solar radiation intermittency: a case study in Corsica Island (France)

Keywords: Solar radiation, Classification, Variability, Territorial compensations, Smoothing effect, Optimization

Territorial compensations of variability are a lever for smoothing the global photovoltaic production. They may allow a higher penetration rate of photovoltaic production in an energy mix, especially for non-interconnected islands. The aim of this paper is to analyse and optimise the smoothing effects of solar resource on a territory. In that view a novel typological classification-based methodology is proposed: intervals from highly intermittent episodes to ones with few fluctuations are discriminated. Then, the potential for compensation is demonstrated through the transformation from classes of high to low variability when summing solar radiation over several sites instead of one alone. The methodology is applied to Corsica. Solar radiation measured at 1 Hz on 8 sites which criss-crosses the island, when summed, demonstrates a reduction in the set of most intense variations compared to individual basis: they represent 5.5 % to 15.5 % of the studied period for the sites taken individually against only 1.7 % for the cumulated radiation. Optimisation of the territorial compensation strategy leads to a near-disappearance of intense variations, with this class representing only 0.13 % of the period. Finally, the optimized dispersion map of solar capacities results from this optimisation.

Introduction

Multicriteria selection of potential photovoltaic (PV) plants in Corsica revealed the challenges associated to the territorial concentration of grid-connected plants [START_REF] Haurant | Multicriteria selection aiding related to photovoltaic plants on farming fields on Corsica island: A real case study using the ELECTRE outranking framework[END_REF][START_REF] Oberti | Photovoltaic Plants Selection on an Insular Grid Using Multicriteria Outranking Tools: Application in Corsica Island (France)[END_REF]. Of primary concern is the abrupt loss in energy generation arising from localized heavy cloud cover. The issue is part of the broader challenge of resource intermittency, which has been extensively studied for the electricity infrastructures with a high penetration of photovoltaic technology [START_REF] Woyte | Voltage fluctuations on distribution level introduced by photovoltaic systems[END_REF][START_REF] Denholm | Evaluating the limits of solar photovoltaics (PV) in traditional electric power systems[END_REF][START_REF] Paatero | Effects of large-scale photovoltaic power integration on electricity distribution networks[END_REF][START_REF] Eltawil | Grid-connected photovoltaic power systems: Technical and potential problems-A review[END_REF]. The risk of grid instability has resulting in a limit being set on the amount of installed capacity that can be based on intermittent renewable sources. This limit can be highly constraining particuliarly for small grids such as non-interconnected islands. For instance, for the islands of La réunion and Corsica (France), this limit was set to 35 % of total electrical power by decrees [START_REF]Décret n°2015-1697 du 18 décembre 2015 relatif à la programmation pluriannuelle de l'énergie de Corse[END_REF][START_REF]Décret n°2017-530 du 12 avril 2017 relatif à la programmation pluriannuelle de l'énergie de la réunion[END_REF].

Today this regulatory constraint is a significant limiting factor for the transition to renewable energy technologies such as photovoltaic facilities, which may only progress if the technical feasibility of safely exceeding the limit can be proven. Small grids in non-interconnected zones such as on islands are especially impacted by this limit. In this context, the intermittence of the solar resource, its territorial variations and its predictability stand as major issues for the energy utility operator. A good knowledge of the solar resource variation over a territory would enable territorial compensation strategies to smooth photovoltaic generation as suggested by [START_REF] Perez | Spatial and temporal variability of solar energy[END_REF] and thereby allow the limit on renewable integration to be increased.

The spatial variability of solar radiation is typically shown through the analysis of inter-sites Pearson correlation coefficients r presented as a function of the distances between the studied sites [START_REF] Gallegos | Spatial variability of the global solar radiation obtained by the solarimetric network in the Argentine Pampa Humeda[END_REF][START_REF] Gueymard | Assessment of spatial and temporal variability in the US solar resource from radiometric measurements and predictions from models using ground-based or satellite data[END_REF][START_REF] David | Spatial and Temporal Variability of PV Output in an Insular Grid: Case of Reunion Island[END_REF][START_REF] Glasbey | Spatio-temporal variability of solar energy across a region : statistical modelling approach[END_REF][START_REF] Perez | Spatial & Temporal Characteristics of Solar Radiation Variability[END_REF]. A more advanced treatment was introduced by [START_REF] Badosa | Scales of spatial and temporal variation of solar irradiance on Reunion tropical island[END_REF] and [START_REF] Hoff | Quantifying PV power Output Variability[END_REF] who defined the discorrelation distance as the distance for which r drops below 0.5. Other authors have proposed statistical indices with which to compare fluctuations of solar radiation at various sites. For instance, [START_REF] David | Spatial and Temporal Variability of PV Output in an Insular Grid: Case of Reunion Island[END_REF] studied the standard deviation of hourly clearsky k cls between 2 sites, obtained with satellite data. [START_REF] Gallegos | Spatial variability of the global solar radiation obtained by the solarimetric network in the Argentine Pampa Humeda[END_REF] discussed the solar radiation difference between two sites by comparing the difference to the mean daily value of solar radiation. [START_REF] Lave | Solar variability of four sites across the state of Colorado[END_REF] compared solar radiation power spectra of several sites, and focused particuliarly on the linear best fit line for time scales less than 1 h.

The spatial variability of solar radiation demonstrated though the discorrelation between sites induces shifts in temporal fluctuations from a site to other.

Thus for a given period, statistical compensations of variations can be expected on a territory. These spatial compensations are exposed and analysed in various studies. [START_REF] Lave | Solar variability of four sites across the state of Colorado[END_REF] examined the average smoothing effect of solar radiation received shown to exist for a small geographic site in Japan by [START_REF] Otani | Study on areal solar irradiance for analyzing areally-totalized PV systems[END_REF]. In that work, the mean value of a network of 9 pyranometers covering a 4 km by 4 km area was found to fluctuate to a significantly lesser degree than the individual readings of each sensor.

Spatial compensations on solar radiation smooth the PV generation that allows a better integration in the energy mix. At PV plants scale first, [START_REF] Marcos | Power output fluctuations in large scale pv plants: One year observations with one second resolution and a derived analytic model[END_REF] analysed the PV output of 6 plants of different sizes, from 0.2 ha to 52 ha, scattered over a 1000 km 2 area in the south of Navarra (Spain). Smoothing of power fluctuations has demonstrated for large PV plants, and this effect is particularly noticeable for small time steps. At larger scale, [START_REF] Urquhart | Optimizing Geographic Allotment of Photovoltaic Capacity in a Distributed Generation Setting[END_REF] performed a multi-objective optimization, maximizing energy and minimizing variability, to allocate 2 MW PV among four sites, with differing amounts at each site, on the island of Lanai, Hawaii. For the optimal allotment, the variability in the form of ramp rates, were reduced by 50 %, whereas energy was reduced by only 3 % over the maximum production case. [START_REF] Jewell | Limits on cloud-induced fluctuation in photovoltaic generation[END_REF] have shown that the dispersal of PV installations over a wide territory can enable a higher degree of penetration of PV systems without jeopardizing grid stability, from 1.3 % for a fully centralized facility to 37 % for a distributed system covering 1000 km 2 .

Also, [START_REF] Collins | Systems-based modeling of generation variability under alternate geographic configurations of photovoltaic (PV) installations in Virginia[END_REF] studied how geographic dispersion of PV installations may minimize electricity supply uncertainty in Virginia (USA). They showed the probabilities of generation shortfalls was lowered for decentralized configurations, since the variability of the production are less variable.

Thus, the smoothing effects on solar radiation over a territory, when it is demonstrated, is a lever to limit the intermittence of the global photovoltaic production. Consequently, the territorial compensations of solar radiation offer a mean to increase the deployment of renewable energy technologies. These compensations are already analysed with statistical tools at different scales and territories in the literature. However no method has been proposed neither to quantified the smoothing effect of the solar radiation, nor to build a roadmap for an optimal dispersion of PV installation over a wide territory.

This paper addresses these two issues through a novel methodology for the The paper is structured with four sections: the overall methodology, based on typological classification and optimisation of territorial compensations of variation, is introduced in section 2. A case study comprising several different sites is described in section 3. Finally, results of the metholodogy applied to the case study are presented in section 4. Solar ressource variations of individual sites, as well as considering them altogether are qualified and classified. Finally, the presence of territorial compensations is demonstrated and then these are optimised for the dataset.

Methodology

Classification methodology

If source intermittency is well understood, periods of deficient solar energy can be mitigated by managing the various energy generation plants at the disposal of the utility operator. As a first step towards developing effective strategies to overcome intermittency, it is important to develop a means to characterise the feature in a rigorous manner using data with a sampling time short enough to capture rapid fluctuations in radiation intensity.

Temporal fluctuation of solar resource is usually characterized by the analysis of statistics on timeseries of solar radiation or photovoltaic generation data.

[12] calculated standard deviation of PV production measured at various times steps. [START_REF] Badosa | Scales of spatial and temporal variation of solar irradiance on Reunion tropical island[END_REF] used boxplots of k cls to describe distributions of this index via their median, 1 st and 3 rd quartiles. [START_REF] Glasbey | Spatio-temporal variability of solar energy across a region : statistical modelling approach[END_REF] applied autocorrelation coefficients presented as a function of different temporal lags. Moreover, ramp rates are also analysed through their probability densities, cumulative probabilities functions [START_REF] Lave | Solar variability of four sites across the state of Colorado[END_REF], power spectral analysis [START_REF] Marcos | Simulating the variability of dispersed large PV plants: Simulating the variability of dispersed large PV plants[END_REF][START_REF] Lave | Solar variability of four sites across the state of Colorado[END_REF] and counting of ramp rates occurrences over a given threshold [START_REF] Marcos | Power output fluctuations in large scale pv plants: One year observations with one second resolution and a derived analytic model[END_REF].

Irradiance measurements for multiple sites spanning a sufficient period with a high sampling frequency leads to a significant volume of data to analyse. Rather than conventional statistical tools, typological classifications can prove advantageous. Within this framework, a classification methodology was developed to identify and describe the nature of variations in solar radiation for individual sites in order to demonstrate the spatial compensations phenomenon.

The typological classification of solar radiation per day or instantaneously has been previously studied for satellite image processing [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF][START_REF] Mefti | Satellite approach based on cloud cover classification: Estimation of hourly global solar radiation from meteosat images[END_REF], the decomposition of total radiation into its direct and diffuse components [START_REF] Notton | Calculation on an hourly basis of solar diffuse irradiations from global data for horizontal surfaces in Ajaccio[END_REF] or the calculation of radiation incident on an inclined plane. Typically, within such studies, periods are classified into intervals of clear sky, partial cloud or overcast conditions, so that different models can be applied to each scenario. The segregation of a data into two or three intervals using instantaneous indicators of cloud cover or sky clarity provides a first means of classification. Thanks to its simplicity, this approach has been widely implemented in previous works [START_REF] De Miguel | Diffuse solar irradiation model evaluation in the North Mediterranean Belt area[END_REF][START_REF] Li | An analysis of climatic parameters and sky condition classification[END_REF][START_REF] Notton | Calculation on an hourly basis of solar diffuse irradiations from global data for horizontal surfaces in Ajaccio[END_REF][START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF][START_REF] Mefti | Satellite approach based on cloud cover classification: Estimation of hourly global solar radiation from meteosat images[END_REF]. However, since it evokes a static measure of daily cloudy conditions, the method is incapable of retaining dynamic information such as the degree of variability of cloud cover. Hence, a highly variable day can be attributed the same indicator of mean clarity as an overcast day.

Other typological methods allow the dynamics of passing clouds to be taken into account: the categorization by the mean clarity index k t was complemented by mathematical quantities to quantify the roughness of the signal such as the fractal dimension [START_REF] Maafi | Preliminary results of the fractal classification of daily solar irradiances[END_REF][START_REF] Harrouni | Classification of daily solar irradiation by fractional analysis of 10-min-means of solar irradiance[END_REF] or the Dirichlet coefficients [START_REF] Soubdhan | Classification of daily solar radiation distributions using a mixture of Dirichlet distributions[END_REF] or wavelet decomposition coefficients [START_REF] Woyte | Fluctuations in instantaneous clearness index: Analysis and statistics[END_REF]. [START_REF] Jeanty | Clustering daily solar radiation from Reunion Island using data analysis methods[END_REF] classified the daily distribution of direct fraction k b time series at reduced time-step. The methodology consists in combining three data mining methods: principal component analysis (PCA), K-means clustering methods and Ward hierarchical method to determine five typical days at La Réunion Island. [START_REF] Muselli | Classification of typical meteorological days from global irradiation records and comparison between two Mediterranean coastal sites in Corsica Island[END_REF] used a set of discrimination parameters derived from hourly clarity index k t to define three types of day in Corsica, applying Ward hierarchical method. The dynamics of solar radiation were taken into account through an integral of the square of the second derivative of k t . Similarly, [START_REF] Badosa | Scales of spatial and temporal variation of solar irradiance on Reunion tropical island[END_REF] identified five different solar radiation regimes at La Réunion, using k-means clustering analysis. They considered the irradiance level with the daily clear sky index and the variability through two parameters:

1. a Total Relative Change (TRC) of the clear sky index between 10 h and In this respect, these methods are inadequate for the electricity utility operator who must be able to forecast the intra-day conditions for energy generation.

To break away from this limitation, days may instead be considered as successions of sunlight regimes to be identified by a suitable classification algorithm.

The typological classification method of solar radiation variability by interval, mentioned in [START_REF] Haurant | Intermittence forecasting of the solar resource in Corsica[END_REF] is presented in a comprehensive manner in this paper.

Variation scale

Regimes in solar radiation distinguish themselves in particular by the form of the variations in irradiance that are observed, which can be characterized by two criteria: their duration and amplitude. The objective of this work is to develop a classification method of the regimes according to these properties.

The variation scale refers to the duration of time intervals for which the variation of solar radiation is evaluated. It is duration between two consecutive instants for which the relative variation δG is calculated, ignoring the variations present in the interval, following the relation:

δG(t) = G(t) -G(t -δt) G(t -δt) (1)

Categories of variations

Different solar radiation variation dynamics are discriminated by the variation scales. In effect, variations that are too slow to be perceived for a given δt can be revealed for larger time intervals (Fig. 1).

For a given day, if the data points with a variation greater than 10 % are selected, one observes that four discontinuous and relatively short groups are revealed for a variation scale of δt = 1 s, corresponding to times when the solar radiation fluctuates considerably. A large set of points is identified for a variation scale of δt = 300 s, and indeed it is observed that low frequency variations occurring between 3:00 pm -5:00 pm are particularly considered.

Finally, nearly all points in the interval 12:45 -5:20 pm are recovered when δt = 1200 s . The scales of variation therefore render visible and quantifiable the intensity of the variability following the rule according to which, the smaller the δt beyond which the variations become visible, the stronger and tighter the variations.

It is therefore possible to construct a typology on the basis of scales of variations according to the following procedure: considering the scales retained, calculate and compare the variations for a threshold in variability. The variations superior to this threshold are associated to the equivalent class, which can be semantically defined as a class of "variations visible at such a scale". Two distinct scales δt 1 and δt 2 give rise to three classifications, including the class of negligible variations. It is useful to choose scales that are separated by several orders of magnitude, in order to categorize qualitatively different variation profiles.

Delimitation of intervals

For most cases, even within a highly active cloudy regime, variations of solar radiation are found to only occasionally exceed the fixed threshold. The to regroup the data is therefore needed in order to limit the number of entities, whist respecting the similarity of variation profiles that have been regrouped.

The two filters have been introduced following this rational. The first threshold represented the maximum duration between two variations so that they are regrouped into the same interval. Labelled b intra , it is the intra-interval threshold. The second threshold compares to a minimum duration of intervals: it is the inter-interval threshold b inter . Using these two thresholds, the following four-steps sequential procedure (Fig. 2)

1. Variable instants t var are identified as instants for which the variations of solar radiation exceed the threshold of variability b var (Step 0 & 1):

t var = {t s.t. δG(t) ≥ b var } (2)
I var,i = I var,i ∪ {t var,i+1 } if t var,i+1 -max t∈Ivar,i (t) ≤ b intra (3)
3. Any intervals that are too short are integrated into adjacent intervals in the inter-interval regrouping step (step 3) :

I var,i+1 = I var,i ∪ I var,i+1 or I var,i-1 = I var,i ∪ I var,i-1 if max t∈Ivar,i (t) -min t∈Ivar,i (t) ≤ b inter (4)
In the case of adjacent intervals being of different classes, the inter-interval regrouping proceeds depending on its class: if the given interval is of the maximal class, it will be concatenated to the nearest neighbouring inferior class, otherwise it will be added automatically to the nearest superior class interval. These conditions are valid because, from the perspective of the utility operator, it is more convenient to choose the less favourable case in order to limit risks. However, in the first case where up-classing is impossible, the marginality or short duration of fluctuations are assumed to dominate their intensity and the interval is downgraded.

Optimisation of spatial compensation of solar radiation variations

The phenomenon of spatial compensation of solar radiation variations represents an original alternative to smoothing photovoltaic generation without the need for additional technologies. To this end, optimization of the phenomenon could serve to attenuate variations and provide a relatively stable photovoltaic energy supply over the whole territory.

The optimization problem shall be to find the set of coefficients

(β Si ) i∈[1,N ] , with ∀i ∈ [1, N ], β Si ≥ 0 and
G tot (t) = N i=1 β Si .G Si (t) (5)
.

So, objective function

F (β Si ) i∈[1,N ] is : F (β Si ) i∈[1,N ] = min (P t,Gtot (C M )) (6)
The coefficients corresponding to the minimum of the cost function indicate the optimal deployment of photovoltaic capacity at the sites included in the study for a minimum of C j variations and hence the smoothest possible generation curve.

In order to determine the optimal set of coefficients β Si an efficient optimisation tool is needed to limit the number of necessary calculations. Ultimately, all combinations (β Si ) i∈ [1,N ] should be tested to find the optimal set: for a coefficient precision of 1×10 -1 , this would entail 9 N tests. Rather than pursing this tedious approach, a genetic algorithm (GA) [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF][START_REF] Goldberg | Genetic algorithms in search, Optimization, and Machine-Learning[END_REF] was employed to accelerate the search.

Case study

The methodology developped in section 2 is applied to Corsica Island as it presents characteristics of an good illustrative case study. Indeed, the island can be considered as a non-interconnected zone, since the importation of electricity from the continental grid represents only a fifth of the total grid capacity. As a consequence, the regulatory limit (35 %) of renewable integration in Corsica was reached in 2012 when all existing PV plants were connected to the local grid [START_REF] Edf | Bilan prévisionnel de l'equilibre offre/demande d'électricité corse[END_REF]. The objective of the present work is to highlight and optimize the territorial compensation of solar resource variation as an effective strategy in order to exceed the regulatory integration limit for intermittent renewables without impacting on the grid. To this end, the spatial and temporal behaviour of solar radiation in Corsica was analysed, using ground-based measurements at a small number of weather stations.

Geography of the studied territory

Corsica is a Mediterranean island lying in the golf of Genova. From north to south the landmass spans 183 km, from latitude N 41°20 ′ 2 ′′ to N 43°1 ′ 31 ′′ .

The island is 83. unstable with a high degree of precipitation (about 700 mm yr -1 ). The geographical distribution of the solar potential across the island is therefore particularly heterogeneous and microclimates are present in the region, correlated with elevation as illustrated by the map of daily average radiation in 2005 obtained with processed Meteosat8 satellite images [START_REF] Haurant | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF].

The geographical dispersion of the studied sites is essential, in order to ensure a representative and coherent dataset for solar radiation across the territory.

The positions of the weather stations were therefore carefully chosen taking into account the microclimate zones on the island. The north-eastern and southwestern parts of the island do not experience variations in solar radiation at the same time. A line of separation exists, running north-west to south-east due to the orientation of the mountain range in the centre of the island. Coastal areas of the north-east and south-west, which represent equivalent microclimates, do not however experience the same weather. It is therefore convenient to install weather stations in both areas. The measurement sites were thus deployed at EDF net interconnection sites, taking into account the above mentioned constraints (Fig. 3), with altitude and geographical coordinates detailed in table 1.

Dataset

Instrumentation

The apparatus used, SP Lite 2 Kipp & Zonen pyranometers, are equipped with a silicon photodiode of sensitivity 60 to 100 µV W -1 m -2 . Their rapid response time to 95 % signal strength within 1 s allowed for the second-wise measurements undertaken in the current study. This time step was chosen in order to allow a thorough study of short term variations in solar radiation.

Measurement period

The data span a period of 18 months, starting the 24 th May 2010 for the stations at Corte and Ajaccio, and a few days later for the other installations. Data acquisition was stopped on the 2 nd October 2011. Over the period, occasional service interruptions lasting from a few hours to a few days were experienced.

Consequently,

• Some days contain several data gaps, rendering them unusable for the analysis. Those days with more than 10 % missing or unusable data were eliminated from the analysis

• The number and distribution of useable days of data in the dataset vary between sites (Table 1)

Applying the requirement that data from all stations were available (necessary for a comparative analysis), a subset of 152 days was kept for the present analysis. This constraint plus the heterogeneous distribution of the useable dates across the seasons represent serious limiting factors for this type of study.

Indeed, the present dataset comprised 62 days in winter, 63 in spring, 21 in summer and only 6 days in autumn. The set of overlapping dates therefore con-tains a seasonal imbalance, which introduces a statistical bias when comparing sites.

Around the hours of sunrise and sunset the effect of nearby obstructions and airmass result in low radiation intensities. Moreover, these phenomena are deterministic whereas the study is focused on stochastic variation due to cloudy cover. It was therefore necessary to define a so-called energy interval, representing 80 % of the solar energy of the day (at clear sky condition), symmetrical about solar noon and excluding the start and end of each day [START_REF] Haurant | De la sélection multicritère de parcs photovoltaïques à la cartographie et l'étude de l'intermittence de la ressource : modélisations appliquées à l'île de Corse[END_REF].

Parametrisation of the classification method

The development of the classification method as well as the definition of intervals of variability and of regimes of solar radiation require the determination of four parameters: The variability threshold was defined by the french grid operator in Corsica, EDF, which judges a variation to be important when it exceeds 10 % of the radiation received during 5 minutes. It was therefore decided to generalize this 10 % threshold for all of the variability scales to be considered.

The scales of variability δt serve to discriminate variations with different features. It suffices to adopt variation scales that are separated by orders of magnitude in order to discriminate different variation dynamics. Those of order a few seconds may be used to identify deep and narrow variations whereas those of a few minutes, as well as intense variations, recover slow and shallower variations.

Consider N tvar the number of instants t var identified, for a given variability scale δt, by the first step of the classification procedure (Fig. 2). Let's define all the study period. The light grey area despicts the range between the 10 th and the 90 th percentiles, dark grey despicts range between the 25 th and the 75 th percentiles and dashed line is the median P t (t var ) = N tvar /N ttot * 100 the ratio of number of variable instants to the total sample of the measurement period N ttot . P t (t var ) can be seen as the cumulative time variable on the whole sample.

On figure 4, we can observe that the ranges between 25 th and the 75 th of P t (t var ) and its derivative drawn a thick envelope curves. It shows that whatever the site, P t (t var ) has the same shape: P t (t var ) increases with the scale of variations δt considered. More precisely, P t (t var ) rises rapidly up until an inflexion of the envelop, for δt = 300 s. Very fewer variable instants are found for δt = 1 s than for δt = 300 s. These points reflect different types of variations to those observed for δt = 300 s : large variation scales can be used to identify slow variations in addition to deep and fast found with small scales. Then, we can see P t (t var ) increases slowly for δt ≥ 300 s and δt ≤ 1200 s: nearly the same instants are considered variable when δt is set in this intervale of values. Finaly, more variations identified at δt > 1200 s. These additionalvariations are not directly associated to meteorological instabilities, but rather may arise from the variation in solar radiation under clear sky conditions.

From this analysis of figure 4, it appears appropriate to undertake an analysis considering two variation scales, δt 1 = 1 s and δt 2 = 300 s, used to define three classes:

• Class 0 (C 0 ): stable conditions, with no visible variations for δt = 1 s or δt = 300 s;

• Class 1 (C 1 ): slow and shallow variations perceptible only for δt = 300 s;

• Class 2 (C 2 ): narrow and deep variations present for both scales, δt = 1 s and δt = 300 s.

The roles of b intra and b inter are equivalent: b intra can be interpreted as a maximal time interval between two variations, and therefore the minimum size of stable intervals; whereas the b inter designates the minimum size of variable intervals. These thresholds therefore represent the minimal sizes of intervals, stable or otherwise, and hence it is appropriate to consider identical values for the two parameters. Overall, the smaller the thresholds, the larger the degree of segmentation of days and the smaller the intervals hat will be obtained (Fig. 4). With this in mind, and in the absence of any objective criterion for the determination of these thresholds, they should rather be set be according to the specific application under study, the subsequent analysis, and the desired temporal resolution needed for that application.

Results and discussion

Characterisation of fluctuations observed at each site

The method introduced in section 2 was implemented to characterise the fluctuations observed at the 8 sites, considering a variability threshold of b var = 10 %, and two scales of variations δt 1 = 1 s and δt 2 = 300 s to allow the evaluation of the three classes defined above. The thresholds intra-and interregrouping were set to b intra = b inter = 900 s, affording a good compromise between temporal resolution for intervals and the number of intervals to be defined. This comparative assessment was undertaken using the 152 days of overlapping data.

The results of this comparative analysis of intervals are shown in figure 6, where each class C j is presented with its associated weight in terms of 1. relative cumulated durations P t (C j ) of intervals, obtained by evaluating the number of measurements associated to each class relative to the sample size.

2. relative number of intervals P N (C j ), obtained by counting the amount of interval of each class normalized by the total amount of intervals. The subset of data corresponding to daylight hours was retained for the comparative analysis. Pair-wise linear correlation coefficients were evaluated for the 152 days of overlapping data for each weather stations and the seven other sites. As can be seen in figure 7, the correlation factor between sites was found to not exceed 0.8 for the available data, indicating that the sites are not strongly correlated. The data exhibit phase differences between the variations observed from one site to the next, which therefore may be exploited to smooth the effects of intermittence and hence take advantage of the spatial compensations of the solar resource.

P • (C j ) [%] S 1 S 2 S 3 S 4 S 5 S 6 S 7 S 8 Site Pt (C1) PN (C1) Pt (C2) PN (C2) Pt (C2) PN (C2)
With the aid of the present data analysis, it has been possible to draw certain observations of the spatial compensation of the solar resource that are nevertheless indicative of the presence of the phenomenon for Corsica.

• Distinct variations in meteorological conditions recorded as a function of geographical position;

• correlation coefficients not exceeding 0.8 between sites (Fig. 7), arising from phase shifts in cloud phenomena; are the least common, covering only 5 -15 % of the measurement time (Fig. 6). When the data are cumulated over all sites, class C1 variations increase at the expense of the other classes (Fig. 9).

P t (C j ) [%] C 0 C 1 C 2

Variation Class

All sites Cumuled

Optimisation results

A simple GA was implemented to determine the coefficients (β Si ) i∈ [START_REF] Haurant | Multicriteria selection aiding related to photovoltaic plants on farming fields on Corsica island: A real case study using the ELECTRE outranking framework[END_REF][START_REF]Décret n°2017-530 du 12 avril 2017 relatif à la programmation pluriannuelle de l'énergie de la réunion[END_REF] for a range of tests, which converged after 70 to 75 generations of 20 individuals.

The obtained coefficients are summarized in table 2. They are the result of 8 applications of the GA on the scale of the study, the calculations being reiterated in order to test the sensitivity of the coefficients and thus exclude the hypothesis of a given solution corresponding to a local minimum. Operationnally, the β i coefficients can be translated as rates of the total PV capacity installed on the territory to be implemented at around each site. Thus, a map of the optimal dispersion of PV installations is shown in figure 10. %. Reciprocally the most unstable site, Bonifacio, was designated the smallest coefficient of 2.81 ± 0.47%. The optimization process had two positive outcomes (Fig. 11):

• The near total suppression of C 2 variations, which after optimisation account for only 0.13% of the measurement period of 6600 s;

• A reduction in the growth in frequency of C 1 variations. The occurrence of this class of variation was limited to a level comparable to the site with the highest frequency of C 1 variations. The ratios of points considered unstable relative to the number of data, for the sample of 152 days, are for any interval in the range dt ∈ [1; 1800] less for the spatially compensated combinations than for the individual sites (Fig. 12(a)).

The optimization improves on the simple sum of irradiances by reducing This is explained by the fact that the optimization processes is focused more towards the larger variations because the criterion to minimize is the number of variations observable at 1 s. In general, the smoothing afforded for this type of variation can only be partial, hence strong variations are downgraded to moderate ones, which are observable for larger scales and thus the ratio for these scales increases.

Conclusion

We have presented in this paper a new typological classification method for solar radiation, in which the principle innovation resides in the decomposition of days into intervals of different classes according the regime of solar radiation and its variability.

This classification method has been applied to the meteorological characterization of sites across the territory of Corsica, and in order to demonstrate the phenomenon of spatial compensation through an understanding of variations: the collection of events for which the variations were sufficiently important to be classified "strong variations" is limited when radiation is integrated across several sites. Therefore, the observable fluctuations at the scale of 1 s, classed C 2 , represent between 5.5 and 15.5 % of the measured time for individual sites whereas they represent 1.7 % for the cumulated intensity.

A genetic algorithm was subsequently implemented in order to optimise the compensation afforded by spatially distributed radiation, and was found able to practically eliminate intervals of strong variation. Variations of class C 2 represented only 0.13 % after optimization.

However, even once optimised, the compensation is not complete for the set of sites under consideration, and some variations were found to persist.

Therefore in the case of a high degree of penetration of intermittent renewable energy resources in Corsica, the utility operator will always need to manage the

  by 4 sites separated by several tens of kilometres (20 km to 180 km), at different time scales. The ramp rate analysis of the average revealed a significant decrease in their mean and maximum magnitude, standard deviation, and kurtosis compared to each site individually. Also, probability distribution functions and cumulative distribution functions of ramp rates indicated that the average of all four sites is less likely to have large fluctuations than each of the other sites individually. Similiarly, [15] compared boxplots of 1 Hz measurements at 7 isolated sites on the island of La Réunion to their aggregated radiation. The averaged radiation variability showed narrower distributions than individual stations, with less hour-to-hour variability. The spatial smoothing of the solar potential was

  analysis and optimization of the territorial compensations of radiation variability. The methodology is based on the classification of solar radiation variability and the quantification of durations of intermittent classes. The optimization consists then in finding the dispersion of PV capacity that minimizes the duration of the most intermittent class. To treat these issues, the main novelties of the paper are the following : • A new typological classification methodology of solar radiation variations • The use of a typological classification for studying the temporal variability of solar resource over measurement sites • The application of the typological classification for analysing and quantifying the smoothing effects on solar radiation over a territory • The optimization of the smoothing effect for minimizing the intermittence of PV production

16 h 2 .

 162 a Random Absolute Relative Change (RARC) index that represents the sum of short-term variability, excluding the variability trend of the day. These recent methods, less trivial and more satisfactory represent nevertheless the limitation of classifying days by averaging over instantaneous indicators, despite the fact that even within a single day it is possible to observe periods of different weather conditions. Once again, the characteristics of each weather phenomenon becomes smoothed and attenuated and, consequently distinctly different days may become embroiled together. Therefore, the categories represent a large number of different radiation profiles. An additional shortcoming is their limitation to the classification of complete days rather than short intervals.

Figure 1 :

 1 Figure 1: Variations revealed by different values of variations scales : δt = 1 s in red, δt = 300 s in black and δt = 1200 s in green

Figure 2 :

 2 Figure 2: Diagram of the method : 0.&1. Identification of variations (bvar = 0.1 ) ; 2. Intrainterval grouping step (b intra = 900 s) and 3. inter-interval grouping step (b intra = 900 s ) (measured irradiance at Corte, 24 th May 2010)

N i=1 β

 i=1 Si = 1 that minimizes P t,Gtot (C M ) the relative cumulated duration of most important variations C M in the weight sum of global irradiance G tot over all the N sites :

Figure 3 :

 3 Figure 3: Geographic localisation of measurement stations S i , micro-climates and elevation contours in Corsica. The greyscale indicates daily irradiation intensity in 2005 over the island. Cyan, blue and red contours highlight the elevation levels at 530 m, 1060 m and 2120 m respectively.

1. A variability threshold b var beyond which the conditions are considered fluctuant 2 .

 2 The values of the variation scales δt used to discriminate the different regimes of solar radiation 3. The thresholds of intra-and inter-intervals regrouping, designated respectively b intra and b inter , which set the limits of the different intervals.

Figure 4 :

 4 Figure 4: Statistics of Pt (tvar) as a function of the scale of variability δt for the 8 sites on

Figure 5 Figure 5 :

 55 Figure 5 presents intervals delimited for different regrouping thresholds, which are observed to as significant an impact on the data as the configurations, when the thresholds are varied from b intra = b inter = 1 s to b intra = b inter = 300 s, b intra = b inter = 900 s and finally to b intra = b inter = 1800 s.As was previously expected, without regrouping (b intra = b inter = 1 s) a large number of small discontinuous intervals are observed (Fig.5(a)). The application of these thresholds serves to eliminate marginal variations. Thus, all unstable points classified as C 1 observed between 10:28 and 10:33 (Fig.5(a)) change to C 0 when regrouping is applied (Fig.5(b), 5(c) et 5(d)). The intervals generated by the instabilities during these times do not exceed the threshold Sinter, and are therefore declassified. These thresholds are responsible for a considerable restructuring of the variation intervals. Therefore,

Figure 6 :

 6 Figure 6: Histograms of relative numbers of intervals P N (C j ) and relative cumulated durations Pt (C j ) of variation class C j for the studied sites S i

Figure 7 :

 7 Figure 7: Variagramm of the irradiance data per paire of sites: standard variations between irradiances measured over paire of sites as a function of distances

•

  Different fluctuations observed by the classification method, indicating different climate and weather from one site to another. Also, we can observe that the degree of fluctuation observed for individual sites (Fig.8(a), 8(b), 8(c), 8(d), 8(e), 8(f), 8(g), 8(h)) is attenuated when the solar radiation is integrated over all sites (Fig.8(i)). For instance, some sites exhibit significant variations in solar radiation during the afternoon: such variations of class C 2 were observed in Lucciana (S 1 ), Calvi (S 2 ) and Bonifacio (S 8 ). These variations are attenuated when considering the average, and none of them are of class C 1 . These observation are consistent with earlier evidence.The distribution of solar radiation over the island is therefore favourable to concept of territorial smoothing of the solar resource, although it remains necessary to develop a more general proof for all days and also a method to quantify the level of compensation achieved. Conventional statistical methods such as histograms of variations are inadequate for the current study. Indeed, we have not found any suitable data representation tool to visualise and quantify the variations. The use of variation classes could serve to prove the existence of the phenomenon, by presenting the class distribution for individual sites and the cumulated signal.

Figure 8 :

 8 Figure 8: Irradiances and variations classes (bvar = 10 %, δt 1 = 1 s, δt 2 = 300 s, and b intra = b inter = 900 s) over the sites and average, on 13 th February 2011

Figure 9 :

 9 Figure 9: Relative cumulated durations of intervals Pt (C j ) of each class C j considering the individual sites in gray histograms (the errorbar depicts the minimal and maximal values all sites considered) and the accumulation of the global radiation over all the sites, illustrating territorial compensation in blue histograms

Figure 10 :

 10 Figure 10: Territorial redistribution of solar electricity generation for an optimal compensation of variations. The sizes of the squares denote the optmial rate of PV global capacity on the site zone and the thickness of the edges represent the uncertainties of the optimization

Figure 11 :

 11 Figure 11: Relative cumulated durations of intervals Pt (C j ) of each class C j considering the individual sites in gray histograms (the errorbar depicts the minimal and maximal values all sites considered) ; the accumulation of the global radiation over all the sites, illustrating territorial compensation in blue histograms and the compensations phenomenon optimization in red histograms

1 ≤Figure 12 :

 112 Figure 12: Proportion of unstable periods to the total data sample for the 8 sites; for the territorial variation compensations (blue) and the optimized case (red)

  5 km wide from east to west, covering longitudes E 8°32 ′ 30 ′′ to E 9°33 ′ 38 ′′ . With a surface area of 8722 km 2 , Corsica is the fourth largest island in the Mediterranean after Sardinia, Sicily and Cyprus. The island is classed as a "mountain in the sea", with relief above 600 m accounting for 39

% of its surface and an average elevation above sea level of 568 m. Some 120 peaks exceed 2000 m altitude (Rome and Giorgetti, 2007). The coastal regions enjoy a Mediterranean climate characterized by warm, dry summers and mild, humid winters, which translates into a higher solar resource potential than most of metropolitan France. In contrast, in the mountains the climate is strongly affected by the altitude. At the highest altitudes Corsica exhibits an alpine climate. Inland, summers are hot and sunny, but in winter conditions are more

Table 1 :

 1 Geographical coordinates of measures sites

	Code	Site	Latitude (φ)	Longitude (λ)	Altitude	N er days
	S1	Lucciana	42°39 ′ 49 ′′	9°25 ′ 28 ′′	60 m	456
	S2	Calvi	42°33 ′ 38 ′′	8°44 ′ 48 ′′	31 m	416
	S3	Corte	42°18 ′ 4 ′′	9°9 ′ 57 ′′	381 m	445
	S4	Ghisonaccia	42°3 ′ 54 ′′	9°22 ′ 14 ′′	65 m	451
	S5	Ajaccio	41°55 ′ 49 ′′	8°45 ′ 23 ′′	2 m	355
	S6	Sainte Lucie	41°41 ′ 59 ′′	9°20 ′ 12 ′′	66 m	405
	S7	Propriano	41°39 ′ 43 ′′	8°55 ′ 2 ′′	17 m	407
	S8	Bonifacio	41°22 ′ 17 ′′	9°12 ′ 10 ′′	46 m	450

Table 2 :

 2 β i coefficients for each site S i , with their uncertainties obtained with GA optimization ran 8 times

	Code	Sites	β S i [%]
	S1	Lucciana	26.04 ± 1.80 %
	S2	Calvi	11.71 ± 0.65 %
	S3	Corte	16.20 ± 2.47 %
	S4	Ghisonaccia	8.40 ± 0.73 %
	S5	Ajaccio	16.51 ± 2.01 %
	S6	Sainte Lucie	5.29 ± 0.33 %
	S7	Propriano	12.55 ± 0.57 %
	S8	Bonifacio	2.81 ± 0.47 %

ensure the stability of the grid supply and guaranty the electricity service.