Dynamical Kernels for Optical Excitations

Juliette Authier 1 and Pierre-François Loos 1, a) Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

We discuss the physical properties and accuracy of three distinct dynamical (i.e., frequency-dependent) kernels for the computation of optical excitations within linear response theory: i) an a priori built kernel inspired by the dressed time-dependent density-functional theory (TDDFT) kernel proposed by Maitra and coworkers [J. Chem. Phys. 120, 5932 (2004)], ii) the dynamical kernel stemming from the Bethe-Salpeter equation (BSE) formalism derived originally by Strinati [Riv. Nuovo Cimento 11, 1-86 (1988)], and iii) the second-order BSE kernel derived by Yang and coworkers [J. Chem. Phys. 139, 154109 (2013)]. The principal take-home message of the present paper is that dynamical kernels can provide, thanks to their frequency-dependent nature, additional excitations that can be associated to higher-order excitations (such as the infamous double excitations), an unappreciated feature of dynamical quantities. We also analyze, for each kernel, the appearance of spurious excitations originating from the approximate nature of the kernels, as first evidenced by Romaniello et

I. LINEAR RESPONSE THEORY

Linear response theory is a powerful approach that allows to directly access the optical excitations ω S of a given electronic system (such as a molecule) and their corresponding oscillator strengths [extracted from their eigenvectors (X S Y S ) ] via the response of the system to a weak electromagnetic field. [1][START_REF] Casida | Time-dependent density functional response theory for molecules[END_REF][START_REF] Petersilka | [END_REF] From a practical point of view, these quantities are obtained by solving non-linear, frequency-dependent Casida-like equations in the space of single excitations and de-excitations 2

R σ (ω S ) C σ (ω S ) -C σ (-ω S ) * -R σ (-ω S ) * • X σ S Y σ S = ω S X σ S Y σ S (1)
where the explicit expressions of the resonant and coupling blocks, R σ (ω) and C σ (ω), depend on the spin manifold (σ = ↑↓ for singlets and σ = ↑↑ for triplets) and the level of approximation that one employs. Neglecting the coupling block [i.e., C σ (ω) = 0] between the resonant and anti-resonants parts, R σ (ω) and -R σ (-ω) * , is known as the Tamm-Dancoff approximation (TDA). In the absence of symmetry breaking, 4 the nonlinear eigenvalue problem defined in Eq. ( 1) has particle-hole symmetry which means that it is invariant via the transformation ω → -ω. Therefore, without loss of generality, we will restrict our analysis to positive frequencies.

In the one-electron basis of (real) spatial orbitals {φ p (r)}, we will assume that the elements of the matrices defined in Eq. ( 1) have the following generic forms: 4

R σ ia, jb (ω) = ( a -i )δ i j δ ab + f Hxc,σ ia, jb (ω) (2a) C σ ia, jb (ω) = f Hxc,σ ia,b j (ω) (2b)
where δ pq is the Kronecker delta, p is the one-electron (or quasiparticle) energy associated with φ p (r), and

f Hxc,σ ia, jb (ω) = φ i (r)φ a (r) f Hxc,σ (ω)φ j (r )φ b (r )drdr (3)
a) Electronic mail: loos@irsamc.ups-tlse.fr

Here, i and j are occupied orbitals, a and b are unoccupied orbitals, and p, q, r, and s indicate arbitrary orbitals. In Eq. ( 3),

f Hxc,σ (ω) = f Hx,σ + f c,σ (ω) (4)
is the (spin-resolved) Hartree-exchange-correlation (Hxc) dynamical kernel. In the case of a spin-independent kernel, we will drop the superscript σ. As readily seen from Eq. ( 4), only the correlation (c) part of the kernel is frequency dependent and, in a wave function context, the static Hartree-exchange (Hx) matrix elements read

f Hx,σ ia, jb = 2σ(ia| jb) -(ib| ja) (5)
where σ = 1 or 0 for singlet and triplet excited states (respectively), and

(pq|rs) = φ p (r)φ q (r) 1 |r -r | φ r (r )φ s (r )drdr (6)
are the usual two-electron integrals. 5 The launchpad of the present study is that, thanks to its non-linear nature stemming from its frequency dependence, a dynamical kernel potentially generates more than just single excitations. Unless otherwise stated, atomic units are used, and we assume real quantities throughout this manuscript.

II. THE CONCEPT OF DYNAMICAL QUANTITIES

As a chemist, it is maybe difficult to understand the concept of dynamical properties, the motivation behind their introduction, and their actual usefulness. Here, we will try to give a pedagogical example showing the importance of dynamical quantities and their main purposes. [6][7][START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF] To do so, let us consider the usual chemical scenario where one wants to get the optical excitations of a given system. In most cases, this can be done by solving a set of linear equations of the form

A • c = ω c ( 7
)
where ω is one of the optical excitation energies of interest and c its transition vector. If we assume that the matrix A is arXiv:2008.13143v2 [physics.chem-ph] 1 Sep 2020 diagonalizable and of size N × N, the linear set of equations (7) yields N excitation energies. However, in practice, N might be (very) large (e.g., equal to the total number of single and double excitations generated from a reference Slater determinant), and it might therefore be practically useful to recast this system as two smaller coupled systems, such that

A 1 b b A 2 • c 1 c 2 = ω c 1 c 2 (8)
where the blocks A 1 and A 2 , of sizes

N 1 × N 1 and N 2 × N 2 (with N 1 + N 2 = N)
, can be associated with, for example, the single and double excitations of the system. This decomposition technique is often called Löwdin partitioning in the literature. [START_REF] Löwdin | [END_REF] Solving separately each row of the system (8) and assuming that ω1 -A 2 is invertible, we get

A 1 • c 1 + b • c 2 = ω c 1 (9a) c 2 = (ω 1 -A 2 ) -1 • b • c 1 (9b)
Substituting Eq. (9b) into Eq. (9a) yields the following effective non-linear, frequency-dependent operator

Ã1 (ω) • c 1 = ω c 1 ( 10
)
with

Ã1 (ω) = A 1 + b • (ω 1 -A 2 ) -1 • b (11)
which has, by construction, exactly the same solutions as the linear system (7) but a smaller dimension. For example, an operator Ã1 (ω) built in the single-excitation basis can potentially provide excitation energies for double excitations thanks to its frequency-dependent nature, the information from the double excitations being "folded" into Ã1 (ω) via Eq. (9b). [START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF] Note that this exact decomposition does not alter, in any case, the values of the excitation energies.

How have we been able to reduce the dimension of the problem while keeping the same number of solutions? To do so, we have transformed a linear operator A into a non-linear operator Ã1 (ω) by making it frequency dependent. In other words, we have sacrificed the linearity of the system in order to obtain a new, non-linear system of equations of smaller dimension [see Eq. (10)]. This procedure converting degrees of freedom into frequency or energy dependence is very general and can be applied in various contexts. [10][11][12][13][14][15][16][17][18][19][20][21] Thanks to its non-linearity, Eq. ( 10) can produce more solutions than its actual dimension. However, because there is no free lunch, this non-linear system is obviously harder to solve than its corresponding linear analog given by Eq. (7). Nonetheless, approximations can be now applied to Eq. (10) in order to solve it efficiently. For example, assuming that A 2 is a diagonal matrix is of common practice (see, for example, Ref. 18 and references therein).

Another of these approximations is the so-called static approximation, where one sets the frequency to a particular value. For example, as commonly done within the Bethe-Salpeter equation (BSE) formalism of many-body perturbation theory (MBPT), 22 Ã1 (ω) = Ã1 ≡ Ã1 (ω = 0). In such a way, the operator Ã1 is made linear again by removing its frequencydependent nature. A similar example in the context of timedependent density-functional theory (TDDFT) 23 is provided by the ubiquitous adiabatic approximation, 24 which neglects all memory effects by making static the exchange-correlation (xc) kernel (i.e., frequency independent). [START_REF] Maitra | Memory: History , initial-state dependence , and doubleexcitations[END_REF][START_REF] Maitra | [END_REF][27] These approximations come with a heavy price as the number of solutions provided by the system of equations ( 10) has now been reduced from N to N 1 . Coming back to our example, in the static (or adiabatic) approximation, the operator Ã1 built in the single-excitation basis cannot provide double excitations anymore, and the N 1 excitation energies are associated with single excitations. All additional solutions associated with higher excitations have been forever lost. In the next section, we illustrate these concepts and the various tricks that can be used to recover some of these dynamical effects starting from the static eigenproblem.

III. DYNAMICAL KERNELS A. Exact Hamiltonian

Let us consider a two-level quantum system made of two electrons in its singlet ground state (i.e., the lowest orbital is doubly occupied). 6 We will label these two orbitals, φ v and φ c , as valence (v) and conduction (c) orbitals with respective one-electron Hartree-Fock (HF) energies v and c . In a more quantum chemical language, these correspond to the HOMO and LUMO orbitals (respectively). The ground state |0 has a one-electron configuration |vv , while the doubly-excited state |D has a configuration |cc . There is then only one single excitation possible which corresponds to the transition v → c with different spin-flip configurations. As usual, this produces a singlet singly-excited state |S = (|vc + |cv )/ √ 2, and a triplet singly-excited state |T = (|vc -|cv )/ √ 2. [START_REF] Szabo | Modern quantum chemistry[END_REF] For the singlet manifold, the exact Hamiltonian in the basis of these (spin-adapted) configuration state functions reads 29

H ↑↓ =           0| Ĥ|0 0| Ĥ|S 0| Ĥ|D S | Ĥ|0 S | Ĥ|S S | Ĥ|D D| Ĥ|0 D| Ĥ|S D| Ĥ|D           (12)
with

0| Ĥ|0 = 2 v -(vv|vv) = E HF (13a) S | Ĥ -E HF |S = ∆ + (vc|cv) -(vv|cc) (13b) D| Ĥ -E HF |D = 2∆ + (vv|vv) + (cc|cc) + 2(vc|cv) -4(vv|cc) (13c) 0| Ĥ|S = 0 (13d) S | Ĥ|D = √ 2[(vc|cc) -(cv|vv)] (13e) 0| Ĥ|D = (vc|cv) (13f)
and ∆ = cv . The energy of the only triplet state is simply

T | Ĥ|T = E HF + ∆ -(vv|cc).
For the sake of illustration, we will use the same numerical examples throughout this study, and consider the singlet ground state of i) the H 2 molecule (R H-H = 1.4 bohr) in the STO-3G basis, ii) the HeH + molecule (R He-H = 1.4632 bohr) in the STO-3G basis, and iii) the He atom in Pople's 6-31G basis set. These three systems provide prototypical examples of valence, charge-transfer, and Rydberg excitations, respectively, and will be employed to quantity the performance of the various methods considered in the present study for each type of excited states. The numerical values of the various quantities defined above are gathered in Table I for each system.

