Elis Marina Sales De Castro

Frédéric Giraud

Low-frequency

published or not. The documents may come


Hugh G A Burton

Pierre-François Loos

Rationale for the Extrapolation Procedure in Selected Configuration Interaction

Selected configuration interaction (SCI) methods have emerged as state-of-the-art methodologies for achieving high accuracy and generating benchmark reference data for ground and excited states in small molecular systems. However, their precision relies heavily on extrapolation procedures to produce a final estimate of the exact result. Using the structure of the exact electronic energy landscape, we provide a rationale for the common linear extrapolation of the variational energy as a function of the second-order perturbative correction. In particular, we demonstrate that the energy gap and the coupling between the so-called internal and external spaces are the key factors determining the rate at which the linear regime is reached. Starting from first principles, we also derive a new non-linear extrapolation formula that improves the post-processing of data generated from SCI methods and can be applied to both ground-and excited-state energies.

I. INTRODUCTION

Selected configuration interaction (SCI), [1][2][3][4] and related methods (such as density-matrix renormalisation group approaches [5][6][7] and others [8][9][10][11][12][13][14][15][16][17] ), have taken a prominent role in modern electronic structure theory. [18][19][20] Their primary purpose is to calculate reference correlation and excitation energies in small molecular systems, 19,[21][22][23][24][25][26][27][28] for which they have demonstrated a remarkable ability to yield highly accurate estimates of full configuration interaction (FCI) results. The numerous variations of SCI all perform a sparse exploration of the Hilbert space by selecting only the most energetically relevant determinants. This natural philosophy emerges from the observation that, among the incredibly large number of determinants in the FCI space, only a tiny fraction of them significantly contribute to the energy. Modern versions of SCI include CIPSI (CI using a Perturbative Selection made Iteratively) 3,18,21,[29][30][START_REF] Caffarel | Using cipsi nodes in diffusion monte carlo[END_REF][START_REF] Garniron | [END_REF][33][34][35][36] adaptive sampling CI (ASCI), [37][38][39][40] semistochastic heatbath CI (SHCI), 23,24,[41][42][43][44][45] iterative CI (iCI), [46][47][48][49][50] Monte Carlo CI (MCCI), 51,52 and FCI quantum Monte Carlo (FCIQMC). [53][54][55][56][57][58] The SCI wave function corresponds to a truncated CI expansion constructed from determinants in some internal (or model) space

I |Ψ var ⟩ = I∈I c I |I⟩ , (1)
with the associated variational energy E var = ⟨Ψ var | Ĥ|Ψ var ⟩, where we assume the normalisation of the variational wave function, i.e., ⟨Ψ var |Ψ var ⟩ = 1. The accuracy of |Ψ var ⟩ can be assessed using the second-order Epstein-Nesbet perturbation correction, computed using the determinants {α} that lie outside the model space a) Electronic mail: hgaburton@gmail.com b) Electronic mail: loos@irsamc.ups-tlse.fr (i.e. in the external space A) as

E PT2 = - α∈A | ⟨Ψ var | Ĥ|α⟩| 2 H αα -E var , (2)
where H αα = ⟨α| Ĥ|α⟩. The exact FCI wave function and energy are indicated by the limit E PT2 → 0 -. Despite the sparse exploration of the Hilbert space, these state-of-the-art methods still rely heavily on extrapolation procedures to produce final FCI estimates. 19,20,23,34 In particular, it is widely observed that E var becomes approximately proportional to E PT2 for small E PT2 , and thus a linear or quadratic extrapolation of E var for E PT2 → 0 is generally used to estimate the exact energy for an unconverged SCI calculation. 23 The precision and reliability of this post-processing extrapolation procedure are critical in order to produce meaningful estimates. However, to the best of our knowledge, no theoretical justification for the linear (or otherwise) relationship between E var and E PT2 has been proposed.

To illustrate this extrapolation procedure, Fig. 1 shows the evolution of the variational correlation energy of benzene as a function of E PT2 , computed in the cc-pVDZ basis and within the frozen-core approximation. These SCI calculations were performed with quantum package using the CIPSI algorithm 34 and the data are extracted from Ref. 35. The FCI estimate of the correlation energy (solid black line in Fig. 1) was estimated to be -862.890 mE h and was obtained by performing a five-point linear fit (dashed black line in Fig. 1) of the CIPSI data. This estimate carries an error of the order of 1 mE h and the fitting error was estimated to be 0.266 mE h . From Fig. 1, it is clear that, for sufficiently small E PT2 , the variational quantity is linear with respect to E PT2 . However, the SCI data deviate significantly from linearity for larger values of E PT2 , which we shall address in detail later on.

In this Communication, we provide a rationale to justify the linear extrapolation of the (zeroth-order) variational energy as a function of the second-order perturbative energy. We adopt a geometric approach that considers the variational wave function as a point on the exact electronic energy landscape, 59 allowing the second-order perturbative correction to be derived from the local gradient and curvature of this energy landscape. Moreover, we investigate a two-state model in which an analytic relationship between E var and E PT2 can be derived, leading to a novel parametrised non-linear formula that facilitates a more robust extrapolation procedure.

II. RATIONALE FOR THE LINEAR EXTRAPOLATION

The linear relationship between E var and E PT2 can be derived from first principles by considering the structure of the electronic energy landscape. While Ref. 59 describes this energy landscape perspective in detail, the salient points are summarised here. Any normalised wave function in the full N -dimensional Hilbert space

|Ψ⟩ = N I=1 v I |I⟩ (3)
can be represented by a vector v subject to the normalisation constraint v † • v = 1, which constrains the wave function to the surface of a hypersphere. The energy is given by the quadratic form

E = v † • H • v (4)
and exact eigenstates of the Hamiltonian correspond to stationary points of E constrained to the surface of the hypersphere, as illustrated in Fig. 2. At any point on the hypersphere, the tangent space T contains the vectors that are orthogonal to v, which can be collected as the columns of an N × (N -1) matrix v ⊥ . These tangent vectors correspond to the states |T ⟩ that are orthogonal to |Ψ⟩ and satisfy

