Sophie Achard
email: sophie.achard@univ-grenoble-alpes.fr

Pierre Borgnat

Irène Gannaz

Marine Roux

Wavelet-based graph inference using multiple testing

Keywords: wavelets, multiple testing, correlation, graph inference

published or not. The documents may come

INTRODUCTION

Brain connectivity analysis is very important in order to understand the mechanism implicated in brain pathologies such as for example schizophrenia, [START_REF] Lynall | Functional connectivity and brain networks in schizophrenia[END_REF] Parkinson, 2 strokes, 3 traumatic brain injury [START_REF] Silva | Disruption of posteromedial large-scale neural communication predicts recovery from coma[END_REF] . . . With the use of non invasive recording of brain functioning it is nowadays easy to collect data on groups of patients covering the whole brain. Usually, the brain is parcellated in regions of interest with the objective to identify the possible pairs of brain regions dependent to each other. [START_REF] Richiardi | Machine learning with brain graphs: Predictive modeling approaches for functional imaging in systems neuroscience[END_REF] This creates a graph where the nodes are the brain regions and the edges represent the dependence between two brain regions. The parcellation of the brain is well documented with several papers proposing different parcellation schemes. [START_REF] Margulies | Resting developments: a review of fmri post-processing methodologies for spontaneous brain activity[END_REF] The choice for the definition of dependence between nodes is more complicated. As temporal signals are extracted from brain regions, the dependence is measured via a statistical estimation facing classical statistical challenges such as number of available samples, long-memory time series, multiple comparisons. . . A poor performance of statistical estimations may have dramatic effect on the resulting extracted graphs and consequently on the graph metrics used to compare groups of patients and controls. [START_REF] Achard | A resilient, low-frequency, smallworld human brain functional network with highly connected association cortical hubs[END_REF][START_REF] Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Patel | A wavelet-based estimator of the degrees of freedom in denoised fMRI time series for probabilistic testing of functional connectivity and brain graphs[END_REF][START_REF] Malagurski | Topological disintegration of resting state functional connectomes in coma[END_REF] In this paper, we provide a whole method from signals to graphs in order to control the misdetection of edges. The dependence is estimated using wavelet correlations. [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF][START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] For a given frequency band, pairs of correlations are estimated using the wavelet coefficients. Then, thanks to the good convergence properties of the estimator, p-values are extracted for each pair. The effect of making error on identification of edges is displayed in Figure 1. Two different paradigms on classical graph metrics are illustrated via global efficiency and local efficiency. Adding wrong edges at random has a dramatic effect on the graph metrics, however, missing true edges is less dramatic. The selection of significant pairs of correlation is then derived so as to control asymptotically the family wise error rate, i.e. the probability of having more than one wrong edge. In practice, we choose to apply the classical Bonferroni procedure on the p-values of wavelet correlations. The paper is organized in 4 different parts. The first part is dedicated to the process to start from brain signals and end up with correlation matrices. The second part gives details on the multiple testing approach in order to extract the graphs. Simulations are proposed in a third part, where we took the classical example of small-world graphs. Finally, our approach is applied to a real fMRI dataset on rats.

FROM SIGNALS TO CORRELATION MATRICES USING WAVELETS

Understanding multivariate time series is the main motivation of our paper. It is now well recognized that the time series observed from brain recordings are well modeled using long-memory time series. [START_REF] Meyer | Wavelet-based estimation of a semiparametric generalized linear model of fmri time-series[END_REF][START_REF] Maxim | Fractional Gaussian noise, functional MRI and Alzheimer's disease[END_REF] Also, from a neuroscience point of view it is common to study the low frequencies part of the signals. Wavelets are perfect tools to both take into account the long-memory and focus on low frequencies. Indeed, wavelets are acting as frequency decompositions of the signals. Besides several results on the good performances of wavelets estimations for longmemory time series have been obtained in a univariate context. [START_REF] Moulines | A wavelet Whittle estimator of the memory parameter of a nonstationary gaussian time series[END_REF][START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF] The estimation of interactions between brain regions is described in this part using wavelet correlation estimations.

Multivariate long-memory models

Let X = {X (k), k ∈ Z, = 1, . . . , p} be a multivariate stochastic process. Each process X is not necessarily stationary. Denote by ∆X the first order difference, (∆X )(k) = X (k)-X (k-1), and by ∆ D X the D-th order difference. For every component X , there exists D ∈ N such that the D -th order difference ∆ D X is covariance stationary. Following, [START_REF] Chambers | The simulation of random vector time series with given spectrum[END_REF][START_REF] Moulines | On the spectral density of the wavelet coefficients of longmemory time series with application to the log-regression estimation of the memory parameter[END_REF][START_REF] Achard | Fractal connectivity of long-memory networks[END_REF] we consider a long-memory process X with memory parameters d = (d 1 , d 2 , . . . , d p ). For any D > d -1/2, we suppose that the multivariate process Z = diag(∆ D , = 1, . . . , p)X is covariance stationary with a spectral density matrix given by for all ( , m) , f

(D ,Dm) ,m (λ) = 1 2π Ω ,m (1 -e -iλ ) -d s (1 -e iλ ) -d s m f S ,m (λ), λ ∈ [-π, π],
where the long-memory parameters are given by d S m = d m -D m for all m. The functions f S ,m (•) correspond to the short-memory behavior of the process. We assume that f S ,m (0) = 1 and that there exist L > 0 and β ∈ (0, 2) such that for all ( , m),

max λ∈[-π,π] |f ,m (λ) -1| ≤ L|λ| β .
The generalized cross-spectral density of processes X and X m can be written as

f ,m (λ) = 1 2π Ω ,m (1 -e -iλ ) -d (1 -e iλ ) -dm f S ,m (λ), λ ∈ [-π, π].
A process X satisfying the above conditions will be denoted in the following as a M (d) process.

Wavelets decomposition

Let (φ(•), ψ(•)) be respectively a father and a mother wavelets. Their Fourier transforms are given by φ

(λ) = ∞ -∞ φ(t)e -iλt dt and ψ(λ) = ∞ -∞ ψ(t)e -iλt
dt. At a given resolution j ≥ 0, for k ∈ Z, we define the dilated and translated functions

φ j,k (•) = 2 -j/2 φ(2 -j •-k) and ψ j,k (•) = 2 -j/2 ψ(2 -j • -k).
Thoughout the paper, we adopt the same convention as in [START_REF] Moulines | On the spectral density of the wavelet coefficients of longmemory time series with application to the log-regression estimation of the memory parameter[END_REF] and, [START_REF] Moulines | A wavelet Whittle estimator of the memory parameter of a nonstationary gaussian time series[END_REF] that is, large values of the scale index j correspond to coarse scales (low frequencies).

Let X(t) = k∈Z X(k)φ(t -k). The wavelet coefficients of the process X are defined by

W j,k = R X(t)ψ j,k (t)dt j ≥ 0, k ∈ Z. For given j ≥ 0 and k ∈ Z, W j,k is a p-dimensional vector W jk = W j,k (1) W j,k (2) . . . W j,k (p) , where W j,k ( ) = R X (t)ψ j,k (t)dt.
The regularity conditions on the wavelet transform are expressed in the following assumptions. They will be needed throughout the paper.

(W5) For all i = 1, . . . , p, (1 + β)/2 -α < d i ≤ M .
For all j ≥ 0, we will denote also by H j the squared gain function associated to wavelet filters at scale j.

In practice, a finite number of realization of the process X, say X(1), . . . X(n), is observed. Assumption (W1) ensures that the wavelets have a compact support. Thus at each scale j only a finite number n j of coefficients are non zero.

Wavelet correlation