The exact values of the singlet single and double excitations, ω ↑↓ 1 and ω ↑↓ 2 , and the triplet single excitation, ω ↑↑ 1 , are reported, for example, in Table II. We are going to use these as reference for the remaining of this study.

B. Maitra's dynamical kernel

The kernel proposed by Maitra and coworkers 30,31 in the context of dressed TDDFT (D-TDDFT) corresponds to an ad hoc many-body theory correction to TDDFT. More specifically, D-TDDFT adds to the static kernel a frequency-dependent part by reverse-engineering the exact Hamiltonian: a single and double excitations, assumed to be strongly coupled, are isolated from among the spectrum and added manually to the static kernel. The very same idea was taking further by Huix-Rotllant, Casida and coworkers, 32 and tested on a large set of molecules. Here, we start instead from a HF reference. The static problem corresponds then to the time-dependent HF (TDHF) Hamiltonian, while in the TDA, it reduces to configuration interaction with singles (CIS). 4 For the two-level model, the reverse-engineering process of the exact Hamiltonian (12) yields

f c,↑↓ M (ω) = | S | Ĥ|D | 2 ω -( D| Ĥ|D -0| Ĥ|0 ) ( 14
)
while f c,↑↑ M (ω) = 0. The expression ( 14) can be easily obtained by folding the double excitation onto the single excitation, as explained in Sec. II. It is clear that one must know a priori the structure of the Hamiltonian to construct such dynamical kernel, and this obviously hampers its applicability to realistic photochemical systems where it is sometimes hard to get a clear picture of the interplay between excited states. [33][34][35] For the two-level model, the non-linear equations defined in Eq. ( 1) provides the following effective Hamiltonian

H σ D-TDHF (ω) = R σ M (ω) C σ M (ω) -C σ M (-ω) -R σ M (-ω) (15)
with

R σ M (ω) = ∆ + 2σ(vc|vc) -(vc|vc) + f c,σ M (ω) (16a) C σ M (ω) = 2σ(vc|cv) -(vv|cc) + f c,σ M (ω) (16b)
yielding, for our three two-electron systems, the excitation energies reported in Table II when diagonalized. The TDHF Hamiltonian is obtained from Eq. ( 15) by setting f c,σ M (ω) = 0 in Eqs. (16a) and (16b). In Fig. 1, we plot det[H(ω) -ω1] as a function of ω for both the singlet (black and gray) and triplet (orange) manifolds in HeH + . (Very similar curves are obtained

�������� (����) �������� (�-����) �������� (����) ��� ��� ��� ��� ��� ��� ��� -� -� � � � � � �� ω ���[�(ω)-ω�] FIG. 1. det[H(ω) -ω1]
as a function of ω (in hartree) for both the singlet (gray and black) and triplet (orange) manifolds of HeH + . The static TDHF Hamiltonian (dashed) and dynamic D-TDHF Hamiltonian (solid) are considered.

for He.) The roots of det[H(ω) -ω1] indicate the excitation energies. Because, there is nothing to dress for the triplet state, the TDHF and D-TDHF triplet excitation energies are equal.

Although not particularly accurate for the single excitations, Maitra's dynamical kernel allows to access the double excitation with good accuracy and provides exactly the right number of solutions (two singlets and one triplet). Note that this correlation kernel is known to work best in the weak correlation regime (which is the case here) where the true excitations have a clear single and double excitation character, 36,37 but it is not intended to explore strongly correlated systems. 38 Its accuracy for the single excitations could be certainly improved in a density-functional theory context. However, this is not the point of the present investigation.

Table II also reports the slightly improved (thanks to error compensation) CIS and D-CIS excitation energies. In particular, single excitations are greatly improved without altering the accuracy of the double excitation. Graphically, the curves obtained for CIS and D-CIS are extremely similar to the ones of TDHF and D-TDHF depicted in Fig. 1.

In the case of H 2 in a minimal basis, because S | Ĥ|D = 0, there is no dynamical correction for both singlets and triplets, and one cannot access the double excitation with Maitra's kernel. It would be, of course, a different story in a larger basis set where the coupling between singles and doubles would be non-zero.

C. Dynamical BSE kernel

As mentioned in Sec. II, most of BSE calculations performed nowadays are done within the static approximation. [START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF][39][40][41] However, following Strinati's footsteps, 22,42,43 several groups have explored this formalism beyond the static approximation by retaining the dynamical nature of the screened Coulomb potential W 6,7,17,44 or via a perturbative approach. [45][46][47][48][START_REF] Loos | Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation[END_REF] Based on the very same two-level model that we employ here, Romaniello et al. clearly evidenced that one can genuinely access additional excitations by solving the non-linear, frequency-dependent BSE eigenvalue problem. For this particular system, they showed that a BSE kernel based on the random-phase approximation (RPA) produces indeed double excitations but also unphysical excitations, 6 attributed to the self-screening problem. 50 This issue was resolved in the subsequent work of Sangalli et al. 7 via the design of a diagrammatic number-conserving approach based on the folding of the second-RPA Hamiltonian. 51 Thanks to a careful diagrammatic analysis of the dynamical kernel, they showed that their approach produces the correct number of optically active poles, and this was further illustrated by computing the polarizability of two unsaturated hydrocarbon chains (C 8 H 2 and C 4 H 6 ). Very recently, Loos and Blase have applied the dynamical correction to the BSE beyond the plasmon-pole approximation within a renormalized first-order perturbative treatment, [START_REF] Loos | Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation[END_REF] generalizing the work of Rolhfing and coworkers on biological chromophores 46,47 and dicyanovinylsubstituted oligothiophenes. 48 They compiled a comprehensive set of vertical transitions in prototypical molecules, providing benchmark data and showing that the dynamical correction can be sizable (especially for n → π * and π → π * excitations) and improves the static BSE excitations considerably. [START_REF] Loos | Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation[END_REF] Let us stress that, in all these studies, the TDA is applied to the dynamical correction (i.e., only the diagonal part of the BSE Hamiltonian is made frequency-dependent) and we shall do the same here.

Within the so-called GW approximation of MBPT, 8,39,52-54 one can easily compute the quasiparticle energies associated with the valence and conduction orbitals. [55][56][57][58] Assuming that W has been calculated at the random-phase approximation (RPA) level and within the TDA, the expression of the GW quasiparticle energy is simply 59

GW p = p + Z GW p Σ GW p ( p ) (17)
where p = v or c,

Σ GW p (ω) = 2(pv|vc) 2 ω -v + Ω + 2(pc|cv) 2 ω -c -Ω (18)
is the correlation part of the self-energy Σ, and

Z GW p =          1 - ∂Σ GW p (ω) ∂ω ω= p          -1 (19)
is the renormalization factor (or spectral weight). In Eq. ( 18), Ω = ∆ + 2(vc|cv) is the sole (singlet) RPA excitation energy of the system, with ∆ GW = GW c -GW v . One can now build the dynamical BSE (dBSE) Hamiltonian 6,22

H σ dBSE (ω) = R σ dBSE (ω) C σ dBSE -C σ dBSE -R σ dBSE (-ω) (20)
with

R σ dBSE (ω) = ∆ GW + 2σ(vc|cv) -(vv|cc) -W c R (ω) (21a) C σ dBSE = 2σ(vc|cv) -(vc|cv) -W c C (ω = 0) (21b)
and where

W c R (ω) = 4(vv|vc)(vc|cc) ω -Ω -∆ GW (22a) W c C (ω) = 4(vc|cv) 2 ω -Ω (22b)
are the elements of the correlation part of the dynamicallyscreened Coulomb potential for the resonant and coupling blocks of the dBSE Hamiltonian, respectively. Note that, in this case, the correlation kernel is spin blind. Within the usual static approximation, the BSE Hamiltonian is simply

H σ BSE = R σ BSE C σ BSE -C σ BSE -R σ BSE ( 23
)
with ). It is worth mentioning that, around ω = ω dBSE,σ 1 , the slope of the curves depicted in Fig. 2 is small, while the other solution, ω dBSE,σ 2 , stems from a pole and consequently the slope is very large around this frequency value. This makes this latter solution quite hard to locate with a method like Newton-Raphson (for example). Let us highlight the fact that, unlike in Ref. 49 where dynamical effects have been shown to produce a systematic red-shift of the static excitations, here we observe both blue-and red-shifted transitions (see values in square brackets in Table III).

R σ BSE = ∆ GW + 2σ(vc|vc) -(vv|cc) -W R (ω = ∆ GW ) (24a) C σ BSE = C σ dBSE (24b) �������� (���) �������� (����) �������� (���) �������� (����) � � � � � -�� � �� �� �� ω ���[�(ω)-ω�] FIG. 2. det[H(ω) -ω1]
In the static approximation, only one solution per spin manifold is obtained by diagonalizing H σ BSE (see Fig. 2 and Table III). Therefore, the static BSE Hamiltonian misses the (singlet) double excitation (as it should), and it shows that the physical single excitation stemming from the dBSE Hamiltonian is indeed the lowest in energy for each spin manifold, i.e., ω dBSE,↑↓ 1 and ω dBSE,↑↑ 1 . This can be further verified by switching off gradually the electron-electron interaction as one would do in the adiabatic connection formalism. 60,61 Enforcing the TDA, which corresponds to neglecting the coupling term between the resonant and anti-resonant parts of the dBSE Hamiltonian (20), does not change the situation in terms of spurious solutions: there is still one spurious excitation in the triplet manifold (ω BSE,↑↑ 2 ), and the two solutions for the singlet manifold which corresponds to the single and double excitations. However, it does increase significantly the static excitations while the magnitude of the dynamical corrections is not altered by the TDA. The (static) BSE triplets are notably too low in energy as compared to the exact results and the TDA is able to partly reduce this error, a situation analogous in larger systems. [62][63][64] Another way to access dynamical effects while staying in the static framework is to use perturbation theory, [45][46][47][48][START_REF] Loos | Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation[END_REF] a scheme we label as perturbative BSE (pBSE). To do so, one must decompose the dBSE Hamiltonian into a (zeroth-order) static part and a dynamical perturbation, such that

H σ dBSE (ω) = H σ BSE H (0) pBSE + H σ dBSE (ω) -H σ BSE H (1) pBSE (ω) (26)
Thanks to (renormalized) first-order perturbation theory, [START_REF] Loos | Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation[END_REF] one gets

ω pBSE,σ 1 = ω BSE,σ 1 + Z pBSE 1 X 1 Y 1 • H σ dBSE (ω = ω BSE,σ 1 ) -H σ BSE • X 1 Y 1 (27)
where