1 = |Ψ⟩⟨Ψ| + T ∈T |T ⟩⟨T | , (5)
|Ψ 0 ⟩ |Ψ 1 ⟩ |Ψ 2 ⟩ Tangent space T FIG.
2. Sketch of the exact electronic energy landscape in a three-dimensional Hilbert space. Eigenstates correspond to stationary points constrained to the surface of the unit sphere. At any point, the tangent space T is spanned by the two vectors (red) that are orthogonal to the position vector.

where 1 is the identity operator. A constrained step s on this landscape is parametrised using a unitary transformation as

|Ψ(s)⟩ = exp T ∈T s T |T ⟩⟨Ψ| -|Ψ⟩⟨T | |Ψ⟩ . (6)
Assuming real wave functions, the components of the constrained energy gradient are then

g T = ∂E ∂s T s=0 = 2 ⟨T | Ĥ|Ψ⟩ , (7)
while the elements of the Hessian matrix of constrained second-derivatives become

Q T T ′ = ∂ 2 E ∂s T ∂s T ′ s=0 = 2 ⟨T | Ĥ -E|T ′ ⟩ . (8)
So far, we have only considered the structure of the electronic energy landscape for an arbitrary wave function in the full Hilbert space. For a SCI variational wave function, the only non-zero elements of v correspond to determinants included in the internal space, with coefficients c I , as defined in Eq. ( 1). The tangent vectors can then be split into two disjoint sets corresponding to the eigenstates within the internal space that are orthogonal to |Ψ var ⟩, denoted J = {|J⟩} J̸ =var , and the determinants in the external space, giving T = J ∪ A, as illustrated in Fig. 3. Since |Ψ var ⟩ is a CI solution within the internal space, such that ⟨J| Ĥ|Ψ var ⟩ = 0, the gradient given by Eq. ( 7) is only non-zero in the direction of the tangent vectors in A.

The local structure of the energy landscape around the variational wave function |Ψ var ⟩ is given by a second-order Taylor series expansion as Optimising this quadratic form with the Newton-Raphson step s = -Q -1 • g gives an estimate of the difference between the exact energy E exact and E var as