Statistical convergence results of variance, covariance and correlation have already been established e.g by 20 and 11 for finite difference processes. Confidence intervals are then derived and statistical hypotheses tests can be constructed. Proposition 2.1 (Convergence of wavelet correlations). Let X(1), . . . , X(n) be n observations of a p-dimensional real valued M (d) process with a spectral density f . Define for all , 1 ≤ ≤ p, {W j,k ( ), k ∈ Z} the wavelets coefficients at scale j of X . Assume that for all , 1 ≤ ≤ p, the process {W j,k ( )} k is a Gaussian process of spectral density D (j) 0;0 (.; , ) so that D (j) 0;0 (f ; , ) > 0 almost everywhere, and the bivariate process {W j,k ( ), W j,k (m)} k is Gaussian with spectral density D (j) 0;0 (.; , m) so that D (j) 0;0 (f ; , m) > 0 almost everywhere. Under the hypotheses (W1)-(W5), the empirical estimator of correlation of wavelets coefficients at scale j is unbiased and is asymptotically Gaussian distributed:

ρ ,m (j) := θ ,m σ (j) σ m (j) = nj k=0 (W j,k ( )W j,k (m)) ( nj k=0 W 2 j,k ( ) nj k=0 W 2 j,k (m)) 1/2 L -→ N ρ ,m (j), Σ cc j (2.1)
where ρ ,m (j) := γ ,m (j) (γ , (j)γ m,m (j)) 1/2 , with γ ,m (j) := H j (λ) f ,m (λ) d λ.

(2.

2)

The matrices Σ cc j are given in formula (5.4) of. [START_REF] Whitcher | Assessing nonstationary time series using wavelets[END_REF] The symbol L -→ corresponds to the convergence in distribution when n tends to infinity.

Brain recordings often present long-range dependence, with non homogeneity of long-memory parameters d. In these cases the tested values {ρ ,m (j), 1 ≤ < m ≤ p, j ≥ 0} depend on the parameters d. We refer to Proposition 1 and Proposition 2 of. [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] More precisely, suppose that at a first order approximation the spectral density satisfies, when λ → 0 + ,

f (λ) ∼ Σ(λ), Σ(λ) = Λ(d) * ΩΛ(d) , with Λ(d) = diag(λ -d e -iπd/2 ) ,
where the * exponent denotes the complex conjugate. This approximation holds for M (d) processes defined in Section 2.2. Let r be the correlation coefficients of Σ(λ), r ,m = Σ ,m (λ)/ Σ , (λ)Σ m,m (λ) = Ω ,m e i π(d -dm)/2 . When ψ is real, Proposition 2 of [START_REF] Achard | Multivariate wavelet whittle estimation in long-range dependence[END_REF] shows that asymptotically for all ( , m), ρ ,m (j) tends to the real part of r when j tends to infinity, that is,

ρ ,m (j) -→ j→∞ Re(r ,m ) = Ω ,m cos(π(d -d m )/2). (2.3)
An interesting perspective is to build an estimation procedure on {Ω ,m } ( ,m) rather then {r ,m } ( ,m) . In particular, notice that when d -d m is close to 1 modulo 2, the correlation of the wavelet coefficients ρ ,m (j) tends to 0 even if Ω ,m = 0.

Fisher transform of wavelet correlation

The exact formula of the asymptotic variance of wavelet correlation 2.1 is complicated and it is tempting to get rid of this computation. It is possible by using the Fisher transform of the wavelet coefficients. Simply as in the usual case, the Fisher transform is stabilizing the variance of the estimator. In our case, this facilitates the computation of the confidence intervals where the variance is now simply a function of the number of samples. [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF] Let z be the Fisher transform: for all x, -1 ≤ x ≤ 1,

z(x) := 1 2 log 1 + x 1 -x = tanh -1 (x).
Proposition 2.2 (Convergence of correlation with Fisher transform). Under the same hypotheses of proposition 2.1, the estimator of empirical correlation verifies

(n j -3)(z( ρ ,m (j)) -z(ρ ,m (j))) L -→ N (0, 1) ;
(2.4)

In the following, we will apply this result for constructing the multiple testing framework to test a large sample of pairwise correlations at a given wavelet scale j, that is,

H 0, m (j) : ρ ,m (j) = 0 against H 1, m (j) : ρ ,m (j) = 0, 1 ≤ < m ≤ p, j ≥ 0.
(2.5)

MULTIPLE TESTING FRAMEWORK

In the previous section, we proposed a procedure which gives correlation matrices, where each entry ρ ,m (j) is the estimated correlation between brain regions and m at scale j. To infer the dependence graph, there is a need to select significant correlations based on a correct statistical approach.

Let V denote the set of the indexes of the brain regions, V = {1, . . . , p}. For a given scale j, define graphs G j as a pair G j = (V, E j ) where the set of edges E j is a subset of V × V . We consider undirected graphs, that is, if ( , m) ∈ E j , then (m, ) ∈ E j . The dependence graph is obtained by defining E j = {( , m), ρ ,m (j) = 0}.

From the estimated correlation matrices obtained in Section 2, our objective is to estimate graphs G j for each scale j. Denote G j = (V, E j ) the estimated graphs. We hence test H 0, m (j) : ρ ,m (j) = 0 against H 1, m (j) : ρ ,m (j) = 0, for all ( , m) ∈ V. The set of edges E j is defined as E j = {( , m) ∈ V × V, H 0, m (j) is rejected}. For each graph G j , estimation G j is obtained applying ν = p (p -1)/2 tests (since sets E j are symmetric). We now detail the procedure of testing.

Multiple hypothesis tests

Consider the specific case of correlations. Let X (n) = {X(1), . . . , X(n)} be independent realizations from a random vector X = (X ( ) ) =1,...,p with values in R p . Define Γ = (ρ ,m ) ,m=1,...,p the correlation matrix of X, with ρ ,m = Cor(X ( ) , X (m) ). Consider the two-sided testing problem:

H 0, m : ρ ,m = 0 against H 1, m : ρ ,m = 0, ( , m) ∈ H, (3.1)
where H is the set of indexes tested. In the following

H = {1 ≤ < m ≤ p}. The number of tests corresponds to the cardinality of H, that is, ν = p (p -1)/2.
The objective of multiple testing procedure is to give a rejection set

R = {( , m) ∈ H, (H 0, m ) rejected},
such that a given error is controlled. We will consider here the type I error called Family Wise Error Rate (FWER), defined as

FWER(R) = P (∃( , m) ∈ R, ρ ,m = 0) .
For a given level α, the objective is to find R such that FWER(R) ≤ α. The rejection set is built on a p-value process associated to the tests. Note that in our particular case, the p-value processes associated to tests (3.1) present a dependence structure. [START_REF] Romano | Control of the false discovery rate under dependence using the bootstrap and subsampling[END_REF] In our proposed graph estimation, the rejection set R gives the estimated set of edges, E. That is, the estimated graph is G = (V, E) where E = R. The FWER criterion corresponds to the probability to have a false additional edge in the estimated graph. As we will illustrate in Section 3.3, controlling this quantity is important in order that usual metrics on the resulting graph remains satisfactory.

Asymptotic control of FWER

The Bonferroni procedure 25 is the most classical example of FWER control without the necessity to have independence of tests. A difficulty of the application of multiple testing in the context of correlations is that the tests are asymptotic. The Bonferroni method has then to be validated in this context.

Method (Bonferroni)

The Bonferroni multiple testing procedure applied on Fisher statistics is defined by

R bonf n,α = ( , m) ∈ H : p n, m ≤ α ν , (3.2)
where