H σ BSE • X 1 Y 1 = ω BSE,σ 1 X 1 Y 1 (28)
and the renormalization factor is

Z pBSE 1 =        1 - X 1 Y 1 • ∂H σ dBSE (ω) ∂ω ω=ω BSE,σ 1 • X 1 Y 1        -1 (29)
This corresponds to a dynamical perturbative correction to the static excitations. The perturbatively-corrected values are also reported in Table III, and it shows that this scheme is very effective at reproducing the dynamical values for the single excitations. Because the value of Z 1 is always quite close to unity in the present systems (evidencing that the perturbative expansion behaves nicely), one could have anticipated the fact that the first-order correction is a good estimate of the non-perturbative result. However, because the perturbative treatment is ultimately static, one cannot access double excitations with such a scheme.

For H 2 , there is no dynamical corrections at the BSE, pBSE or dBSE levels. Indeed, as (vv|vc) = (vc|cc) = 0 (see Table I), we have

W c R (ω) = 0 [see Eq. (22a)].
The lack of frequency dependence of the kernel means that one cannot estimate the energy of the doubly-excited state of H 2 .

D. Second-order BSE kernel

The third and final dynamical kernel that we consider here is the second-order BSE (BSE2) kernel derived by Yang and collaborators in the TDA, 65 and by Rebolini and Toulouse in a range-separated context 66,[START_REF] Rebolini | Range-Separated Density-Functional Theory for Molecular Excitation Energies[END_REF] (see also Refs. 44, 68, and 69). Note that a beyond-TDA BSE2 kernel was also derived in Ref. 66, but was not tested. In a nutshell, the BSE2 scheme applies second-order perturbation theory to optical excitations within the Green's function framework by taking the functional derivative of the second-order self-energy Σ GF2 with respect to the one-body Green's function. Because Σ GF2 is a proper functional derivative, it was claimed in Ref. 65 that BSE2 does not produce spurious excitations. However, as we will show below, this is not always true. Like BSE requires GW quasiparticle energies, BSE2 requires the second-order Green's function (GF2) quasiparticle energies, [START_REF] Szabo | Modern quantum chemistry[END_REF] which are defined as follows:

GF2 p = p + Z GF2 p Σ GF2 p ( p ) (30)
where the second-order self-energy is

Σ GF2 p (ω) = (pv|vc) 2 ω -v + c -v + (pc|cv) 2 ω -c -( c -v ) (31)
and

Z GF2 p =          1 - ∂Σ GF2 p (ω) ∂ω ω= p          -1 (32)
The expression of the GF2 self-energy (31) can be easily obtained from its GW counterpart (18) via the substitution Ω → cv and by dividing the numerator by a factor two. This shows that there is no screening within GF2, but that second-order exchange is properly taken into account. 65,70 The static Hamiltonian of BSE2 is just the usual TDHF Hamiltonian where one substitutes the HF orbital energies by the GF2 quasiparticle energies, i.e.,

H σ BSE2 = R σ BSE2 C σ BSE2 -C σ BSE2 -R σ BSE2 ( 33
)
with

R σ BSE2 = ∆ GF2 + 2σ(vc|vc) -(vv|cc) (34a) C σ BSE2 = 2σ(vc|vc) -(vc|cv) (34b)
To avoid any confusion with the results of Sec. III B and for notational consistency with Sec. III C, we have labeled this static Hamiltonian as BSE2.

The correlation part of the dynamical kernel for BSE2 is a bit cumbersome [65][66][START_REF] Rebolini | Range-Separated Density-Functional Theory for Molecular Excitation Energies[END_REF] but it simplifies greatly in the case of the present model to yield

H σ dBSE2 = H σ BSE2 + R c,σ dBSE2 (ω) C c,σ dBSE2 -C c,σ dBSE2 -R c,σ BSE2 (-ω) (35)
with

R c,↑↓ dBSE2 (ω) = - 4(cv|vv)(vc|cc) -(vc|cc) 2 -(cv|vv) 2 ω -2∆ GF2 (36a) C c,↑↓ dBSE2 = 4(vc|cv) 2 -(cc|cc)(vc|cv) -(vv|vv)(vc|cv) 2∆ GF2 (36b) and R c,↑↑ dBSE2 (ω) = - (vc|cc) 2 + (cv|vv) 2 ω -2∆ GF2 (37a) C c,↑↑ dBSE2 = (cc|cc)(vc|cv) + (vv|vv)(vc|cv) 2∆ GF2 (37b)
As mentioned in Ref. 66, the BSE2 kernel has some similarities with the second-order polarization-propagator approximation 1,71 (SOPPA) and second RPA kernels. 7,32,51,[START_REF] Huix-Rotllant | Improved correlation kernels for linear-response timedependent density-functional theory[END_REF] Unlike the dBSE Hamiltonian [see Eq. ( 20)], the BSE2 dynamical kernel is spin aware with distinct expressions for singlets and triplets. [START_REF] Rebolini | Range-Separated Density-Functional Theory for Molecular Excitation Energies[END_REF] Like in dBSE, dBSE2 generates the right number of excitations for the singlet manifold (see Fig. 3). However, one spurious triplet excitation clearly remains. Numerical results for the two-level models are reported in Table IV with the usual approximations and perturbative treatments. In the case of BSE2, the perturbative partitioning (pBSE2) is simply

H σ dBSE2 (ω) = H σ BSE2 H (0) pBSE2 + H σ dBSE2 (ω) -H σ BSE2 H (1) pBSE2 (ω) (38)
�������� (�����) � � � � � � � -�� � �� �� �� �� ω ���[�(ω)-ω�]
FIG. 3. det[H(ω) -ω1] as a function of ω (in hartree) for both the singlet (gray and black) and triplet (orange and red) manifolds of HeH + . The static BSE2 Hamiltonian (dashed) and dynamic dBSE2 Hamiltonian (solid) are considered.

As compared to dBSE, dBSE2 produces much larger dynamical corrections to the static excitation energies (see values in square brackets in Table IV) probably due to the poorer quality of its static reference (TDHF or CIS). Similarly to what has been observed in Sec. III B, the TDA vertical excitations are slightly more accurate due to error compensations. Note also that the perturbative treatment is a remarkably good approximation to the dynamical scheme for single excitations (except for H 2 , see below), especially in the TDA. This justifies the use of the perturbative treatment in Refs. 65 and 66. Overall, the accuracy of dBSE and dBSE2 are comparable for single excitations although their behavior is quite different (see Tables III andIV). For the double excitation, dBSE2 yields a slightly better energy, yet still in quite poor agreement with the exact value.

Again, the case of H 2 is a bit peculiar as the perturbative treatment (pBSE2) does not provide any dynamical corrections, while its dynamical version (dBSE2) does yield sizable corrections originating from the coupling term C c,σ dBSE2 which is nonzero in the case of dBSE2. Although frequency-independent, this additional term makes the singlet and triplet excitation energies very accurate. However, one cannot access the double excitation.

IV. TAKE-HOME MESSAGES

The take-home message of the present paper is that dynamical kernels have much more to give that one would think. In more scientific terms, dynamical kernels can provide, thanks to their frequency-dependent nature, additional excitations that can be associated to higher-order excitations (such as the infamous double excitations), an unappreciated feature of dynamical quantities. However, they sometimes give too much, and generate spurious excitations, i.e., excitation which does not correspond to any physical excited state. The appearance of these factitious excitations is due to the approximate nature of the dynamical kernel. Moreover, because of the non-linear character of the linear response problem when one employs a dynamical kernel, it is computationally more involved to access these extra excitations.

Using a simple two-model system, we have explored the physics of three dynamical kernels: i) a kernel based on the dressed TDDFT method introduced by Maitra and coworkers, 30 ii) the dynamical kernel from the BSE formalism derived by Strinati in his hallmark 1988 paper, 22 as well as the secondorder BSE kernel derived by Zhang et al., 65 and Rebolini and Toulouse. 66,[START_REF] Rebolini | Range-Separated Density-Functional Theory for Molecular Excitation Energies[END_REF] Prototypical examples of valence, chargetransfer, and Rydberg excited states have been considered. From these, we have observed that, overall, the dynamical correction usually improves the static excitation energies, and that, although one can access double excitations, the accuracy of the BSE and BSE2 kernels for double excitations is rather average. If one has no interest in double excitations, a perturbative treatment is an excellent alternative to a non-linear resolution of the dynamical equations.

We hope that the present contribution will foster new developments around dynamical kernels for optical excitations, in particular to access double excitations in molecular systems.

  al. [J. Chem. Phys. 130, 044108 (2009)]. Using a simple two-level model, prototypical examples of valence, charge-transfer, and Rydberg excited states are considered.

TABLE I .

 I Numerical values (in hartree) of the valence and conduction orbital energies, v and c , and two-electron integrals in the orbital basis for various two-level systems. HF/STO-3G -1.632 802 -0.172 484 +0.943 099 +0.752 526 +0.660 254 +0.145 397 -0.172 968 +0.037 282 He HF/6-31G -0.914 127 +1.399 859 +1.026 907 +0.766 363 +0.858 133 +0.227 670 +0.316 490 +0.255 554 TABLE II. Singlet and triplet excitation energies (in eV) for various levels of theory and two-level systems. The magnitude of the dynamical correction is reported in square brackets.

	System Method	v	c	(vv|vv)	(cc|cc)	(vv|cc)	(vc|cv)	(vv|vc)	(vc|cc)
	H 2 HeH +	HF/STO-3G -0.578 203 +0.670 268 +0.674 594 +0.697 495 +0.663 564 +0.181 258	0	0
			Method
	System Excitation CIS TDHF D-CIS	D-TDHF	Exact
	H 2 HeH + He	ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1	25.78 25.30 25.78[+0.00] 25.30[+0.00] 26.34 44.04 15.92 15.13 15.92[+0.00] 15.13[+0.00] 16.48 29.68 29.42 27.75[-1.93] 27.64[-1.78] 28.05 63.59 63.52 64.09 21.77 21.41 21.77[+0.00] 21.41[+0.00] 22.03 52.01 51.64 51.87[-0.14] 51.52[-0.12] 52.29 93.85 93.84 94.66 39.62 39.13 39.62[+0.00] 39.13[+0.00] 40.18

  as a function of ω (in hartree) for both the singlet (gray and black) and triplet (orange and red) manifolds of HeH + . The static BSE Hamiltonian (dashed) and dynamic dBSE Hamiltonian (solid) are considered.

	It can be easily shown that solving the secular equation
	det H σ dBSE (ω) -ω1 = 0	(25)

yields 2 solutions per spin manifold (except for H 2 where only one root is observed, see below), as shown in Fig.

2

for the case of HeH + . Their numerical values are reported in Table

III