E(s) = E var + s † • g + 1 2 s † • Q • s. (9
∆E = E exact -E var ≈ - 1 2 g † • Q -1 • g, (10)
or, equivalently, ∆E ≈ -

αα ′ ∈A ⟨Ψ| Ĥ|α⟩ ⟨α|( Ĥ -E var 1) -1 |α ′ ⟩ ⟨α ′ | Ĥ|Ψ⟩ , (11)
where we have exploited the fact that the gradient is zero in the direction of the tangent vectors within J . The exact energy landscape is quadratic near an eigenstate, 59 and thus Eq. ( 10) will become the exact energy correction when |Ψ var ⟩ is sufficiently close to an eigenstate. Since the matrix elements ⟨α|( Ĥ -E var 1) -1 |α ′ ⟩ are expensive to compute, we can assume that the Hamiltonian is diagonally dominant in A and take the leading order approximation

Ĥ(0) = II ′ ∈I |I⟩ H II ′ ⟨I ′ | + α∈A |α⟩ H αα ⟨α| , (12)
where Ĥ(0) is the zeroth-order Hamiltonian within the Epstein-Nesbet partitioning. The energy correction then reduces to the second-order Epstein-Nesbet expression obtained in Eq. ( 2) to give

∆E ≈ - α∈A | ⟨Ψ var | Ĥ|α⟩| 2 H αα -E var = E PT2 . (13)
Assuming that |Ψ⟩ is sufficiently close to an eigenstate, such that the quadratic approximation to the exact energy landscape is valid, we can combine Eqs. ( 10), (11), and (13) to obtain E var ≈ E exact -E PT2 , thus rationalising the linear extrapolation of E var as a function of E PT2 .

III. INSIGHTS FROM A TWO-STATE MODEL

In practice, linear or quadratic extrapolation procedures only work using a limited number of points for well-converged calculations (see Fig. 1). As mentioned above, the variational energy generally appears to deviate away from linearity for larger values of E PT2 . The cause of these deviations can be studied using a two-state model that represents the separation of the internal and external spaces in a SCI calculation, and for which the relationship between E var and E PT2 can be analytically derived.

Our model contains individual states |Ψ I ⟩ and |Ψ A ⟩ representing wave functions in the internal and external spaces, respectively, with characteristic energies E I and E A . The Hamiltonian matrix in this basis is then

H = E I t t E A , ( 14
)
where t = ⟨Ψ I | Ĥ|Ψ A ⟩ represents the strength of the coupling between the internal and external spaces, and δE = E A -E I provides a measure of their energetic separation. The exact ground-state energy is

E exact = E I + δE 2 - δE 2 2 + t 2 . ( 15
)
The improvement of |Ψ var ⟩ during the course of a SCI calculation can be modelled by mixing |Ψ I ⟩ and |Ψ A ⟩ to give the parametrisation

|Ψ var (θ)⟩ = cos θ |Ψ I ⟩ + sin θ |Ψ A ⟩ , (16)
with 0 ≤ θ < 2π. The corresponding energy is

E var (θ) = E I + δE 2 (1 -cos 2θ) + t sin 2θ. ( 17
)
Following a Taylor series expansion, the second-order correction is

E PT2 (θ) = - 1 4
(2t cos 2θ + δE sin 2θ) 2 δE cos 2θ -2t sin 2θ .

By solving Eq. ( 18) for θ, we can invert these equations to express E var in terms of E PT2 as

E var = E I + δE 2 -E PT2 - δE 2 2 + t 2 + (E PT2 ) 2 .
(19) which naturally reduces to Eq. ( 15) for E PT2 = 0. This expression reveals that the more general form of E var as a function of E PT2 involves a square-root term that deviates away from linearity, and that this departure from the linear regime is directly related to the energetic separation (δE) and coupling strength (t) between the internal and external spaces. For (E PT2 ) 2 ≪ δE 2 2 + t 2 , we recover the linear behaviour E var ≈ E exact -E PT2 . In a real SCI calculation, we expect |t| ≪ δE 2 . Therefore, the larger the energy separation, the sooner the linear regime is reached, while a large coupling between I and A also leads more rapidly to the linear regime. These features are further demonstrated through the series expansion of E var at small E PT2 ,

E var = E exact -E PT2 - E 2 PT2 2 δE 2 2 + t 2 + O E 3 PT2 , ( 20
)
which shows that the quadratic behaviour is minimal when |δE| or |t| is large. This functional relationship can be illustrated by evaluating E var for various values of E PT2 in the limit of interest, as shown in Fig. 4 for E I = -1, δE = 1 and t = 1. As E PT2 gets larger, E var strongly deviates from linearity and bears a close similarity to previous SCI data. 19,21,[23][24][25][26][27]34 Evar (red markers) as a function of EPT2 for the two-state system with EI = -1, δE = 1 and t = 1. These data deviate from the linear approximation (black dashed line) due to the square-root term in Eq. ( 19).

For real systems, this two-state scenario can always be engineered by using a singular value decomposition of the gradient vector to transform the tangent vectors such that only one direction has a non-zero gradient. In other words, we can transform the external space A such that only a single state |Ψ A ⟩ = α∈A c α |α⟩ couples to |Ψ var ⟩ through the Hamiltonian. Furthermore, since the states within I are decoupled by solving the CI problem, this two-state model can be constructed for both ground and low-lying excited states, the only difference being the precise values of E I , δE, and t. While this transformation is not feasible in practice, it provides a conceptual link between real SCI data and the present two-state model.

IV. A NON-LINEAR EXTRAPOLATION FORMULA

The insights from our two-state model suggest that a non-linear functional form is more suitable for extrapolating SCI data. The separation of the model and perturbation space in SCI means that the relationship between E var and E PT2 is not as straightforward as the model system. However, as more determinants are added to the model space, the variational wave function follows a path towards the exact ground state that is likely to resemble Eq. (19). The concave form of Eq. ( 19) suggests that a linear extrapolation procedure will generally underestimate the exact correlation energy. Therefore, we propose a new non-linear extrapolation formula

E var (E PT2 ; a, b, c) = a+ |c| 2 -bE PT2 - c 2 2 + (bE PT2 ) 2 , (21)
where a, b, and c are fitting parameters. This expression corresponds to a form of quadratic approximant, which has been previously used in the resummation of divergent perturbation expansions. [60][61][START_REF] Goodson | Mathematical Physics in Theoretical Chemistry, Developments in Physical & Theoretical Chemistry[END_REF][START_REF] Mayer | [END_REF] Crucially, Eq. ( 21) reduces to a linear fit for |E PT2 | ≪ c 2b and can reproduce the observed non-linearity for larger E PT2 . The fitted value of a provides the estimate for the FCI result.

The variation of the FCI estimate of the correlation energy of benzene computed in the cc-pVDZ basis using the linear or non-linear extrapolation procedure is compared in Fig. 5. These data show that the non-linear formula can accurately fit the SCI data at larger E PT2 values than the linear procedure. In the left panel of Fig. 6, we study the influence of the number of points included in the linear or non-linear fits, each of them starting from the point associated with the smallest value of E PT2 . The error bars indicate the standard errors associated with the fitting procedure. Although the FCI estimates obtained from the non-linear formula [see Eq. ( 21)] have larger fitting errors for a small number of points, they quickly stabilise and remain relatively consistent compared to those obtained from the linear procedure, which rise much faster. Consequently, a larger number of points can, and should, be employed for the non-linear extrapolation formula, while the linear extrapolation procedure becomes systematically worse when more points are used.

In the right panel of Fig. 6, we consider a fixed-length window of consecutive points for the linear and non-linear fits, but we vary the index of the starting point, with the point at index 1 corresponding to the smallest value of E PT2 (i.e., the higher the index of the starting point, the larger the E PT2 correction). Using 5 and 8 points for the linear and non-linear fits, respectively, we again show that the non-linear procedure is slightly more stable than its linear counterpart with respect to the starting point, although both extrapolation schemes eventually underestimate the correlation energy for larger E PT2 values. These results indicate that, compared to the linear approach, the non-linear extrapolation procedure can provide a more accurate estimate of the FCI result for SCI calculations that are less well converged.

V. CONCLUDING REMARKS

In this Communication, we have proposed a theoretical rationale for the linear extrapolation procedure of the variational zeroth-order energy, E var , as a function of the second-order perturbative correction, E PT2 , commonly employed in SCI methods. Our derivation is based on connecting the SCI variational wave function to the properties of the underlying electronic energy landscape. The accuracy of the extrapolation of E var as E PT2 approaches 0 is critical for the accurate determination of the final FCI estimate.

Our investigations led us to the discovery of a novel non-linear extrapolation formula that more effectively captures the behaviour of E var for larger E PT2 values, thereby enhancing the robustness of extrapolations toward the FCI limit. Based on a two-state model, we derived the analytic form of this non-linear extrapolation formula and examined its mathematical properties. Specifically, we illustrated that the rate at which the linear regime is attained is primarily determined by the energetic gap and the coupling between the internal and external spaces. As a concrete example, we studied the ground-state correlation energy of benzene, which illustrates the versatility and reliability of this new extrapolation procedure.

We anticipate that this study will facilitate SCI calculations on larger Hilbert spaces, while providing confidence in the validity of extrapolation procedures. Furthermore, we expect that our two-state model and our framing of the SCI approach within the electronic energy landscape framework will allow other intriguing mathematical aspects of these methods to be explored in the future.

SUPPLEMENTARY MATERIAL

See the supplementary material for the SCI data of Fig. 1 and the raw data associated with Fig. 6.

Horizon 2020 research and innovation programme (Grant agreement No. 863481).

FIG. 1 .

 1 FIG.1. Frozen-core (variational) correlation energy of benzene as a function of EPT2 computed in the cc-pVDZ basis, as described in Ref.35.

)FIG. 3 .

 3 FIG.3. Sketch of the tangent space construction T = J ∪ A. The variational space I (dashed blue circle) is built from two configurations (blue vectors) and the external space A contains one configuration (red vector). At |Ψvar⟩, the tangent space is spanned by one tangent direction |J⟩, which is locally parallel to the variational space I, and one orthogonal tangent direction in the external space |α⟩.

  FIG.4. Evar (red markers) as a function of EPT2 for the two-state system with EI = -1, δE = 1 and t = 1. These data deviate from the linear approximation (black dashed line) due to the square-root term in Eq.(19).