p n, m = 2 1 -Φ √ n -3 z ρ m X (n) . (3.3)
Proposition 3.1. Let P denote a family of probability distributions of the random vector X. For the two-sided testing problem (3.1) based on the asymptotic p-value process (3.3), the method R bonf n,α provides an asymptotic control of the FWER at level α. That is, for all P ∈ P,

lim n→+∞ FWER R bonf n,α , P ≤ α. (3.4)
Proof. Let P ∈ P. Denote H 0 (P ) the set of indexes such that ρ ,m = 0 under distribution P , and ν 0 (P ) its cardinality. We have

Thus, lim

n→+∞ FWER R bonf n,α , P ≤ ( ,m)∈H0(P ) α/ν ≤ α ν 0 (P ) ν ≤ α.
This proposition is valid for any tests but in this paper, we focus our attention to correlations.

An equivalent formulation of Proposition 3.1 is written in terms of graphs in the following lemma, Lemma 3.1. Let P denote a family of probability distributions of the random vector X. For P ∈ P, let G(P ) = (V, E(P )) be the dependence graph associated to the data, that is,

E(P ) = {( , m) ∈ H, ρ ,m = 0}. The estimated graph is G n,α = (V, E n,α ) where E n,α = R bonf n,α .
Under the same hypotheses as in Proposition 3.1, for all P ∈ P,

lim n→+∞ P( G n,α ⊆ G(P )) ≥ 1 -α. (3.5)

Importance of detection of correct edges on graphs

Classical methods to compare graphs in neuroscience are using graph metrics such as minimum path length, clustering, and many others. [START_REF] Rubinov | Complex network measures of brain connectivity: uses and interpretations[END_REF][START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF] The precision of edge detections have to be quantified in terms of these classical metrics. We conducted simulations to observe the effect of missidentifications of edges on classical small-world graphs. We consider two different ways of missidentification of edges. The first one illustrated on the first row of Figure 1 corresponds to type II error i.e. not all the correct edges are observed, but the one present in the graph are correct. The second one illustrated on the second row of Figure 1 corresponds to type I error i.e. additional edges are added at random to the observed graph, and they are not correct. In the simulations, we considered graphs with 51 nodes. For the first set of simulated graphs (corresponding to the first row of Figure 1), 51 edges of a small-world graph with 204 edges are not observed. For the second set of simulated graphs (corresponding to the second row of Figure 1), 51 edges are added spuriously from a small-world graph with 102 edges. All generated graphs have finally 153 edges. We then compare the obtained graph with fully generated small-world graphs of 51 nodes and 153 edges. The objective is to observe the difference in metrics of the graphs under these two sets of simulations. The Number of edges is determined by the necessity to have the same number of vertices and edges to make the comparison possible. Based on the computation of two graph metrics, global efficiency and clustering, [START_REF] Achard | Efficiency and cost of economical human brain functional networks[END_REF] we can observe that the effect on graph metrics is more important when type I errors are made. Indeed, the difference in the distribution of global efficiency and clustering coefficient is less important between graphs with non observed edges and small-world graphs than between graphs with additional edges and small-world graphs. This justifies also our choice to control the FWER which is an efficient control of type I error.

Figure 1: Influence of missidentifications of edges on Global Efficiency and on Clustering coefficient. 10000 simulations were done on 51 nodes small-world graphs with 153 edges. Missidentifications consist in 33% of spuriousedges or 75% of only correctly observed edges, that is, in each case 51 edges misspecified. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q non obs. edges full 0.45 0.46 0.47 0.48 0.49

Global efficiency

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q non obs. edges full 0.1 0.2 0.3 0.4 0.5

Clustering

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q spurious edges full 0.45 0.46 0.47 0.48 0.49 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q spurious edges full 0.1 0.2 0.3 0.4 0.5

Application to wavelet correlation

In our context, a multiple testing procedure is set for each wavelet scales. At each scale j ≥ 0, the Bonferroni multiple testing procedure applied on Fisher statistics is defined by

R bonf n,α (j) = ( , m) ∈ H : p n, m (j) ≤ α ν , (3.6)
where

p n, m (j) = 2 1 -Φ √ n -3 |z ( ρ m (j))| , (3.7)
and ρ m (j) is defined as 2.1. Similarly to general correlation tests, we have an asymptotic control of FWER.

Proposition 3.2. Let X(1), . . . X(n) be n independent observation of a R p -valued M (d) process. Consider the two-sided testing problem (2.5) based on the asymptotic p-value process (3.7), evaluated with wavelet filters satisfying (W1)-(W5). Introduce a scale j 0 such that n2 -j0 tends to infinity when n goes to infinity.

Under the same assumptions of Proposition 2.1, the method R bonf n,α (j) provides an asymptotic control of the FWER at level α. That is, for all j ≥ j 0 ,

lim n→+∞ FWER R bonf n,α (j) ≤ α. (3.8)
The available number of samples is varying depending on the scales. Because of the asymtotic control of the FWER, when dealing with finite time series, it is not advisable to go for very small scales. The choice of the wavelet scale is complicated and will not be discussed in this paper. [START_REF] Achard | Wavelet-based and fourier-based multivariate whittle estimation: multiwave[END_REF] Besides the problem of size of available samples, the presence of short-memory affects the estimation of the graph. Indeed, for each scale j ≥ 0 the graph G j may be different. As the short-memory appears in the spectral density f , it influences the values of {ρ ,m (j), , m = 1, . . . , p} for low scales j. Note that, based on result (2.3), for sufficiently high scales, the graph G j does not depend on j. Thus, the estimation of the graph G * = (V, E * ) where E * = {( , m) ∈ {1, . . . , p}, Ω ,m ) = 0} may be better to consider. As shown in our previous paper, [START_REF] Achard | Wavelet-based and fourier-based multivariate whittle estimation: multiwave[END_REF] the choice of wavelet scales for the estimation of Ω should include scales not affected by the short-memory and scales with enough remaining points. We illustrate this purpose on simulations in Section 4.

SIMULATIONS

In order to illustrate the different properties obtained theoretically in the previous section, we ran simulations on a small-world graph. The true graph is constructed with 51 nodes and 102 edges with small-world generative models as implemented in the R package igraph. The choice of the number of nodes and edges is motivated by the real data example of Section 5.

Let us denote A the adjacency matrix of this graph. An example is given in Figure 2. Then, let us define the correlation matrix as C = I + ρA, where I is the identity matrix of size p. The parameter ρ is chosen such that C is positive definite. This is the case when |ρ| < 1/|λ s |, where λ s is the smallest eigenvalues of A. Simulations were done on a given graph adjacency matrix. Some extensions to numerical studies with different simulations of adjacency matrices are in progress. They are not displayed here for the sake of clarity.

Good behavior with a small-world generative models

Figure 3 is illustrating what we are expecting from theory. Indeed, the FWER is controlled as soon as the sample size is large enough. It is well known that controlling the FWER is not sufficient to have good results and the power of the procedure must be satisfactory. We also computed the power of our procedure. As illustrated in Figure 3 when the true correlation is too small, it is impossible to get an accurate solution, and no edges in the graph are detected. In order to get a correct solution, different parameters have to be taken into account: the size of the sample, the number of nodes in the graphs and the values of correlations in the correlation matrix. Similar conclusions were obtained e.g. by. [START_REF] Hero | Large-scale correlation screening[END_REF] Figure 3: FWER of edges inference with respect to the sample size for different values of ρ. 10000 simulations were done on 51 nodes small-world graphs. The grey interval around the value α corresponds to the confidence interval of a binomial distribution at level 95%. Schematically, there is now contradiction with the theory F W ER ≤ α even if the observed FWER is greater than α as long as it belongs to this interval. Power of edges inference : proportion of true edges detected with respect to the sample size for different values of ρ. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.02 0.04 0.06 250 500 750 1000 sample size n FWER rho q q q q q q 0.05 0.1 0.15 0.2 0.25 alpha q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Wavelet correlation testing and influence of short-memory on FWER