alongside other variants discussed below. These numbers evidence that dBSE reproduces qualitatively well the singlet and triplet single excitations, but quite badly the double excitation which is off by several eV. As mentioned in Ref.

6

, spurious solutions appear due to the approximate nature of the dBSE kernel. Indeed, diagonalizing the exact Hamiltonian (12) produces only two singlet solutions corresponding to the singlyand doubly-excited states, and one triplet state (see Sec. III A). Therefore, there is the right number of singlet solutions but there is one spurious solution for the triplet manifold (ω dBSE,↑↑ 2

TABLE III .

 III Singlet and triplet BSE excitation energies (in eV) for various levels of theory and two-level systems. The magnitude of the dynamical correction is reported in square brackets.

						Method
	System	Excitation	BSE	pBSE	dBSE	BSE(TDA)	pBSE(TDA)	dBSE(TDA)	Exact
	H 2 HeH + He	ω ↑↓ 1 ω ↑↑ 1 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↑ 2 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↑ 2	26.06 16.94 28.56 20.96 52.46 40.50	26.06[+0.00] 16.94[+0.00] 28.63[+0.07] 21.07[+0.11] 52.12[-0.34] 39.80[-0.70]	26.06[+0.00] 16.94[+0.00] 28.63[+0.07] 87.47 21.07[+0.11] 87.43 52.11[-0.35] 133.38 39.79[-0.71] 133.75	27.02 17.16 29.04 21.13 53.10 40.71	27.02[+0.00] 17.16[+0.00] 29.11[+0.07] 21.24[+0.11] 52.79[-0.31] 40.02[-0.69]	27.02[+0.00] 17.16[+0.00] 29.11[+0.07] 87.47 21.24[+0.11] 87.43 52.79[-0.31] 133.37 40.02[-0.69] 133.75	26.34 16.48 28.05 64.09 22.03 52.29 94.66 40.18

TABLE IV .

 IV Singlet and triplet BSE2 excitation energies (in eV) for various levels of theory and two-level systems. The magnitude of the dynamical correction is reported in square brackets.

						Method
	System	Excitation	BSE2	pBSE2	dBSE2	BSE2(TDA)	pBSE2(TDA)	dBSE2(TDA)	Exact
	H 2 HeH + He	ω ↑↓ 1 ω ↑↑ 1 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↑ 2 ω ↑↓ 1 ω ↑↓ 2 ω ↑↑ 1 ω ↑↑ 2	26.03 15.88 29.23 21.22 50.31 37.80	26.03[+0.00] 15.88[+0.00] 28.40[-0.83] 21.63[+0.41] 51.96[+1.64] 39.26[+1.46]	26.24[+0.21] 16.47[+0.59] 28.56[-0.67] 79.94 21.93[+0.71] 78.70 52.10[+1.79] 121.67 39.59[+1.79] 121.85	26.49 16.63 29.50 21.59 50.69 38.30	26.49[+0.00] 16.63[+0.00] 28.66[-0.84] 21.99[+0.40] 52.34[+1.65] 39.77[+1.47]	26.49[+0.00] 16.63[+0.00] 28.66[-0.84] 79.94 21.99[+0.40] 78.70 52.34[+1.65] 121.66 39.77[+1.47] 121.84	26.34 16.48 28.05 64.09 22.03 52.29 94.66 40.18
		�������� (����)
		�������� (�����)
		�������� (����)

ACKNOWLEDGMENTS

We would like to thank Xavier Blase, Elisa Rebolini, Pina Romaniello, Arjan Berger, Miquel Huix-Rotllant and Julien Toulouse for insightful discussions on dynamical kernels. PFL thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 863481) for financial support.

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article


Lilian Moumaneix

Jenifer Guerrero Parra

Sébastien Fontana

François Lapicque

Claire Hérold

Dr Claire Herold
email: claire.herold@univ-lorraine.fr33372742537

Investigation and mechanism proposal of solvothermal reaction between sodium and 1-(2-Hydroxyethyl)Piperidine as the first step towards nitrogen-doped graphenic foam synthesis

In view to synthetizing 3-dimensional nitrogen-doped graphenic materials, which could be used as oxygen reduction catalyst in membrane fuel cells, a solvothermal-based route has been successfully carried out. However, the solvothermal reaction between metallic sodium and 1-(2-hydroxyethyl)piperidine (HEP), the latter being the source of both carbon and nitrogen, is still little understood. The present work was aimed at investigating the solvothermal process, under different conditions of temperature, pressure and different amounts of sodium. Use of in-situ mass spectroscopy during the three-day reaction revealed the early formation of dihydrogen, as well as carbon oxides, methane and ammonia, in addition to fragment of ethylpiperidine alkoxide (EP-ONa). XRD measurements evidenced the formation of sodium-based compounds, e.g. hydride, carbonate, hydroxide, cyanide. Interestingly, Raman spectroscopy revealed the significant presence of large aromatic molecules as well as a sp 2 carbon network, early precursor of graphene. Analysis of the overpressures and reaction yields led to suggest that the primary compound of HEP reaction with sodium is a large sp 2 carbon-based network entrapping numerous sodium-based molecules as well as a volatile liquid phase. The suggested reactional mechanism provides information to better tailor the solvothermal products, whose pyrolysis at 850 °C lead to very high specific area nitrogendoped carbon materials.

A Introduction

Supercritical fluids (SCF) provide distinctive experimental media as their properties strongly depend on both internal pressure and temperature. Halfway between gas and liquids states, supercritical fluids give great opportunities to have control over many physical properties, e.g. reactivity, density, polarity or viscosity [1]. SCF possess liquid-like densities (0.2-0.9 g.cm -3 ), gas-like viscosities (1-9 x 10 -5 Pa.s -1 ) and intermediary diffusivities (2-7 x 10 -8 m 2 .s -1 ) [START_REF] Erkey | Supercritical Fluid Science and Technology[END_REF]. Technologies based on SCF are daily used in many industrial applications such as the extraction of numerous organic molecules from agricultural products e.g. aroma from plants, caffeine from coffee, cold sterilization, impregnation at heart or plastic packaging treatment. In materials science, phenomena made possible by SCF media can be of a great interest for synthesis or post-treatment methods for numerous materials. With regard to carbon materials, literature contains several examples of SCF applications to process nanocarbon materials including graphenic materials. Nowadays, graphene is widely studied in numerous scientific fields thanks to its electric, thermal, chemical and mechanical properties. Numerous operations can be achieved on carbon materials in SCF conditions e.g. graphite exfoliation [START_REF] Pu | [END_REF]4], synthesis of graphene composites [5], elaboration of three-dimensional graphenic materials [6][7][8] or deposition of metal nanoparticles on graphene materials [9].

In the last decades, various research works had focused on obtaining three-dimensional graphenic materials such as crumpled graphene sheets, graphene frameworks or graphene foams through a solvothermal-based process. A solvothermal reaction is defined as a reaction in a closed system between reagents in the presence of a solvent (which can be part of the reagents) at a temperature higher than the boiling point of the solvent. If both pressure and temperature are above those of the solvent critical point (below respectively), the solvothermal reaction is said in supercritical conditions (subcritical conditions respectively). In addition to allowing reaction occurrence at moderate temperatures, closed system processes also help with the reduction in noxious emissions in environment [10]. A heat treatment is often needed after the reaction to convert the carbon-rich solvothermal product to well-organized graphenic materials. This route of synthesis allows to recover substantial quantity of graphenic material and possesses a good potential for scalability towards industrial processes [10]. In 2004, Kuang et al. [11] reported the elaboration of crumpled graphene nanosheets by the reaction between tetrachloromethane (CCl4) and potassium in an autoclave heated at 60 °C for 10 hours. Eftekhari et al. [12] adapted this synthesis route to produce curly graphene sheets by adding cobalt chloride which acts as a catalyst to transform graphene sheets into carbon nanotubes. In parallel, several works dealt with oxidation reactions of sodium in contact with an alcohol, e.g. methanol, ethanol or butan-1-ol [13][14][15][16][17]. In their study, Choucair et al. [13] synthesized crumpled graphene sheets using ethanol and sodium. The solvothermal reaction was carried out at 220 °C for 72 h, followed by a flash pyrolysis. They proposed a mechanism based on the encapsulation of free-ethanol by saturated metal alkoxide in a clathrate-like structure during the solvothermal reaction to explain the obtained structure. However, complementary investigations would be required to confirm their original theory. Speyer et al. [14][15][16] considered nearly the This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins same process but applied 100 bar upon nitrogen introduction into the reactor to increase furthermore the pressure during the reaction. Graphenic samples showed a noticeable improvement in their crystallinity. These publications and other works not cited here, display the possibility to transform small organic molecules into large sp 2 3-dimensional graphenic networks, using a convenient, highly adaptable technique. Recently, nitrogen-doped carbon materials have attracted great attention due to their potentials as critical metal-free catalysts [18][19][20], their semi-conducting properties or their capacity to provide good anchors for metal nanoparticles nucleation [21,22]. Synthesized by arc-discharge-based synthesis [23], chemical vapor deposition (CVD) [24], thermal treatment [25,26], plasma treatment [27] or metallothermic reduction [28] high quality N-doped carbon materials can also be prepared through a direct solvothermal-based synthesis [29][30][31][32]. However, only a few examples of solvothermal synthesis of Ndoped graphenic materials have been reported so far. Tetrachloromethane (CCl4) is a common precursor used when mixed with lithium nitride (Li3N) [29], Li3N and cyanuric chloride (NCCl)3 [30] or with pyrrole (C4H5N) and potassium [31]. However, the formation of chlorine-derived gases could be avoided by changing the reagents and the reaction conditions. In a previous study [33], a solvothermal-based process using a mixture of cyclohexanol, ethanolamine and sodium, was developed for the synthesis of good quality nitrogen-doped graphenic foam, whose extensive characterization revealed high surface area, good crystallinity and promising catalytic properties. Moreover, this study allowed determination of the optimal temperature of pyrolysis at 850 °C, yielding N-doped graphenic foams with the best compromise between homogeneity, purity and crystallinity, along with appreciable synthesis yields. The present paper has focused on the solvothermal reaction, which has been very little investigated by itself up to now. For this pioneering work, it was preferred to consider the reaction of a single molecule being the sources of C and N atoms: 1-(2hydroxyethyl)piperidine (C7H15NO), subsequently called HEP, and metallic sodium, here combined in view to form, after a subsequent pyrolysis treatment, high-quality nitrogen-doped graphenic foams. The variation of the solvothermal product with the reaction parameters has been followed using complementary techniques for better understanding the reactional mechanism involved here and more generally for carbon materials synthesis through solvothermal-based processes.

B Materials and methods

B.1. Synthesis