FIG. 5 .

 5 FIG.5. Comparison of the linear and non-linear extrapolation procedure for the correlation energy of benzene computed in the cc-pVDZ basis, using the SCI data from Ref.35.

4

 4

  Evolution of the FCI estimate of the correlation energy of benzene computed in the cc-pVDZ basis. Left: A variable number of fitting points is considered in the linear (blue) or non-linear (green) fits, with each fit starting from the point associated with the smallest value of EPT2. Right: A fixed-length window of consecutive points is used for the linear (cyan) or non-linear (magenta) fits. The point of index 1 corresponds to the smallest value of EPT2, and a higher index for the starting point corresponds to a larger EPT2 correction. The error bars indicate the standard errors associated with the fitting procedure.

						-845
	-858
	-860					-850
	-862					-855
	-864
	-866					-860
	-868					-865
	6	8	10	12	14	16	2	4	6	8	10	12	14
	FIG. 6.

ACKNOWLEDGEMENTS

HGAB thanks Downing College, Cambridge for support through the Kim and Julianna Silverman Research Fellowship. PFL thanks financial support from the European Research Council (ERC) under the European Union's

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article and its supplementary material.


Guillaume Le Gall

Martin Thebault

Cyril Caliot

Julien Ramousse

PRINCIPAL COMPONENT ANALYSIS FOR THE CHARACTERISATION OF SPATIOTEMPORAL VARIATIONS OF THE SOLAR RESOURCE IN URBAN ENVIRONMENTS

Keywords: Solar radiations, Urban environment, Spatiotemporal variability, Model Order Reduction, Principal Component Analysis

Urban areas are serious candidates for the production of solar energy but their intrinsic complexity makes it challenging. The heterogeneity in the geometries and radiative properties of the different elements composing the urban fabric, specifically induces important spatiotemporal variations of the distribution of incident solar radiations. Besides, Principal Component Analysis (PCA) has been widely validated as an efficient tool to identify the principal behavioural features of a high-dimensional physical model. This paper proposes a novel approach to analyse and characterise the spatiotemporal variability of the solar resource within an urban context by means of PCA. A theoretical 100 × 100 m 2 asymmetric urban district made of nine cuboids with various heights is studied. The distribution of the incident field of irradiances is modelled via backward Monte-Carlo ray tracing over a full year on the facets of the central building under a clear sky, with a 15 min timestep and 1 m spatial resolution. PCA is subsequently applied to the simulated model to analyse its spatial and temporal variabilities. First results validate modal decomposition as a powerful technique for the analysis of the variability distribution, allowing the identification of the district areas subjected to important spatial and temporal variations of the solar resource. Characteristic scales are clearly represented by orders of decomposition. The contribution of surrounding geometries is also transcribed by particular spatial modes and similar influential variables are encountered across multiple evaluated surfaces but at different modal ranks.

INTRODUCTION

Sunlight is an abundant, inexhaustible and easily accessible natural resource for energy production. Over the last decades, significant efforts have been made to position solar energy as a tangible alternative to traditional non-renewable fuels [START_REF] García-Cuesta | A combination of supervised dimensionality reduction and learning methods to forecast solar radiation[END_REF]. If appropriate measures and support policies are being taken and rapidly implemented by governmental authorities, it could be able to meet over 30% of the global energy needs by 2060 [START_REF]Solar Energy Perspectives[END_REF]. Urban areas are serious candidates for the deployment of solar collectors, since a substantial part of the energy produced is being used directly in cities [START_REF] Freitas | Modelling solar potential in the urban environment: State-of-the-art review[END_REF].

However, the complexity of the urban environment makes the production of solar energy challenging. The presence of buildings, vegetation and other obstacles, and the diversity of their geometrical nature *Corresponding Author: guillaume.le-gall@univ-smb.fr (e.g., shape, size, orientation and inclination) and materials optical properties (e.g., specular and/or diffuse reflectivity, transmissivity) alter solar and infrared radiative transfers [START_REF] Caliot | Model of Spectral and Directional Radiative Transfer in Complex Urban Canopies with Participating Atmospheres[END_REF]. Multiple reflections of solar radiations and complex dynamic overshadowing effects are taking place between structural elements [START_REF] Freitas | Modelling solar potential in the urban environment: State-of-the-art review[END_REF]. Building façades are specifically of important influence here. This participates in the creation of a specific local climate in interaction with the urban infrastructures, contributing to the Urban Heat Island (UHI) effect [START_REF] Johari | Urban building energy modeling: State of the art and future prospects[END_REF]. Another important agent of solar radiations regulation is the atmosphere itself. Light from the sun passing through the atmosphere is subjected to multiple radiative events (absorption, emission, scattering and reflection) before reaching the ground [START_REF] Caliot | Model of Spectral and Directional Radiative Transfer in Complex Urban Canopies with Participating Atmospheres[END_REF]. This interdependency between the local climate and the urban structure eventually adds on the inherent spatiotemporal variability of the solar resource. The geographic location, daily and seasonal cycles, via the apparent movement of the sun relatively to the Earth, and meteorological conditions, especially through the motion of clouds, are mainly responsible for the natural variations of sunlight. All of this leads to a significant spatiotemporal variability of the solar resource at the urban scale. Various approaches and a plethora of increasingly efficient simulation tools have been proposed for the modelling of radiative exchanges taking place in the urban fabric [START_REF] Freitas | Modelling solar potential in the urban environment: State-of-the-art review[END_REF][START_REF] Thebault | A Comparative Study Of Simulation Tools To Model The Solar Irradiation On Building Façades[END_REF]. However, when dealing with multivariate phenomena involving a huge amount of data over a wide range of spatial and temporal scales, like in the present situation, improvements in the numerical capabilities are still usually not sufficient.