As an example, we simulate FIVAR(1,d) processes 31 X, defined as

X(k) = diag (1 -L) -d , = 1, . . . , p u(k),
where L is a lag operator and u follows a R p -valued VAR model:

I p + Φ L u(k) = (k) with (k) =    1 (k) . . . p (k)    i.i.d. with a Gaussian distribution N p       0 . . . 0    , C    .
The covariance matrix C satisfies a small-world model, as described at the beginning of the section. Similarly we consider p = 51 components. The matrix I p is the identity matrix in R p×p and Φ is a matrix in R p×p .

As stated previously, we are interested in the estimation of the graph G * = (V, E * ) where E * = {( , m) ∈ {1, . . . , p}, Ω ,m ) = 0}. In presence of long-memory, the values {ρ ,m (j), 1 ≤ < m ≤ p, j ≥ 0} depend on the long-memory parameters d. Due to equation (2.3) it is possible to apply a bias correction of the multiplicative term cos(π(d -d m )/2). However, for clarity, the influence of the heterogeneity of long-range dependence is not studied in our simulations. Then, we will take d constant, equal here to 0.2 in all our simulations.

To study the influence of short-range memory properties, we will consider three cases for the VAR part of the model: (a) First the matrix Φ in the FIVAR model will be taken equal to the zero matrix, giving simulations with no short-range memory: Φ ,m = 0 for all 1 ≤ , m ≤ p. As we want to recover E * = {( , m) ∈ {1, . . . p}, Ω ,m ) = 0}, the short-memory effect of FIVAR model has a dramatic effect on FWER, Figure 4. This is due to the fact that ρ ,m (j) depends of the short-range memory and thus ρ ,m (j) = 0 does not necessarily implies that Ω ,m = 0. The FWER of (2.5) is still controlled but the inferred graph is far from G * . This problem does not come from the test procedure but from a bias in the statistics estimation. Notice that if the probability to make a false edge detection is not controlled in the first scales, due to short-memory, the number of false positives remains low. Figure 5 illustrates that except for the first scale, only 1 or 2 edges are wrongly put in the rejection set. Finally, when less components are concerned by the VAR part, we observe that much less scales are perturbed by the presence of short-memory.

Additionally, since the procedure is asymptotic, the FWER control is not ensured for small values of n j . This explains why for highest scales (coarsest frequencies) the control is not acquired.

FMRI DATA WITH DEAD RATS

The dataset consists in fMRI recordings on rats. Four rats are scanned dead (rats 1 to 4), and seven rats are scanned alive under aneasthetic (rats 5 to 11). The duration of the scanning is 30 minutes with a time repetition of 0.5 second so that 3, 600 time points are available at the end of experience. After preprocessing as described in, [START_REF] Pawela | Resting-state functional connectivity of the rat brain[END_REF] we extracted 51 time series, each time series being associated with a brain region of rats.

As explained in Section 2, we first estimate the correlation using a wavelet decomposition. Relatively to previous discussions, we will focus on wavelet scale 4 corresponding to the frequency interval [0.06 ; 0.12] Hz. The number of wavelet coefficients available is n 4 = 122. As seen in simulations of Section 4.2 the expected power is not far from 1 with such a number of coefficients (see Figure 6). Examples of correlation matrices obtained are given in Figure 7. Even if dead rats correlation matrices are much closer to identity matrices than correlation matrices on alive rats, some non zero values appear which need a test of significance.

Next, we define the rejection set of (2.5) applying Bonferroni's multiple testing procedure on Fisher statistics. As we have to deal with 50 × 51/2 = 1275 tests for each of the 11 rats, a multiple testing procedure is necessary. Without a multiple testing correction on p-values, the inferred graph is highly biased, as illustrated by the clearly non-empty graph for a dead rat displayed in Figure 8. 4: Empirical probability to make a false edge detection with respect to the scale. The horizontal plain line corresponds to the threshold α = 0.05. 1000 simulations were done with FIVAR models, with (a) the VAR part is null (b) all the components are concerned with a non null VAR part (c) 10 of the 51 components are concerned with a non null VAR part.

(a) q q q q q q q q q 0.00 0.25

q q q q q q q q q q 0.00 0.25 26) scale (number of coefficients) FWER (c) q q q q q q q q q q 0.00 0.25 (a) q q q q q q q q q q 0.00 0.25

q q q q q q q q q q 0.00 0.25

q q q q q q q q q q 0.00 0.25 It is not possible in general to evaluate the accuracy of on real dataset because of the lack of ground truth. However, since there is actually no cerebral activity in a brain of a dead rat, all null hypotheses to be tested are true nulls and the number of rejected hypotheses must be equal to zero in this case. Indeed, we obtain that the number of significant correlations is zero or near to zero for dead rats (see Table 1). We can still observe a few remaining connections in the graphs from very close brain regions. This is due to the fact that the fMRI scanner is introducing bias in the neighborhood parts of the brain.

For alive rats, the proportion of rejected null hypotheses is always higher than 5%. Note that there is a high variability with respect to rats. The highest number of detections is obtained for the rats 5 and 8 with more than 33% of edges. For the others, the number of detections is between 15% and 25%.

We applied a bootstrap procedure on time series using moving windows in time (with overlapping). Results are given in Figure 9. We observe an important difference between alive and dead rats. Dead rats have a mean proportion of rejected null hypotheses equal to zero or really near to zero. The highest mean of rejected Figure 6: Proportion of true edges detected with respect to the scale. 1000 simulations were done with FIVAR models, with (a) the VAR part is null (b) all the components are concerned with a non null VAR part (c) 10 of the 51 components are concerned with a non null VAR part.

(a) q q q q q q q q q q 0.00 0.25

q q q q q q q q q q 0.00 0.25 (c) q q q q q q q q q q 0.00 0.25 Figure 8: Illustration of the necessity to make mutiple correction for a given dead rat.

Without multiple testing q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q With multiple testing q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Some examples of estimated graphs for dead rats and alive rats are given in Figure 10.

CONCLUSION

A whole range of graphical models have been proposed for graph inference, using either correlation or partial correlation as a measure of dependence. Most of these methods relies on assumptions on the structure of the estimated graphs: sparsity for gLasso, conditions on the degree distribution, independence between blocks for stochastic block model. In this paper we propose to use a multiple testing strategy to avoid making such assumptions. Our approach is valid and control the FWER for any graphs' structures. However, quantifying the power of the approach is more complicated and will be the objective of future works.

We considered wavelet-based graph inference by multiple testing. We highlighted some specificity of this approach. The number of available coefficients at each scale is the first limitation due to the asymptotic testing.

Figure 10: Examples of graphs obtained at scale 4 for some alive rats and some dead rats. Bonferroni's correction was used with Fisher statistics.

Alive rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Alive rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Alive rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Dead rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Dead rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Dead rat q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Next the memory properties influences the estimated graph. Long-range dependence implies a specific discrepancy (with respect to graph G * ) which can be corrected if needed. Short-range memory has for consequence that the inferred graph is dependent to the high frequency scales. A perspective is to develop a multiple testing approach to infer the graph G * with a correction of all memory effects.