A mixture of 0.17 mol of 1-(2-hydroxyethyl)piperidine (ACROS organics, 99%) was put in contact with a certain amount of metallic sodium (Merck Millipore) under inert atmosphere (nitrogen Air Liquide, Alphagaz 2) in a 0.5 L Parr autoclave. 1-(2hydroxyethyl)piperidine was used as received. Sodium rods were washed in petroleum ether to remove the paraffin oil, then scraped before use to remove the passivation layer formed even under protective atmosphere. Theoretical critical conditions of 1-(2-hydroxyethyl)piperidine have been estimated around 385 ± 29 °C and 37 ± 7 bar using computational models [34][35][36]. Thus, the solvothermal reaction has been carried out at 300 °C, 350 °C and 380 °C, further increase in temperature being avoided for technological reasons. The solvothermal reaction pressure has been varied between 100 bar and 200 bar, a range which have been proved in previous works suitable for synthesis of high-quality graphenic materials [33]. Once the autoclave sealed, nitrogen (Air Liquide, Alphagaz 2) was introduced so that the internal pressure at the reaction temperature could reach its desired level at 100 bar or 200 bar. The reactor was then heated up at 8 °C.min -1 to the targeted temperature and the internal pressure reached the desired value. The mixture was let to react for 72 h, then allowed to cool down to room temperature. The internal pressure inside the autoclave was measured and the solvothermal product was recovered and stored under nitrogen atmosphere. Because solvothermal products are poorly stable in air, characterizations were thus carried out as fast as possible, keeping the product within inert atmosphere as long as possible.

B.2. Characterizations

Raman spectroscopy investigations were conducted with a Renishaw inVia Quontor equipped with a 532 nm laser operating in the spectral range 100-3600 cm -1 in synchroscan mode. HEP was simply poured inside a watch glass, the laser beam being focused just under the surface of the liquid. Solvothermal products were kept inside sealed silica quartz cells to avoid any contact with air. X-Ray diffraction (XRD) were performed on a Bruker D8 Advance diffractometer using Mo(Kα1) radiation, λ = 0.70930 Å. Analysed samples were introduced in 1.5 mm wide Lindemann glass capillaries under nitrogen (nitrogen Air Liquide, Alphagaz 2) atmosphere. In situ mass spectrometry was carried out using a solvothermal reactor directly connected to a mass spectrometer Pfeiffer Vacuum Prisma. The pressure in the entire device was controlled at the desired level using a set of valves. Measured intensities were normalized to the N2 peak (channel 28) in order to perform proper comparisons throughout the analysis without influence from slight changes in the measured gas flow.

Please do not adjust margins

Please do not adjust margins Please do not adjust margins Please do not adjust margins Molecular masses in the range 0-120 g.mol -1 were investigated, no peaks being identifiable above this limit. Compounds identification was achieved using the National Institutes of Standards and Technology (NIST) library as a reference.

C Results and discussion

C.1. Pressure and reaction yield investigations

Regarding the stoichiometry generally considered for the reaction between sodium and an alcohol (equation 1), the reaction yield and released pressure can be calculated. EP stands for 1-ethylpiperidine.

HEP + Na → EP-ONa + ½H2 Equation 1Because sodium is always introduced in the same or higher quantity than HEP, the latter one will be considered as being limiting in equation 1. Considering a total reaction, the final molar amount of EP-ONa is then equal to the initial molar HEP amount. The amount of H2 produced is equal to half the one of HEP initially introduced, i.e. 0.087 mol here. The reaction yield was defined as the ratio between the amount of solvothermal product recovered over the maximum amount of EP-ONa that could be obtained, taking into account the excess of sodium introduced. The theoretical pressure released by the formation of H2 has been estimated using the ideal gas law. Considering the volume of the solvothermal reactor, the theoretical partial pressure of hydrogen at room temperature (293 K) was calculated as 4.24 bar (4.24 x 10 5 Pa). The overpressures displayed in table 1 and figure 1 result from the difference between the pressure of N2 injected at room temperature before the heating, and the final pressure, measured at room temperature, after the reaction. In this study, the impact of three parameters has been investigated, namely solvothermal reaction temperature and pressure as well as the molar ratio between sodium and HEP. Table 1 gathers the operating conditions considered and elements on the products obtained. The uncertainty on the pressure measures was estimated at ± 1 bar. As full recovery of the solvothermal product from the reactor can be troublesome because of its high viscosity, the loss in solvothermal product was estimated to be lower than 1 g, leading to an uncertainty in the reaction yield near ± 5 wt.%. Please note that runs SP-T/350 and SP-R/1.5 as well as SP-T/380 and SP-P/200 are strictly identical but have been given different names for the sake of clarity. The reaction yields and overpressures are reported in figure 1. Theoretical H2 partial pressure as expected from equation 1 has been added to evidence the difference in pressure between theory and experiments. The obtained results are interesting in many aspects. Firstly, the overpressure before opening the solvothermal reactor at ambient temperature, was found in all cases larger than the theoretical hydrogen overpressure, expressing the fact that the postulated reaction given in equation 1 would not be the only chemical process: other reactions leading to the formation of gas species would then occur. Secondly, the overpressure value can be correlated to that of reaction operating conditions. An increase in the solvothermal reaction temperature causes an elevation of the overpressure, from 4 ± 1 bar to 12 ± 1 bar at 300 °C and 380 °C respectively, with a slight decrease between 350 °C and 380 °C (cf. figure 1 (a)). This increase in pressure could be due to the promoting of gas creating reactions at higher temperatures compared to the lowest one. On the contrary, an increase in the Na/HEP ratio at constant temperature results in a reduction in the overpressure, from 14 ± 1 bar to 11 ± 1 bar for a Na/HEP ratio of 1.0 and 2.0 respectively (cf. figure 1 could reduce the overpressure by reacting with the gas phase species. A decrease in overpressure is also observed with a larger solvothermal reaction pressure, from 18 ± 1 bar to 12 ± 1 bar at 100 bar and 200 bar respectively (cf. figure 1 (c)). The less significant gas phase formation can simply be explained by Le Chatelier's principle and has been previously reported in literature as a compression effect [37], being mostly favourable to the creation of denser structures. Finally, the reaction yields always ranged from 60 wt.% to 70 wt.%. However, in a sealed reactor every loss of matter must result in the formation of another product, being in solid, liquid or gaseous state. Different hypotheses have been expressed and tested to explain the observed loss of matter. A direct relation between the overpressure and the low yields has been firstly considered. This overpressure could be due to the release of light molecules, e.g. H2, CH4, H2O, coming from the degradation of the solvothermal product. However, the low molar masses of such molecules do not suffice to explain the observed mass loss. Formation of heavier molecules could also be imagined, such as oxidation products CO, NO, CO2 or NO2. In this case, a relationship between the overpressure and the mass loss could be established. However, the initial amount of oxygen atoms inside the reactor is insufficient to produce such amounts of carbon or nitrogen oxides. Formation of volatile organic compounds could moreover be envisaged but it is unlikely that their volatility would correspond to overpressure at room temperature over 10 bar. Regarding the obtained results and calculations, the direct correlation between the gas and solid phases does not suffice to fully explain the losses of mass in the solid phase. After cooling down and measurement of the gas pressure, solvothermal products were recovered and calculations of reaction yields were achieved however, only after opening of the reactor and subsequent release of part of the gas phase: a loss of matter due to the sudden decompression is likely to happen and may explain the yields usually below 70 wt.%. This phenomenon will be further discussed later in this paper. Solvothermal products can also be visually characterized by their physical aspects and more especially their oily aspect. As this introduced parameter is difficult to quantify, also because of their instability in air, a qualitative appreciation has been given using three different denominations. The first one is "very oily" and has been applied for samples presenting a liquid phase, often presenting a brownish colour. The second one is "oily" and has been given to samples containing a unique solid phase but susceptible to leave oily traces on the glass flask. The last one is "dry" and correspond to powdery samples leaving no oily trace in their flasks. Photographs illustrating these denominations are presented in figure 2. The aspect of solvothermal product is given in table 1, along with the reaction conditions. Solvothermal products pass from very oily to oily and then dry by increasing the temperature or the Na/HEP ratio. At 300 °C, liquid compounds are present in high quantities in comparison to the solids formed. By increasing the temperature to 350 °C, this liquid phase appears partly replaced by a solid phase, being in little significant amount upon further increase to 380 °C. This phenomenon could be either explained with thermal decomposition or thermally favoured conversion reactions. A slight change of colour, from brown to perfectly black, has also been observed with the increase in temperature, indicating that the solvothermal product obtained at 380 °C is likely carbon-richer than the one obtained at 300 °C, more especially with a larger density of conjugated double bonds. Similar observations have been made upon increase in Na/HEP ratio, from 1.0 to 2.0, except for the colour remaining dark brownish for the powdery sample SP-R/2.0. The less marked oily aspect can be attributed to the higher number of possible reactions with Na, which is present in greater excess as the Na/HEP ratio increases. Thus, the liquid phase is more likely to react with the available sodium to give solid sodium compounds such as NaOH or Na2CO3. No visible change has been observed with the increase in solvothermal reaction pressure, this parameter being possibly less critical in the range investigated than solvothermal reaction temperature or Na/HEP ratio. Pressures above 100 bar do not seem to be determinant to form dry compounds. However, experiments not presented here showed that a 100-bar pressure leads, after a pyrolysis step, to poor quality N-doped graphenic foams.

(b)). Addition of sodium

C.2. Gas phase characterization: Mass spectrometry investigations

In situ mass spectrometry has been used to investigate the gas phase throughout a typical solvothermal reaction, carried out at 380°C and 200 bar with a Na/HEP ratio of 1.5. These conditions have been chosen as they are known to produce the graphenic material with the highest quality after thermal post-treatment (data not shown here). Channels of interest have been plotted over reaction time in figure 3. The log scale in the plots has been used only because of the large range in order of magnitude of the signal. The in-situ analysis was preceded by a blank test, without reagents introduced in the reactor, which allowed to take into account the impurities in the glove bag atmosphere, such as traces of O2, CO2 and H2O. Values from the blank analysis have been added in the form of orange dotted lines in figure 3 (a-e). Analysis of collected data shows a fast change in the gaseous phase composition within the first hours of the reaction. A large increase in the intensity of channel 2 (cf. figure 3 (a)), related to H2, coupled with the appearance of peaks at m/z = 98, 99, 84 and 85 (cf. figure 3 (f)), associated with solvothermal alkoxide EP-ONa fragments, demonstrate the occurrence of the reaction between Na and HEP (equation 1). Moreover, since the Please do not adjust margins Please do not adjust margins intensities of both H2 and EP-ONa channels tend to stabilize after a few hours, it is likely that this reaction be completed within this time lapse or reach a steady level. Signal intensity of channels 15, 17 and 44 (cf. figure 3 (b,c,e)) related to CH4, NH3 and CO2 respectively, increased rapidly in the first hours of the synthesis. This shows the occurrence of secondary reactions, which are likely responsible for the increase in the final gas pressure discussed in the previous part. CH4 intensity increased around 40-fold in the first hours, then followed by a steady but far slower increase with time: the intensity was three times higher at the end of the run (4300 min) than at 220 min. The steady rise in CH4 could be related to the thermal degradation of the solid phase, rich in carbon atoms, throughout the synthesis. Additionally, methanation reactions are highly suspected to occur, as high pressures and the temperature domain investigated are known to favour CO and CO2 conversion to CH4 [38]. Methane formation usually occurs through a set of four main gas-phase reactions, involving reverse water-gas shift (equation 2), CO methanation (equation 3), reverse dry reforming (equation 4) and CO2 methanation (equation 5) [39]. Formation of CO could not be followed since its molecular weight (m/z = 28) is the same as that for N2, a major component in the gas phase.