A solution to deal with such complex physical problems is to apply Model Order Reduction (MOR) techniques [START_REF] Schilders | Model Order Reduction: Theory, Research Aspects and Applications[END_REF]. MOR has been widely employed in various fields for models defined in high-dimensional spaces (e.g., quantum mechanics, fluid dynamics) or for multiscale problems with very large characteristic spatial or temporal scales (e.g., solid mechanics, acoustics) [START_REF] Chinesta | A Short Review on Model Order Reduction Based on Proper Generalized Decomposition[END_REF]. In multiple situations, a high-dimensional physical phenomenon can be reasonably approximated by the behaviour of a few numbers of its dominant features [START_REF] Schilders | Model Order Reduction: Theory, Research Aspects and Applications[END_REF]. The benefit of MOR in such cases is thus twofold:

1. The essential features are identified, which inherently facilitates the interpretability of the global model's behaviour;

2. The restricted set of features can be used to create a reduced order model of the physical problem and drastically lessen the computational effort for future modelling.

Among the variety of available MOR approaches, Principal Component Analysis (PCA) [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF] has been extensively used and validated in multiple areas, including renewable energies [START_REF] García-Cuesta | A combination of supervised dimensionality reduction and learning methods to forecast solar radiation[END_REF][START_REF] García-Hinde | Evaluation of dimensionality reduction methods applied to numerical weather models for solar radiation forecasting[END_REF]. Also known as Principal Orthogonal Decomposition (POD) [START_REF] Berkooz | The Proper Orthogonal Decomposition in the analysis of turbulent flows[END_REF] in computational fluid dynamics, Hotelling Transform [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF] and Karhunen-Loève Transform [START_REF] Loève | Probability theory II[END_REF] in image processing or the Empirical Orthogonal Function (EOF) [START_REF] Lorenz | Empirical orthogonal functions and statistical weather prediction[END_REF] in climate science, it consists of identifying the most representative orthogonal features of a given dataset [START_REF] Tu | Dynamic Mode Decomposition: Theory and Applications[END_REF], namely directions in which data varies. PCA generates a set of uncorrelated features (principal components) being linear combinations of the original variables and summarising at best the dataset, i.e., accounting for the maximum variance [START_REF] García-Hinde | Evaluation of dimensionality reduction methods applied to numerical weather models for solar radiation forecasting[END_REF]. It is robust to sets of data dealing with multiple interdependent variables and with a high dimensionality with respect to their number of samples [START_REF] García-Cuesta | A combination of supervised dimensionality reduction and learning methods to forecast solar radiation[END_REF].

The aim of this paper is to present and investigate the applicability of Model Order Reduction methods as an efficient means to analyse and characterise the spatiotemporal variability of the solar resource within an urban context. In that objective, a numerical simulation of the distribution of annual solar radiations, in terms of irradiances, has first been conducted on a simple theoretical three-dimensional urban geometry of heterogeneous morphology. A subsequent statistical analysis has been performed by means of PCA on the simulated field of irradiances. The main methodology is detailed in section 2. The principal resulting features have finally been analysed with respect to their influence on the model's spatial and temporal variabilities (sections 3 and 4).

METHODOLOGY

This section details the methodology and corresponding methods adopted in this work, as well as their underlying rationales and assumptions.

Numerical model

The input geometry has consisted of an assembly of nine three-dimensional buildings, namely one central and eight peripheral structures (Fig. 1). The motive for this configuration was to keep the domain size reasonable for maintaining an acceptable trade-off between speed and accuracy in the simulation process. The main design parameters, i.e., Shape Factor, Floor Area Ratio, Site Coverage and Average Building Height [START_REF] Natanian | A parametric approach to optimizing urban form, energy balance and environmental quality: The case of Mediterranean districts[END_REF], have been defined in the aim of mimicking a theoretical urban district with an open mid-rise typology [START_REF] Stewart | Local Climate Zones for Urban Temperature Studies[END_REF]. The morphological diversity has been restricted here to a heterogeneity in the buildings location and height. They have been arranged in a disorganised manner across the 100 × 100 m 2 district area and modelled by cuboids of 20 × 20 m 2 bases with a height ranging from 15 to 35 m. For the sake of simplicity, no semi-transparent (e.g., vegetation) or highly specular (e.g., glazings) elements, which would alter the reflected, transmitted and absorbed components of incoming radiations at the district scale, have been introduced. All surfaces have been considered Lambertian, i.e., perfect diffusers, with a reflectance of 0.3. Thus, both the directional and spectral dependencies of the full geometry on incident radiations have been omitted.

Similar atmospheric profiles to Caliot et al. (2022) [START_REF] Caliot | Starter Pack and validation scripts for the Monte-Carlo-based urban radiative transfer model described in "Model of Spectral and Directional Radiative Transfer in Complex Urban Canopies with Participating Atmospheres[END_REF] have been employed to describe the atmosphere, where the optical properties of its constituting gases are defined using the ECRAD radiation scheme for the European Centre for Medium-Range Weather Forecasts (ECMWF) [START_REF] Hogan | ECRAD: A new radiation scheme for the IFS[END_REF]. No clouds have yet been consider at this stage, resulting in typical clear sky conditions for the location of Paris (48 • 51 ′ 12.28 ′′ N; 2 • 20 ′ 55 ′′ E).

Annual solar irradiances modelling