( W1 )

 W1 The functions φ(•) and ψ(•) are integrable, have compact supports, R φ(t)dt = 1 and ψ 2 (t)dt = 1; (W2) There exists α > 1 such that sup λ∈R | ψ(λ)|(1 + |λ|) α < ∞, i.e. the wavelet is α-regular; (W3) The mother wavelet ψ(•) has M > 1 vanishing moments. (W4) The function k∈Z k φ(• -k) is polynomial with degree for all = 1, . . . , M -1.

  FWER R bonf n,α , P = P   ( ,m)∈H0(P ) {p n, m ≤ α/ν}   .

Figure 2 :

 2 Figure 2: Adjacency matrix of a small-world graph with 51 nodes and 102 edges.

( b ),

 b Next, we will consider a matrix were all components are concerned by a VAR part: Finally, we provide simulations with a matrix Φ such that only a part of the components are concerned with short-memory: with only the upper 10 × 10 block matrix with non-zero values.

Figure

  Figure 4: Empirical probability to make a false edge detection with respect to the scale. The horizontal plain line corresponds to the threshold α = 0.05. 1000 simulations were done with FIVAR models, with (a) the VAR part is null (b) all the components are concerned with a non null VAR part (c) 10 of the 51 components are concerned with a non null VAR part.

Figure 5 :

 5 Figure 5: Empirical proportion of false edge detections in the rejection set with respect to the scale. 1000 simulations were done with FIVAR models, with (a) the VAR part is null (b) all the components are concerned with a non null VAR part (c) 10 of the 51 components are concerned with a non null VAR part.

Figure 7 :

 7 Figure 7: Examples of two correlation matrices between wavelet coefficients obtained at scale 4 for an alive rat (left) and a dead rat (right).

Figure 9 :

 9 Figure 9: Proportion of significant edges with Bonferroni procedure with overlapped time windows bootstrap, for dead rats (left) and alive rats (right).

  correlations with Bonferroni procedure with 50 (highly) overlapped windows bootstrap hypotheses is indeed 4 for 1275 tests. On the contrary for alive rats we observe always more than 10% of null rejected.

Table 1 :

 1 Number of rejected null hypotheses (or estimate edges) obtained by Bonferroni's procedure on Fisher statistics, evaluated at scale 4.

ACKNOWLEDGMENTS

The authors would like to thank Etienne Roquain for very interesting scientific discussions. We also would like to thank Emmanuel Barbier and Guillaume Becq for providing us the data of the resting state fMRI on the rats.


Jonathan Laurent

John Bozek

Marc Briant

Pierre Çarçabal

Denis Cubaynes

Aleksandar Milosavljević

Ralph Püttner

Niloufar Shafizadeh

Marc Simon

Benoît Soep

Aleksan- Dar Milosavljević

Benoit Soep

Gildas Goldsztejn
email: gildas.goldsztejn@universite-paris-saclay.fr

B Soep

Consistent characterization of the electronic ground state of iron(ii) phthalocyanine from valence and core-shell electron spectroscopy

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

Introduction

Since their discovery, metal phthalocyanines (MPc) have been widely studied because of their potential applications in a broad range of domains. From catalysis [START_REF] Leznoff | Phthalocyanines, Properties and Applications[END_REF] , optoelectronics [START_REF] Walzer | [END_REF] , spintronics 3 to photodynamic therapy 4 to name a few. Among them, Iron(II) phthalocyanine (FePc) is viewed as a potential molecular magnet [START_REF] Ishikawa | Phthalocyanine-Based Magnets[END_REF] . However, despite numerous studies in the last six decades, the characterization of its electronic ground state (GS) remains highly controversial.

The Pc is an aromatic macrocycle that ressembles to a porphyrin molecule, except for the substitution of aza-bridge carbons by nitrogen atoms. This therefore leads to the existence of two different types of nitrogen atoms, those surrounding the iron which are directly involved in the chemical bonds with the metal and the peripheric "aza-bridge" nitrogen atoms. In the center of Pc either a metal (Fe in the present case) or two hydrogens (free-based phthalocyanine) are located and the central ring is surrounded by four aromatic rings. The molecule is represented in Fig. 1.

Fig. 1 Iron phthalocyanine posesses a D 4h symmetry and is composed of a porphyrin-like aromatic macrocycle where peripheric nitrogen atoms are in place of aza-bridge carbons. The macrocyle is surrounded by four aromatic rings and a Fe atom is located in its center. The labels "a" to "h" indicate different distances that will be refered to in the discussion .

The GS electronic properties of FePc are linked to the coupling J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

of the Fe 3d orbitals with the π system of the phthalocyanine. This results in three electronic states with different distributions reflected in their spin states. For Fe II ferrous heme, the spins S=0, 1 and 2 are possible due to the population of the d orbitals within the phthalocyanine cage.

One consensus is that the GS of FePc is in a spin-intermediate state with S=1 based on measurements of its magnetic susceptibility as a function of the temperature performed by Dale et al. [START_REF] Dale | [END_REF] and further reinforced by their Mössbauer-effect study 7 where the authors conclude the ground state to be of 3 E g symmetry.

FePc belongs to the D 4h point group where the degeneracy of the five 3d orbitals is lifted into four levels b 1g (x2 y 2 ), a 1g (z 2 ), b 2g (xy) and e g , which is a doublet (xz, yz). With the exception of Thole et al. 8 who described FePc's GS as a quintet/triplet spinmixed state, the literature also agrees on the fact that the b 1g orbital is too high in energy to be occupied. This is essential since the affinity of FePc to ligands is generally determined by the respective b 1g (d x 2 -y 2 ) and a 1g (d z 2 ) populations. However, all possible arrangements for the three remaining orbitals can be found from experimental and/or computational studies.

The 3 A 2g GS with the configuration |a 2 1g b 2 2g e 2 g b 0 1g has been proposed by Stillman and Thomson 9 from magnetic circular dichroism measurements in solution. Note that we adopt a convention where the orbital ordering relates their relative energy. The same GS was predicted by Liao and Scheiner 10 and Sumimoto et al. 11 through density functional theory (DFT) studies. Recently, Greulich et al. 12 described FePc's GS as lying close to the transition between 3 E g and 3 A 2g states after performing X-ray absorption spectroscopy (XAS) and X-ray magnetic circular dichroism measurements at the Fe L-edge in thin films, supported by crystalfield multiplet calculations.

An early reference from Barraclough et al. 13 described the GS as 3 B 2g (i.e. |e 4 g b 1 2g a 1 1g b 0 1g configuration) after single crystal magnetic anisotropy measurements, followed by Brena and collaborators 14 who compared valence photoemission spectra and DFT calculations as well as in reference 15 where the authors performed XAS and XMCD spectra at the Fe L-edge of FePc adsorbed on Ag(001) and compared these with multiplet calculations.

As already mentioned, several articles found the GS to be a 3 E g state 7,[16][17][18][19][20][21] , either in the configuration |b 2 2g e 3 g a 1 1g b 0 1g or |a 2 1g e 3 g b 1 2g b 0 1g . Finally, various references also conclude that the GS is a mixed state of 3 E g -3 B 2g 22-24 or 3 E g -3 A 2g 25-27 symmetries.

Such a diversity of results comes from diverse biases. Experimentally, it was found that for thin films the electronic state is highly dependent on the substrate (see for instance 12,28 ) and, in particular, to charge transfer between the phthalocyanine and the substrate. An alternative approach is to probe MPc in thick films as performed by Bidermane et al. 21 who showed that the C 1s, N 1s and Fe 2p emission spectra were highly similar to those measured in the gas-phase. However, they did not arrive at the same conclusion as Brena et al. 14 who performed valence spectroscopy of FePc in the gas-phase and interpreted the GS on the basis of DFT calculations. Additionally, relying on one simulated spectrum only may not guarantee deriving the proper GS of such complex transition organometallic compound, independently of the theoretical method used. This is mainly due to the fact that the orbitals are quite close in energy and, therefore, a substantial amount of low-lying excited states can give qualitatively satisfying results for the simulation of the spectra.