CO2 + H2 ⇌ CO + H2O Equation 2 CO + 3H2 ⇌ CH4 + H2O Equation 3 2CO + 2H2 ⇌ CH4 + CO2 Equation 4 CO2 + 4H2 ⇌ CH4 + 2H2O Equation 5 Channel 18, related to H2O, has also been plotted (cf. figure 3 (d)). No evident trend but a slight decrease with time could be observed. However, water was detected in the same order of magnitude as in the blank measurement: it therefore appears not certainly representative of real changes in the gaseous phase. However, water formation followed by its fast reaction with Na excess -resulting in low steady H2O content -could not be excluded. Besides, thermodynamic calculations made on the simple (N2-H2 system) showed that a significant part of the introduced nitrogen can react with dihydrogen to form ammonia, from approx. 3.6 bar at 300 °C to 2.7 bar at 380 °C, under 200 bar of total pressure [START_REF] Kaiser | Techniques de l'ingénieur[END_REF]. Although conducted considering a simpler Please do not adjust margins Please do not adjust margins Please do not adjust margins Please do not adjust margins

. X-Ray diffraction investigations

In order to further understand what species or bond could be formed during the solvothermal reaction, X-ray diffraction (XRD) has been carried out on sample SP-T/350 (same as SP-R/1.5). This sample being median in terms of solvothermal reaction temperature and Na/HEP ratio, every possible solidstate phases are expected to be present on its diffractogramshown in figure 4. The major part of diffraction peaks has been indexed as sodium compounds such as sodium carbonate Na2CO3, sodium hydride NaH, sodium cyanide NaCN and sodium hydride cyanamide Na5H(CN2)3. No peaks related to a crystalline carbon phase can be observed due to the lack of organization in the carbon-based network at this early stage of the synthesis. XRD data allow to highlight what kind of environment is created during the solvothermal process. It is observed here that sodium does not only react with HEP but could also be involved in many secondary reactions. Sodium compounds are good indicators of the formation of the various molecules generated by HEP degradation. For instance, the presence of Na2CO3 likely reveals the formation of H2O and CO2, according to equations 6 to 8.

Na + H2O → (Na + OH -) + ½H2 Equation 6 2OH -+ CO2 → CO3 2-+ H2O Equation 7 2Na + + CO3 2-→ Na2CO3 Equation 8 H2 once formed, can also react with sodium to form sodium hydride revealed by XRD, according to equation 9 suggested formerly by Speyer et al. [15].

Na + ½H2 → NaH Equation 9 Sodium cyanide can be formed through different pathways. Each of them involves the earlier formation of simple molecules such as NaNH2 or HCN, after equations 10 and 11 [START_REF] Rubo | Ullmann's Encyclopedia of Industrial Chemistry[END_REF][START_REF] Ray | [END_REF]:

NaNH2 + C → NaCN + H2 Equation 10HCN + NaOH → NaCN + H2O Equation 11NaNH2 can be formed by the reaction of sodium with ammonia [START_REF] Greenlee | Inorganic Syntheses[END_REF], while HCN can be produced via the reaction of methane with ammonia, respectively according to the Degussa (equation 12) and Andrussow (equation 13) processes [START_REF] Formentin | [END_REF].

Na + NH3 → NaNH2 + ½H2 Equation 12 CH4 + NH3 → HCN + 3H2 Equation 13The thermal decomposition during solvothermal reaction of sodium hydride cyanamide Na5H(CN2)3, revealed by the XRD indexation, might lead to the formation of NaCN during the solvothermal process. In the literature, the decomposition of sodium cyanamide Na2CN2 has been shown to yield NaCN, Na and N2 [START_REF] Sole | Proc. R. Soc. Lond. A[END_REF], according to equation 14.

Na2CN2 → 2NaCN + Na + N2 Equation 14 It is suggested here that a similar reaction could occur during the solvothermal process starting from Na5H(CN2)3. However, the formation of Na5H(CN2)3 during the solvothermal reaction is not well understood yet.

C.3.2. Raman spectroscopy investigations

Raman spectroscopy has been used here to better characterize the solvothermal product in terms of density of π-conjugated bonds. For the sake of comparison, spectra of HEP and sample SP-T/350 have been recorded in the range 600 cm -1 -3100 cm -1 using a 532 nm wavelength laser (cf. figure 6 and7). Two interesting features in the form of two large bands appear on the solvothermal product spectrum at 1374 cm -1 and 1605 cm -1 , corresponding to the well-known carbon bands D and G respectively [START_REF] Sasaki | [END_REF][47][48]. A high content of organized sp 2 carbon atoms is thus likely to be present after the solvothermal reaction as indicated by the G band, the intense and broad D band indicating the presence of numerous defects such as amorphous carbon, aliphatic fraction or defects in the graphenic network. Strong fluorescent background indicates the existence of high molecular weight aromatic molecules, with a high number of π-conjugated bonds. The more conjugated the π-system, the higher the Raman shift towards high wavelengths. For the sake of clarity, the fluorescent background has been removed using computational using computational tools, giving the spectrum referred as solvothermal product BR (standing for background removed). The HEP spectrum in figure 7 shows many characteristic bands in the studied range, none of them being present in the solvothermal product BR spectrum. From this observation, the conversion of HEP during the solvothermal process appears total.

D Proposed overall mechanism

The first part of this work shows the broad source of information obtained on the solvothermal reaction between sodium and 1-(2-hydroxyethyl)piperidine using complementary analytical techniques. It has been first observed that after each solvothermal reaction, an overpressure and a loss of mass were noticed. It has been shown that a direct correlation between these two phenomena cannot explain the observations. Therefore, occurrence of several reactions has been suggested from the results of mass spectroscopy and XRD investigations, mainly linked to thermal degradation of HEP and EP-ONa into smaller molecules such as CO2, H2O or CH4, as well as reactions between the excess of sodium and various species. Raman spectroscopy has revealed the early formation of sp 2 hybridized carbon atoms, in both graphitic domains and large aromatic molecules. The second part of this work aims to make bridges between these results to propose an overall mechanism of the solvothermal reaction. Based on the in-situ mass spectrometry observations, the expected reaction between HEP and Na (equation 1) is likely to start as soon as the reagents are put in contact. The stabilization of the channels related to H2 and EP-ONa after a few hours indicated nearly terminated reaction. Moreover, the total conversion of HEP to EP-ONa has been confirmed by Raman spectroscopy as no traces of the initial reagent have been detected in the final product. The large increase in H2 concentration (related to equation 1) within the first hours goes along with a significant elevation of the CO2, CH4 and NH3 levels. The formation of these products can be assigned to thermal degradation of EP-ONa, methanation reactions or reactions within the gaseous phase, e.g. ammonia formation from N2 and H2.

D.1.2. Second reactional stage

Throughout the whole solvothermal reaction, sodium can react with numerous molecules, e.g. H2, H2O, CH4, NH3, leading to the formation of numerous solid molecules such as NaH, Na2CO3, NaCN, NaOH or NaNH2. These secondary reactions are of first importance as they allow the conversion of the liquid phase, supposedly formed from a mixture of H2O, NH3 and various CxHy molecules, and a part of the gaseous phase into a solid phase. As a matter of fact, it has been observed that for Na/HEP ratio equal to unity, the conversion of liquid products is far from complete, as shown by the very oily aspect of the product recovered. On the contrary for higher Na/HEP ratio, the fraction of liquid converted to a solid phase increases. In the meanwhile, the solid phase, mainly composed of carbon, nitrogen and oxygen atoms, undergoes rearrangement reactions promoting the formation of large aromatic molecules, with a high number of π-conjugated bonds, as well as sp 2 carbon networks. Elimination reactions are likely to occur as they become predominant over substitution reactions at high temperature and under alkaline conditions [49]. Condensation reactions, aromatization reactions and cyclization reactions are likely to take place during the solvothermal process. Please do not adjust margins Please do not adjust margins At the end of the solvothermal reaction, before the opening of the reactor, the solid phase composed of a carbon-based network is likely to entrap solid sodium compounds as well as liquid or gaseous small molecules (cf. figure 8). The gaseous phase contains most of the initial N2 injected before heating up the reactor, as well as H2, NH3, CO2, CH4 and probably other light CxHy molecules. The overall gases formed during the reaction are responsible for the increase in pressure, as reported in table 1. When the reactor is opened, the sudden change in pressure, from around 75-85 bar to atmospheric pressure, is likely to cause the vaporization or sublimation of one part of the entrapped light molecules, which are ejected from the reactor. The presence of these molecules is not taken into account neither when measuring the pressure nor when measuring the mass amount. This phenomenon could explain the calculated yields below 70 wt.%.

D.2. Description of the solvothermal products

Conclusions