A new engine, named htrdr-urban, has recently been proposed and validated for the accurate simulation of solar and infrared radiative transfers in an urban context with complex geometries and in the presence of a cloudy atmosphere [START_REF] Caliot | Model of Spectral and Directional Radiative Transfer in Complex Urban Canopies with Participating Atmospheres[END_REF][START_REF] |méso|star> | htrdr[END_REF]. By combining backward Monte-Carlo ray tracing with Null-Collision Algorithms (NCA), the computational effort becomes quasiindependent of the complexity of the three-dimensional scene. This tool has thus been used here to evaluate the shortwave irradiances across the measurement grids at each timestep, by solving the monochromatic Radiative Transfer Equation (RTE) for the defined simulation setup. No heat transfers by conduction or convection have been considered, meaning that the simulation boundary conditions at the domain borders are entirely established by the geometry and surfaces albedo of the numerical model.

The simulation has been performed for one building of interest located at the centre of the studied neighbourhood, highlighted in green in Figure 1. The distribution of incident radiations on the considered structure has been evaluated in terms of horizontal or vertical irradiances (W.m -2 ) for a set of measurement points across its visible surfaces, i.e., for its four vertical façades (east, west, north, south) and roof. Each surface has been discretised into a regular array of points. Similarly to previous studies using the same geometry [START_REF] Thebault | A Comparative Study Of Simulation Tools To Model The Solar Irradiation On Building Façades[END_REF][START_REF] Govehovitch | Numerical Validation of the Radiative Model for the Solar Cadaster Developed for Greater Geneva[END_REF], a spatial resolution of t ϕ = 1 m has been adopted here to allow for a reasonable observation of variations in the spatial patterns of the solar resource. Values have been computed for a full year, from 21 st December 2020, 08:30 to 20 th December 2021, 15:00, with a timestep t α = 15 min.

Model Order Reduction

Principal Component Analysis has been employed to decompose the simulated field of irradiances along its principal features both in the temporal and spatial domains. The rationale behind this decomposition builds on the spatiotemporal separated representation suggested by Pierre Ladevèze a few decades ago [START_REF] Chinesta | A Short Review on Model Order Reduction Based on Proper Generalized Decomposition[END_REF], by assuming that the incident irradiance distribution on a given surface can be expressed as:

E (x, t ) = k≥1 γ k ϕ k (x)α k (t ) (1)
with ϕ k and α k the spatial and temporal contributions, respectively, and γ k a real number vouching for the field amplitude at rank k. In that aim, the raw dataset of simulated irradiances E (x, t ) has initially been arranged in a 2D real matrix of size n × p:

E =    features e 11 • • • e 1p . . . . . . . . . e n1 • • • e np    n,p        samples (2)
Its rows (samples) and columns (features) refer to the n annual timesteps and to the p measurement points, respectively. The number of discretisation points being here significantly higher in the temporal than in the spatial domain (n ≫ p), it results in a very "thin" input matrix E.

When performing PCA, it is a common requirement that data are standardized, i.e., centred around a null mean and scaled to have a standard deviation of one, prior to the reduction procedure. Translating the values by subtracting their means ensures that the first principal component actually reflects the maximum variance of the dataset and not its mean value [START_REF] Tu | Dynamic Mode Decomposition: Theory and Applications[END_REF], which would otherwise be likely to make it statistically misleading. Scaling is usually employed to avoid allocating too much weight to some features over the others and reduce noise in the higher-order components [START_REF] Alexandris | Remote sensing of burned areas via PCA, Part 1; centering, scaling and EVD vs SVD[END_REF]. PCA is in that case performed on the correlation instead of the covariance matrix. Within the scope of this work, scaling would not be of high relevancy, since all features are by definition on the same scale. Hence, data centring only has been carried out during the pre-processing stage by deducting to each column j of E its mean value µ α, j .

PCA can either be performed via Eigenvector Decomposition (EVD) of the covariance (or correlation) matrix for the input data or by applying Singular Value Decomposition (SVD) to the raw dataset directly. In EVD-based PCA, the principal components are accessed via eigenvectors of the covariance matrix (principal directions) and their related amount of information is given by its eigenvalues. When SVD is employed, the principal directions, i.e., eigenvectors, correspond to the singular vectors of the input matrix, while the singular values are equal to the square root of the eigenvalues (standard deviations). The latter aspect justify why PCA is generally more accurate and numerically robust, yet slightly slower, when based on SVD rather than on EVD for a higher number of samples than features. SVD-based PCA is therefore preferred for high-dimensional datasets [START_REF] Alexandris | Remote sensing of burned areas via PCA, Part 1; centering, scaling and EVD vs SVD[END_REF] and has been used thereafter. Application of SVD to the centred matrix Ė provides the following decomposition:

Ė = U.Σ.V ⊤ (3)
with U a unitary matrix of size n ×n and Σ a n × p matrix whose main diagonal contains the singular values σ α, j associated with the singular vectors represented by the rows of the p × p matrix V ⊤ 1. Having here n ≫ p, only the first p columns and rows of U and Σ, respectively, related to the p non-null singular values σ α, j , have been kept to avoid unnecessary huge sizes for these matrices. Equation (3) can thus be rewritten:

Ė′ = U ′ .Σ ′ .V ⊤ (4)
where U ′ is now of size n × p and Σ ′ of size p × p. The principal components, i.e., the projection of the data onto the principal directions, have then been obtained by multiplying Ė′ by V:

Ė′ V = U ′ .Σ ′ .V ⊤ .V = U ′ Σ ′ = pc α,1 • • • pc α,p n,p (5)
with their related variances, i.e., the eigenvalues λ α, j of the covariance matrix of Ė′ , being directly accessed from the singular values σ α, j :

λ α, j = σ 2 α, j n -1 (6)
Finally, the principal component pc α, j of rank j is simply selected from the columns of Ė′ V = U ′ Σ ′ , while its rows give the samples coordinates in the projection space. To visualise the corresponding mode Êα,j , i.e., the reconstruction of the original dataset from this principal component only, pc α, j has further been multiplied by its related principal direction v * j (j -th column in V ⊤ ), and its mean value µ α, j summed back:

Êα,j = pc α, j .v * j + µ α, j I (7)
The above demonstration is valid within the scope of a temporal decomposition, with the n annual timesteps and the p grid points accounting for the samples and the features, respectively, in the PCA calculation. The idea being to investigate both the temporal and spatial variabilities of the field of irradiances, E has also been transposed so that its rows now depict the p evaluation points (samples) and its columns the n related timesteps (features). The input dataset E ⊤ for spatial decomposition would thus be a "thick" matrix of size p ×n, with µ ϕ,i the mean values over its temporal features (columns). The role of U and V would be reversed with respect to the illustrated case. Though p ≫ n in this configuration, the size of the dataset employed within the scope of this work was still reasonable for the SVD procedure to remain fast enough.

The first modes resulting from the temporal ( Êα,j ) and spatial ( Êϕ,i ) decompositions have then be examined to assess the variabilities of the incident field of irradiances (section 3).

1V ⊤ refers here to the transpose of the real matrix V

RESULTS

Results from the modal decomposition of the incident field of irradiances on the different facets of the test building are presented here. The spatial and temporal decompositions for the first four modes on the east façade (Fig. 2) and for the first modes on each of the orther surfaces, i.e., west, north, south and roof, (Fig. 3) have been selected by the authors for visualisation. Their amplitudes have all been normalised between -1 and 1 and centred around 0 using a diverging colour scale to facilitate the analysis. Only modes up to the fourth rank are investigated, since they capture, by definition, most of the variability (cf. 3.3). It is not excluded that higher order modes account for interesting behaviours, yet they will not be analysed within the scope of this preliminary study.

Temporal decomposition

The temporal decomposition for each mode of interest is represented on the left column in Fig. 2 and 3 by the normalised irradiance received at a given time in a day (24 hours) as a function of a speficic day in the year (365 days). As it could have been expected, a first common trend to the temporal modes across the different surfaces is that they depict the variabilities in the distribution of irradiance at different time scales, with the level of details provided by each of them increasing with their order. Namely, Êα,1 , which accounts for the highest explained variance (EV ), generally provides information about the global annual irradiance distribution and highlights the seasonal repartition (winter/summer) (Fig. 2, 3 -Êα,1 ). Subsequent higher-dimensional modes represent finer variations, from the daily sharing between the morning and the afternoon, including midday in some cases, and up to an explanation at the hour scale (Fig. 2 -Êα,2-4 ). Êα,4 for the east, west and south façades further shows the influence of equinoxes and solstices on the daily irradiance distribution (Fig. 2 -Êα,4 ).

Rationally, the distribution for each temporal mode is similar yet symmetric with respect to midday for the east and west façades (Fig. 2,3). All of the first four dominant modes for both surfaces also share very close respective variances. Similarly, modes for the south and north sides seem to be complementary (Fig. 3), taking into account the daily course of the sun. However, the amount of information is interestingly more concentrated within the first mode for the north than for the south orientation. In general, this trend has been observed across surfaces of different orientations, with notable differences in the explained variances provided by modes of a same decomposition rank k.

Spatial decomposition

The spatial decomposition for the same modes as detailed in section 3.1 is represented on the right column in Fig. 2 and 3 by the normalised irradiance received at a given position on the considered facet. Overall, results for the spatial decomposition correlate with the above observations in the temporal domain. However, more similarities have been found between the spatial distributions of the first dominant modes from the different surfaces. All the first four modes for the east and west façades now depict almost exactly the same allures (not shown here). Their respective explained variances are also very similar and close enough to the values of the corresponding modes for the south façade. The total variance spreads more equally over the four dominant modes here than for the temporal decomposition, the amount of information held by the first mode alone being lower in the present situation (Fig. 3).

Interestingly, similar modal distributions are encountered across multiple vertical surfaces but at different ranks k, e.g., Êϕ,3 for the north side is similar to Êϕ,2 for the east and west orientations (not shown here).

The influence of surrounding geometries on the distribution of the field of irradiances is more clearly highlighted here. This is especially illustrated by the slight variations in the distributions of modes of same rank k between the east and west surfaces (Fig. 2, 3), which would not have been observed in a completely symmetric area under clear sky conditions. The differences in the heights and positions of neighbouring buildings on the western and eastern sides of the considered structure may have induced these variations.

Scree Plots

Finally, scree plots, which express the explained variance as a function of principal components, have been computed for both the temporal and spatial decompositions for each facet. For all surfaces and domains, it has been observed that a substantial part of the total explained variance is contained in the first few components, i.e., modes. In the case of the east-facing façade, modes after the tenth already account for less than 0.5 % of the total explained variance each (Fig. 4).

DISCUSSION & CONCLUSIONS

A novel approach for the analysis of solar radiations variabilities in urban environments by means of Principal Component Analysis has been proposed. The field of irradiances has been simulated onto the façades and roof of a cuboidal building within an urban context of heterogeneous structures. First results show that modal decomposition provides a very efficient tool to analyse the first dominant modes of the variability distribution in both the temporal and spatial domains.

Overall, a reasonable consistency has been observed between results from the temporal and spatial decompositions. Lower-dimensional modes, i.e., with the highest variances, have been found to account for the variations on wide scales while modes at higher ranks describe the variations on smaller scales. In the complex urban environment where a plethora of parameters are likely to play a role in the variability of the solar resource, MOR methods would greatly facilitate the identification of the most influential parameters at characteristic spatiotemporal scales. Especially, the contribution of specific surrounding geometries to the disturbance of the incident field of irradiances has been pinpointed by certain spatial modes. The orientations of the evaluation surfaces have clearly been represented in the decomposition, with complementarities and/or similarities between facets on the same axes (north/south or east/west). Results for the east-and west-facing façades have presented great consistency across the different decomposition modes. In general, a different amount of information was provided by modes of same dimensions for the north and south façades, and compared to the east/west axis. This suggests that specific characteristic scales contribute more to the global spatiotemporal variations for some orientations than for others. The distribution of certain spatial modes has also been found to be common to several of the vertical surfaces, but at different orders of decomposition. This implies that same parameters are involved in the spatial variability of the field of irradiances across the surfaces, yet to a different extent depending on their orientation. These findings would therefore be of high relevancy to identify areas subjected to high spatial and temporal variabilities of the solar resource.