From the theoretical viewpoint, it has been found that the GS determined by DFT calculations depends highly on the chosen functional and basis set, as in the reference 29 where all simulated valence photoemission spectra are in fair agreement with the experimental one. Regarding the multiplet calculations, which have been shown to be well suited to calculate absorption spectra at the Fe L 2,3 -edge [START_REF] De Groot | Advances in Condensed Matter Science[END_REF] , many different sets of parameters have been reported for the crystal-field. Since these parameters are interdependent and often used to fit a particular data set, it is difficult to find parameters that would fit all the different experimental spectra.

In the present study, we decided to have a comprehensive experimental and theoretical approach to identify the GS of an isolated FePc molecule. We recorded outer-valence, shallow-core and core-shell photoemission spectra in the gas-phase and simulated our results with both DFT and multiplet calculations.

Methods

Experimental

The experiments were performed at the PLEIADES beamline 31,[START_REF] Miron | [END_REF] of the synchrotron radiation facility SOLEIL which is equipped with two undulators (HU256 and HU80) that allow to have access to a tunable photon energy between ≈ 10 -1000 eV. This permits to probe a broad range of edges, from valence spectra close to ionization thresholds -the first IP of the molecule is at ≈ 6 eV-, to C 1s and N 1s electrons (at ≈ 280 and 400 eV, respectively) and Fe 3p, 3s and 2p edges (≈ 60, 90 and 700 eV respectively). This allows recording an extensive set of spectra, in view of reaching a consistent characterization of the electronic GS able to reproduce the corpus of data. The spectra were recorded using the wideangle lens VG-Scienta R4000 electron spectrometer installed at the beamline at a fixed position and the polarization vector of the incoming light was set at the so-called magic angle of 54.7°so that there is negligible influence of the photoelectron angular distribution on the photoionization cross sections. FePc has been graciously provided by Porphychem as a sublimated product. It is mostly free of synthesis pollutants. For the valence electron spectroscopy, data sets were recorded with photon energies between 52 and 110 eV and the electrons with a binding energy in the 6 to 20 eV range were detected with a step width of 20 meV. For the Fe 3p photoelectron spectra, we centered the binding energy of the electrons around 60 eV and used a step width of 100 meV; we varied the photon energy between 80 to 120 eV. We also measured resonant Auger spectra at the Fe 2p edge using photon energies between 703 and 735 eV with a step width of 500 meV. The electrons were recorded in the binding-energy range of 5 to 110 eV with a step width of 200 meV. Such a large binding energy range allows extracting most of the emitted electrons after resonant Auger decays and give, after integration over all binding energies, access to a partial electron yield (PEY) spectrum, which is comparable to absorption spectra. The total experimental resolution, composed of the convolution between the monochromator and spectrometer resolutions, was chosen to represent a good compromise between intensity of the signal and its resolution. In case of the valence spectra the total instrumental resolution amounted ≈ 60 meV, while in case of the Fe 3p XPS spectra it was ≈ 500 meV and for the core resonant Auger spectra where the spectral contributions are inherently large, it was ≈ 600 meV for the photon bandwidth.

The molecules were evaporated at ≈ 700 K under vacuum in a home-made oven designed to fit the Scienta's specifications and both the temperature of the oven and pressure in the interaction chamber were monitored during the data acquiqition to ensure stable target densities. The oven was set on a 3D-axis manipulator and its position was optimized on the signal.

DFT calculations

The DFT calculations were performed with the Gaussian 16 software package [START_REF] Frisch | Gaussian~16 Revision C.01[END_REF] . The structures have been optimized using the hybrid Becke 3-parameter, Lee, Yang and Parr (B3LYP) and the B97D3 functionals. We performed two sets of calculations with and without imposing the D 4h symmetry. We verified that the electronic structures as well as the c 2 coefficients for the molecular orbitals and their energies were almost identical in both cases. Therefore, to enable the comparison with both our multiplet calculations and the literature, we retained the results of the imposed-symmetry calculations in the following. We used different basis sets in order to evaluate their effect on the optimized structure and found that the B3LYP functional was robust with the choice of the basis set, i.e. whatever the chosen basis set we found i) very comparable geometric structure, ii) that the triplet state is the most stable electronic structure, followed by the quintuplet (≈ 0.9 eV higher) and the singlet (≈ 1.3 eV higher in energy). In the following, we present the results with the correlation-consistent polarized triple-zeta (cc-pVTZ) basis set used for all atoms. To simulate the valence photoelectron (PE) spectra, we calculated the partial and total density of states (DOS) by using the optimized structure. The result using the B3LYP functional matches slightly better the experimental spectrum than with the B97 functional, as can be seen in the figure S2 of the supplementary material. The partial DOS (PDOS) were calculated using the keyword pop=full iop(3/33=1,3/36=-1) in Gaussian on a single point calculation on the optimized structure and the percentages of each groups of atoms to the molecular orbitals were calculated using the Mulliken population analysis and the software GaussSum 34 . To simulate the PE spectra, we followed the procedure described by Brena et al. 14 and used the Gelius model 35 which consists in weighing each atomic orbital of the PDOS by the atomic subshell photoionization cross sections; for these cross sections we used the theoretical values from Yeh and Lindau [START_REF] Yeh | Atomic Data and Nuclear Data Tables[END_REF][START_REF] Yeh | Atomic Calculation of Photoionization Cross-Sections and Asymmetry Parameters[END_REF] . Finally, due to the hybrid character of the sp orbitals of C and N atoms, we decided to weight their contributions in the PDOS by the population given in the DFT calculation instead of separating them into 2s and 2p contributions as in Ref 14 . To simulate the PE spectra, we convoluted the theoretical peaks using Gaussian functions with four different full width at half maximum (FWHM) according to their binding energies, i.e. FWHM = 135 meV for energies below 6.5 eV, FWHM = 300 meV for energies between 6.5 and 8.2 eV, FWHM = 500 meV for 8.2 ≤ E ≤ 9.5 eV and FWHM = 1.5 eV above 9.5 eV. This allows to take into account that at higher energies the contributions of 2-holes 1-electron final states corresponding to simultaneous ionization and excitation become more important.

Multiplet calculations

In Fig. 2, we described the crystal-field parameters in a singleelectron picture. In a free ion, all 3d orbitals have the same energy. In contrast to this, in the O h (octahedral) symmetry there is a splitting of the d-level into two groups of orbitals: e g , which is a doublet comprising the z 2 and x 2 -y 2 orbitals, and t 2g , which is a triplet constituted by the xy, xz and yz orbitals. In this symmetry, the only important parameter is 10Dq which relates to the energy difference between the e g and t 2g orbitals. When the symmetry is further lowered to a tetragonal symmetry D 4h , the orbitals are split into 4 levels, namely b 1g (x 2 -y 2 ), a 1g (z 2 ), b 2g (xy) and e g , which is doubly degenerated (xz, yz). To describe the crystal field, two more parameters are needed, the so-called Ds and Dt. The four levels are linked by a set of formulas [START_REF] De Groot | Advances in Condensed Matter Science[END_REF][START_REF] Zimmermann | [END_REF][START_REF] Reddy | Electronic (Absorption) Spectra of 3d Transition Metal Complexes[END_REF] :