The solvothermal reaction between metallic sodium and 1-(2hydroxyethyl)piperidine yields a carbon-rich product, whose pyrolysis at 850 °C has been previously showed to lead to appealing 3-dimensional N-doped graphenic materials. In the present study, the impact of the solvothermal temperature and pressure, as well as the molar ratio between Na and HEP have been investigated. An increase in temperature, from 300 °C to 380 °C, has been shown to promote the conversion of the liquid phase into solid and gaseous compounds. It can be assumed that at high temperature, elimination reaction as well as other rearrangement reactions are more predominant than at lower temperatures. A rise in the Na/HEP ratio, from 1.0 to 2.0, has been revealed to contribute to various reactions between the excess of sodium and small molecules, e.g. CH4, NH3 or H2O. Change in pressure, from 100 bar to 200 bar, has not been associated with meaningful variations in the solvothermal product. However, as it has been mentioned in a previous work [33], the pyrolysis of solvothermal products from 100-bar reactions lead only to poorly organized graphenic materials. This could mean that the solvothermal pressure has an impact on the structure and microstructure of the solvothermal product, but this could not be investigated since the synthetized products are unstable under electronic beam. Complementary characterization techniques have been carried out to understand how the solvothermal reaction proceeds. Insitu mass spectrometry has shown the formation of H2, CH4, NH3, and CO2 throughout the synthesis, linked with the reaction between HEP and Na, methanation reactions as well as possible thermal degradation. The presence of these small molecules has been successfully related to XRD results, revealing the formation of various sodium compounds such as Na2CO3, NaH or NaCN. Raman spectroscopy has demonstrated the early formation of large sp 2 carbon network as well as molecules containing numerous π-conjugated bounds. The various sources of data collected have led to an overall mechanism proposal which could hopefully help to better understand this type of reactions under strongly alkaline conditions, in a temperature and pressure domain still little known in organic synthesis. The schematic development still remains incomplete, in particular because of the complexity and the diversity of the phases present. Further work would require more thorough analysis of the oils and tars formed, as done for middle temperature range conversion of biomass for instance. Besides, the link between the solvothermal mixture and the highly structured products formed after pyrolysis is currently being investigated and will be thoroughly treated in a further paper. The first results tend to show that the pressure inside the reactor plays a role on the porous structure of the graphenic foams, low pressures being more favourable to the apparition of higher porosity, at the cost of a lower crystallinity. In the other hand, the reaction temperature seems linked with the chemical composition of the graphenic foams, as well as the crystallinity and the porosity. A break in trends tends to occur between the properties of the graphenic foams elaborated from the solvothermal reaction temperature of 350 °C and 380 °C. This observation could be a marker of the transition from a subcritical reaction to a supercritical reaction. Further investigations will be carried out to clarify this hypothesis.

Figure 1 :

 1 Figure 1: Theoretical gas overpressure (grey line), experimental gas overpressure (green dashed line) and yield (dark blue dotted line) with the evolution of the solvothermal reaction temperature (a), Na/HEP ratio (b) and solvothermal reaction pressure (c)and/or hinder the reactions leading to the formation of gaseous

Figure 2 :

 2 Figure 2: Photographs of solvothermal products presenting a very oily aspect (a), oily aspect (b) and dry aspect (c)

Figure 3 :

 3 Figure 3 : Time variation of channel signals associated with H2 (a), CH4 (b), NH3 (c), H2O (d), CO2 (e) and solvothermal alkoxide EP-ONa(f); the orange dotted lines are the signal for the blank test; the black dashed line has been added on the EP-ONa graph to better visualize the apparition of this product; solvothermal conditions: 380 °C, 200 bar, Na/HEP 2.0

Figure 4 :

 4 Figure 4: X-ray diffractogram of solvothermal product SP-T/350, indexed with sodium carbonate (PDF 37-0451), sodium hydride (PDF 76-0171), sodium cyanide (PDF 37-1490) and sodium hydride cyanamide (PDF 89-1367) system than that in the solvothermal reactor, these calculations demonstrated the likely formation of ammonia during the process. C.3. Solid phase characterization C.3.1. X-Ray diffraction investigations

Figure 5 :

 5 Figure 5: Raman spectrum of the solvothermal product

Figure 6 :

 6 Figure 6: Raman spectra of HEP (blue) and solvothermal product BR (black)

D. 1 .

 1 Reaction pathway D.1.1. First reactional stage

Figure 7 :

 7 Figure 7: Drawing of the proposed structure for the solvothermal solid product before opening the reactor

Table 1 :

 1 Sample identifications, reaction parameters and experimental results

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 7

Acknowledgements

This work was supported partly by the French PIA project "Lorraine Université d'Excellence", reference ANR-15-IDEX-04-LUE.

Notes and references

Conflicts of interest

There are no conflicts to declare.


F Bessler

Lapicque

A V Shirsath

S Raël

C Bonnet

L Schiffer

W Bessler

F Lapicque

Electrochemical pressure impedance spectroscopy for investigation of mass transfer in polymer electrolyte membrane fuel cells

Keywords: Impedance spectroscopy, Pressure oscillations, Fuel cells, Mass transfer phenomena, Liquid water

Electrochemical pressure impedance spectroscopy for investigation of mass transfer in polymer electrolyte membrane fuel cells

Introduction to transfer phenomena investigations in fuel cells

Polymer electrolyte membrane fuel cells (PEMFC) in addition to being an attractive system for energy conversion, represent a huge research area for electrochemists and engineers, combining electrode mechanisms and kinetics to transport or transfer phenomena in the various porous media and flow channels of the cell. Progress in the understanding of the physical phenomena involved, actually allowed progressive, however huge improvement in fuel cell performance and efficiency. The porous materials in PEMFCs comprise the catalyst layer formed by a nanostructured carbon material supporting the catalyst clusters, and the gas diffusion layer in contact to the bipolar plates, which allow gas distribution in the fuel cell. In most current technologies, gas diffusion layers are carbon-based and consist of a macroporous substrate with 10-30 µm pores with a microporous layer directly in contact with the catalyst layer.

Modelling of the fuel cell core was initiated in the 1980s through the description of porous electrodes, where the catalyst is in the form of a thin film on the surface of the pores, or as agglomerates in a parallel-piped porous structure [START_REF] Watanabe | Controlled micropore distribution of porous gasdiffusion electrodes for a molten carbonate fuel cell by a slurry sintering method[END_REF][START_REF] Springer | Electrical impedance of a pore-wall for the floodedagglomerate model of porous gas-diffusion electrodes[END_REF]. Numerous works published since then could refine the pioneering approaches for more accurate representation of the catalyst layer [START_REF] Springer | Modeling and experimental diagnostics in polymer electrolyte fuel cells[END_REF][START_REF] Springer | Characterization of polymer electrolyte fuel cells using AC impedance spectroscopy[END_REF][START_REF] Bultel | Modified thin film and agglomerate models for active layers of P.E. fuel cells[END_REF], before more complete theoretical works could integrate physical phenomena occurring in the non-reactive gas diffusion layer.

Electrochemical impedance spectroscopy (EIS) for mass transfer phenomena in PEMFC

EIS featuring the ratio of complex voltage over current fluctuations over a large range of frequencies can show the occurrence of ohmic phenomena for frequency on the order of 1 kHz, charge transfer at lower frequencies, then mass transfer phenomena for frequencies at the Hz level or below. As for other electrochemical systems, EIS spectra of fuel cells in the Nyquist plot often consist of two main loops -with possible overlapping -, the first one in the high frequency range expressing finite charge transfer rates while the second one is usually believed to be related to control by mass transfer phenomena [4], [START_REF] Brug | The analysis of electrode impedances complicated by the presence of a constant phase element[END_REF][START_REF] Agarwal | Critical issues associated with interpretation of impedance spectra[END_REF][START_REF] Barsoukov | Impedance spectroscopy: theory, experiments and applications[END_REF][START_REF] Orazem | Electrochemical impedance spectroscopy[END_REF]. The most simple model for mass transfer derives from Fick's second law, whose complex impedance solution involves a resistance, the value of which expressing the significance of mass transfer control, and a time constant indicating on the frequency range for this loop. This time constant, which is given as the ratio of the square of a characteristic length and the effective diffusivity of the electrochemical species, indicates on the frequency range of the mass transfer phenomenon.

The above ideal description of Nyquist EIS spectra does not hold in all cases, because of the occurrence of more complex processes as discussed here.

Mass transfer processes in fuel cells cover various physicochemical phenomena in the different cell components, involving gaseous reactants or the water produced, as exemplified below.

 In the catalyst layer, in addition to the refined modelling approaches of the catalyst layer, one can mention use of EIS for evaluation of improved catalyst layer deposition [START_REF] Conde | Mass-transport properties of electrosprayed Pt/C catalysts layers for polymer-electrolyte fuel cells[END_REF], or fundamental investigations on the importance of the nature of the pore filling liquid [START_REF] Ji | Electrochemical impedance spectroscopy evidence of dimethyl-solicon-oil enhancing O2 transport in a porous electrode[END_REF] or the specific roles of catalyst clusters and carbon support [START_REF] Meyer | Electrochemical impedance spectroscopy of catalyst and carbon degradations in proton exchange membrane fuel cells[END_REF]. Moreover, in addition to its transport in the catalyst layer pores, the reacting gas dissolves in the thin ionomeric layer attached to the catalyst cluster for possible charge transfer, accompanied by H + transfer: the set of process is often viewed by EIS as an overall single phenomenon, with signature in the high-frequency loop. The lumped nature of the charge and mass transfer in the catalyst layer often results in a single loop in the Nyquist plot, making its interpretation uneasy or little accurate. Aside mass transport, adsorption of intermediate species such as hydrogen peroxide or platinum oxides on the electrode surface [START_REF] Orazem | Electrochemical impedance spectroscopy[END_REF][START_REF] Roy | Application of complementary analytical tools to support interpretation of polymer electrolyte membrane fuel cell impedance data[END_REF], is known to result in an inductive small loop below 1Hz  Whereas gas transport in the channels of the distribution (bipolar plates) is often a fast process, it is often believed to be rate-controlling in the gas diffusion layer structure. Linked to the gas permeability and the pore size, gas transport rates in the MPL -often reported as helping in the water management in the cell [START_REF] El-Kharouf | Gas diffusion media and their degradation[END_REF][START_REF] Nishiyama | Water transport characteristics in the gas diffusion media of proton exchange membrane fuel cell -Role of the microporous layer[END_REF][START_REF] Lapicque | A critical review on gas diffusion micro-and macroporous layer degradations for improved membrane fuel cell durability[END_REF] are to differ strongly from those in the macroporous substrate. Modeling the overall transfer rates for a microporous layer-coated substrate by a single time constant approach can reveal as inefficient, because of very different pore size distributions and permeability exhibited by the two layers [START_REF] Malevich | Investigation of chargetransfer and mass-transport resistances in PEMFCs with microporous layer using electrochemical impedance spectroscopy[END_REF].

 Whereas drying conditions have to be avoided, presence of liquid water (flooding) can be problematic as hindering gas transport, and even hazardous for fuel cell operation and durability. EIS has been proven an efficient diagnosis for occurrence of both undesired situations. In essence, drying is observed by increased impedance in the whole frequency range, whereas only the low frequency loop in the Nyquist plot is enlarged by occurrence of flooding [START_REF] Merida | Characterisation of proton exchange membrane fuel cell (PEMFC) failures via electrochemical impedance spectroscopy[END_REF][START_REF] Fouquet | Model based PEM fuel cell state-of-health monitoring via ac impedance measurements[END_REF][START_REF] Antonacci | Feasibility in conbining electrochemcial ilmpedance spectroscopy and synchrotron X-ray radiography for determining the influence of liquid water on polymer electrolyte membrane fuel cell performance[END_REF]. Besides, liquid flow in the cell is per se discontinuous, as slug flow in the channels and being an "eruptive" transport in the carbon layers [START_REF] Haines | Water transport across biological membranes[END_REF][START_REF] Manke | Investigation of water evolution and transport in fuel cells with high resolution synchrotron R-Ray radiography[END_REF]. Because of the little deterministic nature of liquid transport, Roy and Orazem [START_REF] Roy | Analysis of flooding as a stochastic process in polymer electrolyte membrane (PEM) fuel cells by impedance techniques[END_REF] developed a stochastic approach in EIS, flooding being revealed by far larger EIS data scattering at a fixed frequency.