Scree plots have further shown that the first few principal components were sufficient to capture a significant part of the information in our physical model (above 96.5 % and 84, 3 % of cumulative EV held by the first 10 components for the temporal and spatial decompositions, respectively), paving the way for the definition of a reduced order model that would be very close to the original. These results therefore demonstrate that the use of PCA, and MOR methods in a wider extent, can be very promising approaches to analyse the spatiotemporal variabilities of the irradiance field, which would be otherwise very complex given the high spatial and temporal definitions.

This work has led to promising findings but still presents some limitations, which should be pointed out.

The main limitation comes from the numerical model itself. The geometry has been defined simple on purpose, with a limited number of variables. The objective was here to validate the applicability and efficiency of MOR methods for the analysis of solar radiations variabilities in a heterogeneous environment. The consideration of more complex (e.g., shapes, inclination of surfaces), diverse (e.g., vegetation) and numerous geometries, with non-trivial reflective properties (e.g., specular, semi-transparent), would have likely impacted the results and the difficulty of their analyses. The influence of clouds has neither been considered, since a simple clear sky model has been used here. Besides, the analysis has been restricted to the first four modes for each orientation. Though the authors have sensibly based their decision from the results of preliminarily computed scree plots, higher-dimensional modes could have provided more insights on specific aspects.

Future suggested work would be therefore to apply the proposed methodology to more complex numerical models, with geometrical and physical properties closer to a real urban neighbourhood. Especially, the variations in common morphological factors like the domain size or the buildings dimensions and arrangement, which would have an apparent impact on the distribution of incident irradiances, are intended to be investigated by the authors. Non-Lambertian surfaces, e.g., highly specular glazings or semi-transparent materials, are also common features of modern urban environments, and their effect on the incoming radiations is worth being examined via the proposed method. The influence of different sky types should also be explored, to consider the effect of clouds in the simulation of solar radiations (e.g., movement, radiative properties). The spatiotemporal resolutions for the simulation could also be finer, e.g., minutely and sub-meterly in the temporal and spatial domains, respectively. Regarding the modal decomposition, other MOR methods should be investigated. Although the relevancy and efficiency of PCA has been evidenced in this study, other techniques may prove relevant for this specific problem. For instance, Proper Generalized Decomposition (PGD) is based on separated representations, e.g., of space and time, and would therefore be a prime candidate in such situations [START_REF] Chinesta | A Short Review on Model Order Reduction Based on Proper Generalized Decomposition[END_REF]. Similarly, Dynamic Mode Decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], which could be sensibly apprehended as a combination of PCA and Fourier transform, has been proposed specifically to analyse and capture dynamics of evolving physical fields, e.g., fluid flows [START_REF] Tu | Dynamic Mode Decomposition: Theory and Applications[END_REF]. The use of deep neural networks for the selection and extraction of the optimal features is also of current interest [START_REF] Higashiyama | Feature extraction of numerical weather prediction results toward reliable wind power prediction[END_REF], though still mainly employed for the post-processing of reduced order models. The nature of the decomposition method remains crucial, as on it depends the meaning of the decomposition modes, so does the underlying physics [START_REF] Tu | Dynamic Mode Decomposition: Theory and Applications[END_REF].

Fig. 1

 1 Fig. 1 3D (a) and top view (b) of the considered heterogeneous district. Each building is coloured in light grey and referenced in terms of its height (m). The building of interest is highlighted in green. A buffer band of 10 m width surrounds the district area (salmon-pink).

2 -Ê φ, 3 -Ê φ, 4 -EASTFig. 2 Ê φ, 1 -xFig. 3

 234213 Fig. 2 Temporal (left column) and spatial (right column) decompositions of the incident field of irradiances on the east façade of the test building. The first four modes Ê1-4 are displayed row by row in growing order. Spatial axes are defined with respect to the coordinate system used in Fig. 1.

Fig. 4

 4 Fig. 4 Scree plot for the temporal decomposition on the east façade of the test building.

ACKNOWLEDGMENT

This work has been supported by the French National Research Agency, through the Investments for Future Program (ref. ANR-18-EURE-0016 -Solar Academy). The research unit LOCIE is member of the INES Solar Academy Research Center.


Claire Dumas

Régine Basséguy

Alain Bergel

Electrochemical oxidation of acetate by Geobacter sulfurreducens : Influence of anode materials

Electrochemical oxidation of acetate by Geobacter sulfurreducens : Influence of anode materials

Dumas C.,Basséguy R., Bergel A. (1) (1) Laboratoire de Génie Chimique-CNRS, Basso Cambo, 5 rue Paulin Talabot, BP 1301, 31106 Toulouse, France

The recent discovery of micro-organisms that could transfer directly electrons to insoluble electrode [1] had encouraged the researches on microbial fuel cells. The interfaces anode/micro-organism are explored to enhance the bacterial electron transfer and increase the power densities produced: some modified anodes have been designed, binding for example electron mediators (Mn 4+ , neutral red …) to woven graphite [2]. Moreover the porosity of different forms of carbon has been found to influence current productions [3]. The purpose of this work was to compare the efficiency of graphite, DSA® and stainless steel as anode materials for growing Geobacter sulfurreducens under imposed potential conditions. Current densities obtained were seven times higher than values presented in literature [4] and depended on anode material. Analyses of the electrode surface with techniques like optical interferometry, SEM or confocal microscopy revealed significant differences in roughness and bacterial coverage, which may explain the variety of current densities obtained.