E b 1g = +6Dq + 2Ds -Dt (x 2 -y 2 ) E a 1g = +6Dq -2Ds -6Dt (z 2 ) E e g = -4Dq -Ds + 4Dt (xz, yz) E b 2g = -4Dq + 2Ds -Dt (xy)
However, these energy shifts are only correct for a single electron in the d-shell because the formulas do not take into account the coupling between the electrons, i.e. the situation is much more complex when the four levels are occupied by several electrons.

To perform the multiplet calculations and simulate the PEY and PE spectra, we used the Charge Transfer Multiplet program for X-ray Absorption Spectroscopy (CTM4XAS version 5.5) [START_REF] Stavitski | [END_REF] and for the determination of the ground state, we used the version 1.9 of the CTM4DOC 41 (differential orbital covalency) software. Depending on the chosen parameters, the CTM4XAS program allows to take into account the core-hole induced effects on the electronic structure, which are neglected when considering the core/inner-shell excitation or ionization as a single electron process. These effects include its potential, the spin-orbit coupling, the core-hole induced charge transfer effect, the core and valence holes exchange and multipole interactions [START_REF] Stavitski | [END_REF] . To further detail these multiplet calculations, the Hamiltonian is described as a sum of three contributions, namely i) a free ion Hamiltonian regrouping the kinetic energy of the electrons, the electron-nuclear, the electron-electron and the spin-orbit interactions, ii) a crystal field Hamiltonian and iii) an hybridization Hamiltonian. Note that the contributions ii) and iii) characterize the ligand field experienced by the centered metal atom. Once the Hamiltonian is described, the wavefunction can be calculated and further pro- jected on a basis consisting of all possible configurations, each affiliated with a coefficient α i (eigenvalues). The numbers given on page 6 in front of each configuration are the α 2 i values for the 6 dominant configurations (∑ α 2 i = 1, if all possible configurations are taken into account).

A large number of parameters' combinations and approaches can be found in the literature to match various spectra. Here we chose the strategy to use the same parameters to fit all spectra (L-edge PEY and PE spectra and M-edge PE spectrum). We reduced the Slater radial integrals to 68% of their Hartree-Fock value which is in agreement with the Refs. [START_REF] Dale | [END_REF]12,24,25,27 This lowering is related to the fact that the bonding between Fe and its neighbours is partly covalent and partly ionic. In general, when this factor is comprised between 0.7 and 0.9 the bond is considered predominantly with an ionic character, while when it is between 0.5 and 0.6, it has a covalent character. Here, we are in an intermediate situation adding complexity to the characterization of the GS. The optical parameters found as a best fit for our experimental spectra are 10Dq = 2.66 eV, Ds = 0.625 eV and Dt = 0.234 eV which are in relatively close agreement with the references 12,18,23,25 . To find the accurate charge transfer parameters, we fixed the charge transfer energy to ∆ = 1.63 eV. This value comes from the interpretation of a ligand-to-metal charge transfer band in the visible absorption spectrum of gas-phase FePc 42 and is consistent with the interpretation of the absorption spectrum of FePc in dichlorobenzene 9,11 .

To determine the values of U pd , the core-hole potential, and U dd , the Hubbard value, we fitted our X-ray PE spectra. It resulted in the charge transfer parameters: ∆ = 1.3 eV, U dd = 5 eV and U pd = -2 eV. ∆ describes the energy needed to transfer one electron from the ligand to the metal atom and U dd as well as U pd represent the Coulomb energy between two 3d electrons as well as a 3d electron and a 2p core electron, respectively. Finally, we used T(b1) = T(a1) = 2 eV and T(b2) = T(e) = 1 eV, where T are the metal-ligand hopping integrals.

Although the discussion of the influence of the substrate on the electron spectra, as compared to the gas-phase measurements, is out of the scope of the present study, we tried to fit our data with the parameters given in the works of Greulich et al. 12 and Stepanow et al. 23 that deal with films of FePc on GeS(100) and Au(111) substrates respectively. We observed that for the first one cited, the agreement with our spectra is only slightly worse than with the parameters used in the present study, which goes along with the choice of Greulich et al. for this substrate as to avoid interactions with FePc. There is, however, a clear difference when using the parameters used in Stepanow et al. where, as they noted, the spin of the substrate couples with the Fe ion.

Results and discussion

Valence spectroscopy and DFT calculations

The experimental spectrum taken at hν=55.64 eV is shown in Fig. 3 (red dots), along with the simulated DOS (black line) and the individual atomic intensities are plotted below as sticks. The simulated spectrum has been normalized over the most intense experimental point. The other features, marked by a star, are dominated by water contribution (blue stars), other peaks (red stars) are either due to other contaminants of the sample or to simultaneous ionized and excited states which are not accounted for in our simulation. We show in Fig. S1 of the supplementary material, the experimental spectrum of the background taken in the same experimental conditions, i.e. oven on but without the FePc. The good agreement between simulated and experimental spectra is an indication that the structure of FePc is properly reproduced by the DFT calculations. Indeed, in Table 1, we summarized the values of our calculations for the distances labeled from "a" to "h" in Fig. 1 and compared them to experimental values obtained with X-ray diffraction 43 . The values presented in Sumimoto et al. 11 and Brena et al. 14 and ∆r max is the maximum deviation between calculated and experimental values, both in Å.

Table 1 Summary of the distances "a" to "h" in Fig. 1 given for the optimized structures of the present work, of Sumimoto et al. 11 and of Brena et al. 14 , as well as the experimental values from Ref. 43 . δ r is the average deviation to the experimental values (see text) and ∆r max is the maximum deviation between a single calculated and experimental values, both in Å.

Expt. 43 Although the three DFT calculations presented in table 1 show comparable agreements with the experimental values and geometric structures being very close to each other, they predict three different descriptions of the electronic GS. We find the GS to be 3 E g (|b 2 2g e 3 g a 1 1g b 0 1g ), Sumimoto et al. 11 found it to be a 3 A 2g (|b 2 2g e 2 g a 2 1g b 0 1g ), while Brena et al. 14 gave a 3 B 2g (|b 1 2g e 4 g a 1 1g b 0 1g ) GS. The origin of these differences in the ground-state configurations is probably due to the fact that the energy levels of the three lowest d orbitals (according to the present calculation, the b 2g orbital is the lowest in energy, but it is separated by only 0.4 eV and 0.44 eV from the e g and the a 1g orbital, respectively, while the b 1g orbital is 2.65 eV higher in energy) are close to each other so that different theoretical approaches may result in different sequences for the three lowest levels.

Moreover, in agreement with Sumimoto et al. 11 , our DFT calculations give similar optimized geometrical structures for the singlet and quintuplet states or even by using a different functional, namely the B97D3 (see table S1 in the supplemental material), as those given in table 1. Furthermore, all calculations give a reasonably good agreement with the experiment spectrum when simulating their DOS spectra (see Figs. S2-S3 of the supplementary material). This is also in agreement with Marom et al. 29 who performed DFT calculations with many different functional and basis sets. In detail, they found for MgPc and FePc that the various simulated DOS spectra match well the valence spectra, although the obtained GS were different.

To summarize, we found that the simulated DOS spectrum related to the DFT calculation done with the B3LYP functional gives the best agreement with the valence spectra, compared to the B97D3 functional. It also describes the triplet state as being more stable than the quintuplet and singlet electronic ones. However, since the agreement of DOS spectra stemming from all the abovementioned DFT calculations with our experimental spectrum is relatively good, it did not seem advisable to us to simply rely on this result and chose to confirm it with the use of shallow-core and core-shell spectroscopies.

Core-hole spectroscopy and multiplet calculations

To determine the electronic structure of the d-orbitals, it is more efficient to probe it directly at the Fe site. This can be done by core electron emission.