EIS is usually not conducted alone to characterise the operation or the state of health of fuel cells: other electrochemical techniques e.g. voltage measurement at fixed current density (cd), cyclic or sweep voltammetry are routinely employed in complement to EIS, as well as -more recently -estimation of the limiting cd [START_REF] Conde | Mass-transport properties of electrosprayed Pt/C catalysts layers for polymer-electrolyte fuel cells[END_REF][START_REF] Chowdhury | Transport resistances in fuel-cells catalyst layers, Polymer electrolyte fuel cells 17 (PEFC 17) Book Series[END_REF]. Mechanical, analytical, spectroscopic and microscopic techniques can help in establishing reliable diagnosis on a fuel cell [START_REF] Schneider | Fast locally resolved electrochemical impedance spectroscopy in polymer electrolyte fuel cells[END_REF][START_REF] Rozeai Niya | Study of proton exchange membrane fuel cells using electrochemical impedance spectroscopy technique -A review[END_REF]. In particular, the significant issue of liquid water in the cell structure can benefit from spectroscopic techniques e.g. NMR [START_REF] Yuan | Impedance technique in PEM fuel cell diagnosis-A review[END_REF], neutron imaging and synchrotron X-ray. Synchrotron X-ray can provide local saturation in liquid water with image resolution below 10 µm in gas diffusion layer and electrode structures [START_REF] Antonacci | Feasibility in conbining electrochemcial ilmpedance spectroscopy and synchrotron X-ray radiography for determining the influence of liquid water on polymer electrolyte membrane fuel cell performance[END_REF], [START_REF] Cho | Analysis of the transient response and durability characteristics of a proton exchange membrane fuel cell with different micro-porous layer penetration thicknesses[END_REF][START_REF] Chevalier | Inoperando measurements of liquid water saturation distributions and effective diffusivities of polymer electrolyte membrane fuel cell gas diffusion layers[END_REF][START_REF] Liu | Microporous Layer Degradation in Polymer Electrolyte Membrane Fuel Cells[END_REF], in particular revealing far higher saturations over the ribs than over the channel [START_REF] Ge | Resolving the gas diffusion layer substrate land and channel region contribution to the oxygen transport resistance of a partially-saturated substrate[END_REF].

EIS although powerful may appear insufficient to explore mass transfer phenomena because of the often limited information drawn from spectra in the low frequency domain. Besides, in some cases, an inductive loop below a few Hz often corresponds to adsorption processes, which can interfere with the signature of mass transfer phenomena.

Impedance spectroscopy in electrochemistry is not restricted to correlation of current and potential: other physical variables have successfully used such as pulsation of velocity in electrohydrodynamic impedance [START_REF] Deslouis | Oxygen as a tracer for measureements of steady and trubulent flows[END_REF], or the sine weight fluctuations monitored by a quartz balance upon potential modulation [START_REF] Bourkane | Investigation of gold oxidation in a sulphuric medium: electro-gravimetric transfer function technique[END_REF]. Moreover, it has recently been suggested to extend dynamic diagnostics in cells with gaseous reactants or products to the pressure at the cell outlet. This led to the concept of electrochemical pressure impedance (or pressure-induced electrochemical) spectroscopy (EPIS) introduced by Hartmann et al. [START_REF] Hartmann | Pressure dynamics in metal-oxygen (metal-air) batteries: a case study on sodium superoxide cells[END_REF]. Development of this promising technique for investigations of low frequency phenomena occurring in a fuel cell is introduced below.

Pressure-induced impedance spectroscopy for fuel cell behaviour in the low frequency range

Gas pressure is known to be an asset to performance of fuel cells: increase in the pressure has a positive effect on the zero-current cell voltage, exchange current density and mass transfer rates [35,[START_REF] Zhang | Understanding the effects of backpressure on PEM fuel cells reactions and performances[END_REF]. The dynamics relationship between pressure and the cell voltage was shown in [START_REF] Niroumand | Pressure-voltage oscillations as a diagnostic tool for PEFC cathodes[END_REF][START_REF] Grübl | Electrochemical pressure impedance spectroscopy (EPIS) as diagnostic method for electrochemical cells with gaseous reactants: a model based analysis[END_REF] for the case of sodium-air batteries, with various impedances defined between two variables amongst voltage, current and pressure, the third one being kept constant. The potential of EPIS was soon confirmed by simulation [START_REF] Grübl | Electrochemical pressure impedance spectroscopy (EPIS) as diagnostic method for electrochemical cells with gaseous reactants: a model based analysis[END_REF], before EPIS -here defined as the ratio of voltage fluctuations over back pressure oscillations, ZV/p -was successfully performed on a fuel cell using a loudspeaker to generate pressure oscillations below 100 Hz [START_REF] Engebretsen | Electrochemical pressure impedance spectroscopy applied to the study of ^polymer electrolyte fuel cells[END_REF], although the experimental spectra were presumably subject to considerable noise.

EPIS: the techniques used

Because of the above limitations of EIS for frequencies below 1 Hz for distinct analysis of mass transfer phenomena, EPIS has been investigated through experiments and simulations in view to evaluating its potential as accurate diagnosis tool with distinct analysis of rate-controlling mass transfer processes below 1 Hz.

Experimental EPIS spectra were recorded in a 100 cm 2 single FC (UBzM) installed in a fully equipped test bench, for operation at controlled back pressure, with sustained fluctuations of the cathode outlet and at constant current [START_REF] Shirsath | Pressure-induced electrochemical spectroscopy in polymer electrolyte fuel cells: role of back pressure in the analysis of mass transport phenomena[END_REF]. Air was humidified at various levels in thermostated double-walled vessel. Stoichiometry factors for H2 and air O2 at 1.2 and 2.5 respectively, and the cell was operated at 55°C with a steady pressure at 1.15 bar. The whole system was controlled a D-space card with Matlab-Simulink environment. Fluctuating response of the potential was recorded. Back pressure amplitude was carefully selected depending on the frequency in the range 1 mHz -1 Hz to allow linear behaviour of the system through low distorted voltage signal, and appreciable voltage amplitude [START_REF] Shirsath | Pressure-induced electrochemical spectroscopy in polymer electrolyte fuel cells: role of back pressure in the analysis of mass transport phenomena[END_REF]. In complement, an EPIS model based on the in-house code DENIS [START_REF] Grübl | Electrochemical pressure impedance spectroscopy (EPIS) as diagnostic method for electrochemical cells with gaseous reactants: a model based analysis[END_REF][START_REF] Bessler | A new framework for physically based modeling of solid oxide fuel cells[END_REF] with a (1D+1D) model has been developed, with Navier-Stokes equation for gas transport in the channels and diffusion-convection equation in the gas diffusion layers [START_REF] Bessler | A new framework for physically based modeling of solid oxide fuel cells[END_REF][START_REF] Schiffer | Model-based analysis of electrochemical pressure impedance spectroscopy (EPIS) for PEM Fuel Cells[END_REF]. The electrode is considered as flat at the membrane interface, and water transport in the membrane covers electroosmotic drag and back-diffusion.

Preliminary results and conclusions

Comparison of EIS and EPIS is shown on an example in Figure 1 in both Bode and Nyquist plot. The modulus of ZV/p exhibits a sigmoidal variation with the frequency, being at a fraction of µV/Pa at a few mHz and attaining 5 µV/Pa at 1 Hz. In agreement with data in [START_REF] Niroumand | Pressure-voltage oscillations as a diagnostic tool for PEFC cathodes[END_REF], the phase shift angle decreased regularly from nearly 0° at 1 mHz to approx. 190° at a Hz. On the contrary, corresponding EIS spectra ZV/I exhibit conventional profiles (Fig. 1c andd), in particular showing very little changes below 1Hz.

Water management investigations in the test bench showed that for low air humidification (RH=20%), no liquid water is formed in the cell, whereas for fully hydrated air, the production rate of liquid formed is nearly proportional to the cd. At low cd, EPIS modulus increases with the excitation frequency, with little effect of the hydration level (Fig. 2a). On the contrary, at high cd (1 A/cm 2 ), the presence of liquid water is visible in EPIS spectra below 50 mHz, whereas the impedance is far less affected by the relative humidity for higher frequencies (Fig. 2b). In all cases, the phase shift angle followed similar variations in all hydration and cd conditions.

The pressure dynamics within the cell can be influenced by other system components. We have therefore investigated the dynamic ratio of fuel cell inlet over fuel cell outlet pressure, 𝑍 ( )/ ( ) . In the model we have included a 0-D humidifier upstream of the cell. Boundary conditions were constant mass flow into the humidifier and oscillating back pressure at fuel cell outlet. Figure 3 shows a comparison of experimental and simulated 𝑍 / spectra. The data show a drop in modulus at around 0.1 Hz, associated with a drop in phase angle from 0° to -90°. By varying the humidifier volume in the model, we could show that this response is dominated by the humidifier. We conclude that the pressure response of balance-of-plant components need to be included into the EPIS interpretation.

These preliminary data show the potential of EPIS for transport/transfer investigation in fuel cells in the low frequency domain, which can be a complementary technique with a technology based on simple mechanical devices (loudspeakers or pressure valves), far cheaper than spectroscopic equipments.

Figure 1 .

 1 Figure 1. Comparison of EPIS measurements at 0.4 A/cm 2 , air humidified at 55%; T(fuel cell) = 55°C; (a) Bode plot (b) Nyquist plot, with EIS measurements (c) Bode plot (d) Nyquist plot. The solid lines in Nyquist plots b) and d) are only for illustrative purpose, not related to modelling

Figure 2 . 2 Figure 3 .

 223 Figure 2. Modulus of the voltage-over-pressure impedance at two cd levels, depending on the relative humidity (RH) of the air injected; (a), top, 0.2 A/cm 2 ; (b), bottom 1 A/cm 2

Acknowledgments

The work was partially funded by the German Research Federation (DFG) (no. 391377080) and French National Research Agency in the frame of Epistel project (ANR-2017)