In Fig. 4, we plotted the PEY measured around the L 2,3 -edges as red dots. The black solid line is the simulated spectrum using the charge transfer and crystal field parameters described above. To compare with the experimental spectrum, we used two Lorentzian functions of 0.1 (L 3 ) and 0.6 eV (L 2 ) and convolved by a Gaussian function of 0.6 eV to simulate the experimental resolution. The result of our simulation is in fairly good agreement with the experimental spectrum. In particular, it reproduces well the shoulder at low binding energy (≈706 eV) and the spin-orbit energy difference between L 2 and L 3 is also well reproduced. The main discrepancy resides in the tail region which may be due to vibrational progressions (see for instance Ref. 44 ). We did not measure the Fe 2p photoelectrons spectrum of FePc, however in Fig. 5, the experimental spectrum from Bidermane et. al. 21 , reproduced with permission, for thick films (red dots) was used to test the quality of our parameters on these data as well. The result of our simulation is displayed on the same graph as a black solid line and the intensities of each contribution is shown below as black sticks. To match the experimental spectrum, we convoluted the spectrum with two Lorentzian functions of 0.1 and 0.6 eV and a Gaussian function of 0.3 eV. The agreement with the PE spectrum is very good and reproduces well the shoulder at ≈707 eV. In contrast to the cited reference, where the authors included a final state different from the initial one in order to better fit their results, we chose to find a consistent set of parameters that suits all shallow-core and core-shell spectra. In Figure . 6 we present the experimental PE spectrum of the shallow-core level Fe 3p as red dots and the simulated spectrum by a black solid line. Here, we also find a very good match between our CTM calculations, which is convoluted with Lorentzian functions of 0.1 and 0.4 eV of FWHM and a Gaussian function of 0.2 eV, and the experimental spectrum. The Fe 3p photoelectron spectrum, which to the best of our knowledge has never been published before, turned out to be particularly important to determine the parameters that describe the ligand field splitting. where the last two lines correspond to the effect of charge transfer. Therefore the ground state is predominantly of 3 E g symmetry, as was found from our DFT calculation. In both approaches we found the same orbital order, i.e. the b 2g orbital is the lowest in energy, followed by e g , a 1g and b 1g . We find charge transfer to be responsible for ∼26% of the GS, which is a much larger contribution than that reported in Ref. 12 .

As in the article of Carlotto et al. 45 , our calculations show that the energy levels of the three different symmetries 3 E g , 3 B 2g and 3 A 2g are close to each other. However, we found differences in the sequence of states and obtained splittings. In detail, in the Ref. 45 , the GS was found to be the 3 A 2g state separated to the 3 E g and 3 B 2g states by 0.05 and 0.06 eV respectively. In the present study the 3 E g is lower in energy by 0.14 and 0.37 eV than the 3 B 2g and 3 A 2g states, respectively.

To further corroborate our assessment from the CTM calculations, we also simulated the Fe 2p XAS and XPS spectra for the states 3 B 2g and 3 A 2g . The simulated XAS and XPS spectra of the 3 A 2g state do not agree with the experimental spectra. For the 3 B 2g state, the situation is more complex since in this case the simulated XAS and Fe 3p XPs spectra of the 3 B 2g state are very similar to the 3 E g state, which we identified as ground state. In contrast to this, the 2p XPS spectra (Fig. S4 of the supplementary material) shows clear differences in the energy range from 707 to 708 eV, with the simulated spectrum of the 3 E g state being better in line with the experimental spectrum. Secondly, the intensity ratio between the L2 and L3 edges is poorly reproduced for the 3 B 2g state. In summary, only the simulated spectra of the 3 E g state match all experimental data. This is an additional evidence that one needs a consistent analysis of different data sets in order to identify the electronic GS of such complex systems.

The configuration given by the CTM and the DFT calculations, respectively |b 1.73 2g e 3.13 g a 1.25 1g b 0.15 1g and |b 1.97 2g e 2.93 g a 1.01 1g b 0.62 1g , compare well to each other. In case of the DFT calculation, the result shows that 6.53 electrons are located on the Fe 3d orbitals, this is also in relative agreement with the 6.26 electrons found by the CTM calculation. However, a difference resides in the fact that in case of the DFT calculation almost all of the charge is transferred from the lone pairs of the N atoms go the b 1g orbital, while it is predicted by the CTM calculation to be evenly distributed between a 1g and b 1g orbitals. This difference between both results probably lies in the fact that for the latter, we define the metal atom in a certain degree of oxidation (+2 in the present study) and the molecular environment is considered only as a perturbation to the free ion's Hamiltonian.

Conclusion

We have performed PE and PEY measurements at various edges of the FePc molecule and complemented them by DFT and CTM calculations in order to simulate these spectra and, in turn, identify the GS of the isolated molecule. Our approach consisted in finding a consistent description able to reproduce the corpus of available data. We observed that valence photoelectron spectroscopy is less sensitive to the exact GS, since different structures and spin-states can provide simulated spectra that all reasonably match the experimental one. This result can readily be understood by the fact that with valence photoelectron spectroscopy, the entire molecule is probed. In contrast to this, inner-shell PE spectra that are direct probes of the centered-metal atom, which is a crucial step to identify the GS. By comparing the PE spectra at the 2p and 3p edges as well as the PEY around the 2p resonances and by fitting them with the same crystal field and charge transfer parameters, we were able to identify that the GS of FePc to be of 3 E g symmetry and find that the charge transfer between the Fe atom and its N atoms neighbours has to be taken into account.

Fig. 2

 2 Fig. 2 Representation of the symmetry effect on the 3d orbitals. For a spherical symmetry the 5 3d orbitals have the same energy, in an octahedral symmetry (O h ), there is a degeneracy splitting into two different levels e g and t 2g separated in energy by ≈ 10 Dq. When the symmetry is further lowered to the D 4h point group, the d orbitals are futher separated into 4 different levels named b 1g , a 1g , e g and b 2g . In this case, two distorsion parameters Ds and Dt describe the splitting between the levels.

Fig. 3

 3 Fig. 3 Valence photoelectron spectrum measured at hν = 55.64 eV (red dots) along with a simulation using the density of states which is represented as a black solid line. The intensities of each individual molecular orbital is indicated by the stick lines. The blue stars indicate contributions of the background which is mostly due to water contributions. Red stars correspond to feature that may belong either to other contaminants of the sample or to simultaneous ionized and excited states not accounted for by the DOS simulation.

Fig. 4

 4 Fig. 4 The partial electron yield (red dots) in the photon energy range from 703 to 724.5 eV. The parameters used to plot the simulated XAS spectrum (black solid line) are detailed in the main text. The intensities of the simulated spectrum are indicated with black sticks.

Fig. 5

 5 Fig. 5 XPS spectrum at the Fe L 2,3 edges measured (red dots) by Bidermane et al. 21 and reproduced here to test our CTM parameters. The result of the calculation is plotted as a black solid line and the intensities of each contribution are plotted as stick lines underneath.

Fig. 6

 6 Fig. 6 Fe 3p XPS spectrum measured at hν = 100 eV. The red points correspond to the experimental spectrum, while the black solid line is the calculated spectrum. The individual theoretical intensities are plotted as stick lines below the graph.

  are also displayed in the table and show a comparable agreement with our values. The two last lines of the table are the average deviation to the experimental values calculated as δ r =

	i="h" ∑ i="a"	|r i,exp -r i,DFT | 8

| 1-8

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

G.G. wants to thank Prof. Dr. Frank de Groot for invaluable help on the usage of the CTM softwares. We are thankful to Dr Habermeyer from Porphychem for providing a sample of sublimated FePc.

Notes and references

Author Contributions

Conflicts of interest

There are no conflicts to declare.