Federica Agostini

Ivano Tavernelli

Giovanni Ciccotti
email: giovanni.ciccotti@roma1.infn.it

Nuclear Quantum Effects in Electronic (Non)Adiabatic Dynamics

Keywords: XX, XX, XX No code given

Trajectory-based approaches to excited-state, nonadiabatic dynamics are promising simulation techniques to describe the response of complex molecular systems upon photo-excitation. They provide an approximate description of the coupled quantum dynamics of electrons and nuclei trying to access systems of growing complexity. The central question in the design of those approximations is a proper accounting of the coupling electron-nuclei and of the quantum features of the problem. In this paper, we approach the problem in the framework of the exact factorization of the electron-nuclear wavefunction, re-deriving and improving the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm recently developed to solve the exact-factorization equations. In particular, a procedure to include quantum nuclear effects in CT-MQC is derived, and tested on a model system in different regimes.

Introduction

Simulating the excited-state dynamics of complex molecular systems is a challenging task that concerns a growing community of theoretical chemists, chemical physicists and material scientists. The dynamical response of systems like graphene quantum dots [1], the provitamin D 3 molecule [2], or an organic chromophore [3], after photoexcitation can only be described in simulation in an approximate way. The reason is indeed the exponential growth of the computational cost with the number of degrees of freedom for solving the quantum-mechanical problem exactly. In order to realistically predict molecular structures that are products of photo-activated reactions, these approximations have to correctly capture the coupling between electronic and nuclear motion, that is source of phenomena such as nonadiabatic transitions and quantum decoherence, or also account for the quantum nature of the nuclear degrees of freedom.

Methods for excited-state molecular dynamics range from quantum wavepacket simulation schemes [4][5][6][7][8][9] and semiclassical approaches [10][11][12][13][14][15][16][17], to trajectory-based techniques with travelling Gaussian wavepackets [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33], with quantum trajectories [34][35][36][37][38][39][40][41][42], or as quantum-classical approaches [43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62]. Among these, the most promising to gain access to molecular systems of increasing size are Send offprint requests to: clearly trajectory-based techniques, based on a simplified description of nuclear dynamics in terms of trajectories, and exploiting the possibility of on-the-fly calculation of electronic-structure properties. In this paper, we will focus on the trajectory-based scheme derived from the exact factorization of the electron-nuclear wavefunction [63,64]. As extensively discussed in the literature [65,66,52,[67][68][69], the exact factorization represents a suitable starting point for developing quantum-classical approximations to the quantum-mechanical electron-nuclear problem. In fact, the (molecular) time-dependent Schrödinger equation is decomposed into coupled electronic and nuclear evolution equations, that naturally lend themselves to the developments of approximations that (i) can be designed to affect only the nuclear subsystem, (ii) can be used to perform the quantum-to-classical limit- [START_REF] Vleck | Proc. Nat. Ac. Sci[END_REF][START_REF] Wyatt | Quantum dynamics with trajectories. Introduction to quantum hydrodynamics[END_REF] on the nuclear dynamics governed by a timedependent Schrödinger equation. The work proposed here, still within the exact factorization, is an attempt to improve the approximations involved in (ii).

It has been shown [53,54] that the nuclear time-dependent Schrödinger equation can be mapped onto the coupled Hamilton-Jacobi and continuity equations describing, respectively, the evolution of the phase of the nuclear wavefunction and the nuclear density. Hamilton-Jacobi equation could be solved with the method of characteristics, generating trajectories produced by a force term accounting for the effect of the electrons (either adiabatic or nonadiabatic), via time-dependent scalar and vector potentials [65,66,[START_REF] Curchod | [END_REF], and a contribution to the force that depends on the, so-called, quantum potential [34-36, 38, 37, 39,40, 71, 73, 74]. The quantum potential is a well-known quantity in the field of Bohmian dynamics, as it allows to reproduce quantum effects within a trajectory-based description of the dynamics. It depends, however, on the nuclear density, or more precisely on the second spatial derivative of the modulus of the wavefunction, thus it couples the Hamilton-Jacobi equation for the phase to the continuity equation for the density. An approximate quantum-classical solution of the coupled electronic and nuclear equations of the exact factorization, from which the coupled-trajectory mixed quantum-classical (CT-M-QC) algorithm [53] has been derived, discards the quantumpotential term in the force. This approximation simplifies the numerical procedure, as it decouples the Hamilton-Jacobi equation from the continuity equation.

We reconsider this approximation, and suggest a new solution in two steps. First, we propose to compute the quantum-potential contribution to the force within a modified version of the CT-MQC algorithm, in which the nuclear density is reconstructed as the sum of travelling Gaussians centered at the positions of the trajectories. The merit of this improvement is that it still does not require the explicit solution of the continuity equation, and thus it does not affect heavily the computational cost of a CT-MQC calculation. Second, we apply the method of characteristics to the Hamilton-Jacobi equation, finding Hamilton's equations generating the characteristics. Although the idea behind this procedure has been already sketched in previous work [54], we provide here the, still lacking, detailed derivation.

A first example, in which the effect of the quantum potential is observable, illustrates an adiabatic (in the electronic ground state) tunneling process. It yields results not only in qualitative, but also in quantitative, agreement with exact calculations. Instead, in our second example, a strong nonadiabatic process in which the role of the quantum potential is less important, the new approximation is able to improve the trajectory-based calculations.

The paper is organized as follows. In Section 2 we recall the exact factorization, and discuss within our approximation the procedure to solve the time-dependent nuclear Schrödinger equation, to determine the phase and the density of the nuclear wavefunction. The illustrative one-dimensional two-state model, solved in the two cases discussed above, is presented in Section 3 together with simulation details. Section 4 discusses the numerical results and validates the proposed approximation. Conclusions are contained in Section 5.

Exact factorization of the electron-nuclear wavefunction

The exact factorization of the electron-nuclear wavefunction [63,64] provides a prescription for decomposing the time-dependent Schrödinger equation (TDSE) for a system of interacting electrons and nuclei into the coupled dynamics of the subsystems, i.e., the electronic and the nuclear.

The time-dependent molecular wavefunction, Ψ (r, R, t), is the solution of the TDSE ĤΨ = i∂ t Ψ , with Hamiltonian Ĥ(r, R) = Tn (R) + ĤBO (r, R), containing the nuclear kinetic energy, Tn , and the electronic Born-Oppenheimer (BO) Hamiltonian, ĤBO , defined as the sum of the electronic kinetic energy and of the interaction potentials. Here, the symbols r, R indicate all electronic and nuclear coordinates, respectively. The full wavefunction can be exactly written as the product

Ψ (r, R, t) = χ(R, t)Φ R (r, t), (1)
where χ(R, t) can be considered a genuine nuclear wavefunction, yielding the exact nuclear many-body density and current density, and Φ R (r, t), the electronic function, depends parametrically on the nuclear configuration. Note that the partial normalization condition dr|Φ R (r, t)| 2 = 1 ∀ R, t has to be imposed to guarantee that |χ(R, t)| 2 reproduces at all times the nuclear density.

Inserting the factorization (1) into the TDSE, and using the partial normalization condition [75,76], one finds the coupled evolution equations for the two components of the molecular wavefunction, namely

ĤBO + Û coup en -Φ R (r, t) = i ∂ t Φ R (r, t) (2) Nn ν=1 [-i ∇ ν + A ν ] 2 2M ν + χ(R, t) = i ∂ t χ(R, t), (3)
where the new quantities introduced will be discussed below. The derivation of these equations can be found in Ref. [64] and is detailed for full understanding in Appendix A. Nuclear masses are indicated by the symbol M ν , with the index ν running over the N n nuclei. In the electronic equation ( 2), the operator Û coup en [Φ R , χ] [67] couples the electronic evolution to the nuclear dynamics as it depends on the nuclear wavefunction,

Û coup en [Φ R ,χ] = Nn ν=1 1 M ν [-i ∇ ν -A ν ] 2 2 + -i ∇ ν χ χ + A ν • (-i ∇ ν -A ν ) . ( 4
)
Note that in the second term the gradient in (-i ∇ ν χ)/χ only acts on the nuclear wavefunction. The scalar potential, or time-dependent potential energy surface (TD-PES) [65,77,66,78,[START_REF] Curchod | [END_REF][79][80][81][82], (R, t), and the time-dependent vector potential [START_REF] Curchod | [END_REF]69,[83][84][85][86], A ν (R, t), are defined by

(R, t) = Φ R (t)| ĤBO + Û coup en -i ∂ t |Φ R (t) r , (5)
and

A ν (R, t) = Φ R (t)| -i ∇ ν Φ R (t) r , (6)
respectively, where • r stands for an integration over the electronic coordinates. In the nuclear TDSE (3), the timedependent potentials fully account for electronic nonadiabatic effects, i.e., excited-state effects, on nuclear motion. It has been shown [68] that in the limit of infinite nuclear masses, or in the limit m/M → 0 [87], Eqs. ( 2) and ( 3) reduce, as it would be good to expect, to the standard BO (adiabatic) equations. The (R, t)-dependent phase transformations

χ(R, t) → χ(R, t) = e -i θ(R,t) χ(R, t) (7)
Φ R (r, t) → ΦR (r, t) = e i θ(R,t) Φ R (r, t), (8)
clearly leave unchanged the product form of the molecular wavefunction, Eq. ( 1). Under such transformation, the time-dependent potentials transform as

˜ (R, t) = (R, t) + ∂ t θ(R, t) (9) Ãν (R, t) = A ν (R, t) + ∇ ν θ(R, t). ( 10
)
These transformation relations are easy to prove, after some algebra, if the expressions of χ(R, t) and ΦR (r, t) are replaced into the definitions of the time-dependent potentials given in Eqs. ( 5) and ( 6). The evolution equations are form-invariant under a gauge transformation. Uniqueness of the solution of Eqs. ( 2) and ( 3) is guaranteed up to a chosen gauge.

Aiming at developing a strategy to efficiently solve the coupled electronic and nuclear equations of motion, we write the nuclear wavefunction in polar form, i.e., χ(R, t) = exp[(i/ )S(R, t)] |χ(R, t)|. By inserting this expression into the TDSE, applying the differential operators, and separating its real and imaginary parts, we derive an evolution equation for the phase

∂ t S = - ν [∇ ν S + A ν ] 2 2M ν -- ν -2 2M ν ∇ 2 ν |χ| |χ| , (11)
= - ν [∇ ν S + A ν ] 2 2M ν -( + Q) , (12)
which implicitly defines the quantum potential Q(R, t), and another for the density

Γ (R, t) = |χ(R, t)| 2 ∂ t Γ = - ν ∇ ν • Γ ∇ ν S + A ν M ν . (13)
Eq. ( 11) can be viewed as a Hamilton-Jacobi equation, if we take the right-hand side as a Hamiltonian, whose kinetic term is ν [∇ ν S + A ν ] 2 /(2M ν ), and the potential term is

V(R, t) = (R, t) + ν -2 ∇ 2 ν |χ(R, t)| 2M ν |χ(R, t)| . ( 14
)
The quantity Q(R, t) in Eq. ( 12) depending on the modulus of χ is the, so-called, quantum potential. Eq. ( 13) is a continuity equation for the nuclear probability density, whose probability current density is

J ν (R, t) = Γ (R, t) ∇ ν S(R, t) + A ν (R, t) M ν . ( 15
)

Method of characteristics to solve the Hamilton-Jacobi equation

The solution of the Hamilton-Jacobi equation ( 11) can be constructed by employing the method of characteristics. Rearranging Eq. ( 11), we introduce the function [88]

F (∇S, S t , R, t) = S t + H (∇S, R, t) = 0 ( 16
)
where we have indicated with ∇S the 3N n -dimensional gradient of S(R, t) and with S t its partial time derivative. Note that the Hamiltonian explicitly depends on time via the TDPES and the quantum potential. The solution S depends on the variables R, t, while the Hamilton-Jacobi equations depends non-linearly on the first derivatives of S with respect to these variables. Before giving the characteristic equations of the Hamilton-Jacobi equation, we denote ∇S = P, such that H(∇S, R, t) = H(P, R, t), and we introduce the following partial derivatives

∂ Rν F (P, S t , R, t) = ∂ Rν H (P, R, t) (17) ∂ t F (P, S t , R, t) = ∂ t H (P, R, t) (18) ∂ Pν F (P, S t , R, t) = ∂ Pν H (P, R, t) (19) ∂ St F (P, S t , R, t) = 1 (20) ∂ S F (P, S t , R, t) = 0, (21)
with ν = 1, . . . , 3N n , as in Appendix B. The characteristic equations yield the "evolution" of R(s), t(s), P(s), S t (s), S(s) as functions of the parameter s. We have (see Eqs. ( 68), ( 70) and ( 72))

Ṙν (s) =∂ Pν H (P(s), R(s), t(s)) (22) ṫ(s) =1, ( 23
) Ṗν (s) = -∂ Rν H (P(s), R(s), t(s)) (24) Ṡt (s) = -∂ t H (P(s), R(s), t(s)) (25)
Ṡ(s) = 3Nn ν=1 P ν (s)∂ Pν H (P(s), R(s), t(s)) -H(P(s), R(s), t(s)). (26)
In our case the characteristic equations to be explicitly solved are

Ṙν = P ν (t) + A ν (R, t) M ν (27) Ṗν = -∇ ν H n (P(t), R(t), t) (28)
Ṡ = Nn ν=1 P ν (t) • P ν (t) + A ν (R, t) M ν -H n (P(t), R(t), t), (29)
while Eq. ( 23) is trivial and Eq. ( 25) is just a consistency relationship.

The symbol H n has been used to indicate the Hamiltonian as written in the right-hand side of Eq. ( 11),

H n (P(t), R(t), t) = Nn ν=1 [P ν (t) + A ν (R, t)] 2 2M ν + (R, t) + Q(R, t), (30)
thus the right-hand side of Eq. ( 28), considering that

A ν • (∇ ν A ν ) = ∇ ν (A ν • A ν )/2, is Ṗν = -∇ ν Nn ν =1 P ν (t) • A ν (R, t) M ν + A 2 ν (R, t) 2M ν + (R, t) -∇ ν Q(R, t). (31)
In sight of the gauge freedom ( 7) and ( 8), we can choose a phase θ(R, t), satisfying the equation

Nn ν=1 P ν • ∇ ν θ M ν - (∇ ν θ) 2 2M ν + ∂ t θ = - Nn ν=1 P ν • A ν M ν + A 2 ν 2M ν + ( 32
)
such that the term in parenthesis in Eq. ( 31) can be set to zero. Notice that the quantum potential Q(R, t) is not a gauge potential, since it depends on the modulus of the nuclear wavefunction, which is not a gauge-dependent quantity. Eq. ( 31) differs from the expression of the force that has been derived in previous work [53,54,[START_REF] Min | [END_REF]55], since there we have worked in the Lagrangian frame and we have, a bit abusively, called conjugate momentum the linear momentum M ν Ṙν .

In principle, the solution of the Hamilton-Jacobi equation, S(R, t), is obtained by integrating the ODEs for an infinite number of initial conditions, to determine the value of S at time t and position R, thus obtaining S(R, t).

Computing the nuclear density and the quantum potential

Solving the Hamilton-Jacobi equation (11) employing the method of characteristics described in Section 2.1, together with the continuity equation (13), is fully equivalent to solving the nuclear TDSE (3). However, our aim here is mainly to propose and justify an approximation which decouples the equation for the phase (11) from that for the denisty (13).

The approximation we suggest tentatively for the probability density associated to the nuclear wavefunction appearing in the expression of the quantum potential is

|χ (R, t)| 2 = 1 N tr Ntr I=1 Nn ν=1 G σ I,ν R ν ; R (I) ν (t) (33)
where

G σ I,ν (R ν ; R (I)
ν (t)) are normalized Gaussians centered at the classical trajectory R (I) (t) and the variance

σ I,ν , if n (I)
tr is the number of trajectories whose R (J) (t) falls in a sphere of given radius (typically 1 Å in Ref. [START_REF] Min | [END_REF]) centered at R (I) (t), is given by [54]

σ I,ν = D 2 I,ν -D 2 I,ν n (I) tr . (34)
Here,

D I,ν = 1 n (I) tr n (I) tr J=1 R (I) ν (t) -R (J) ν (t) (35)
D 2 I,ν = 1 n (I) tr n (I) tr J=1 R (I) ν (t) -R (J) ν (t) 2 . ( 36
)
Each Gaussian in Eq. ( 33) is associated to a constant weight, 1/N tr , that guarantees the normalization of the density. To avoid misunderstanding, let us stress that this choice does not mean equivalent weights associated to the trajectories, since the initial conditions are sampled from the given nuclear probability distribution at time t = 0. This means that we expect to have more trajectories in regions of higher probability at time 0. Of course there is no guarantee that this initial distribution will be conserved at long times, but it can be considered reasonable at least at relatively short times.

It is worth noting that this Ansatz, Eq. ( 33), was already used in Refs. [53,54,[START_REF] Min | [END_REF] to solve the electronic part of the problem, but then it was not used to compute the quantum potential in Eq. (31). The novelty here is that we want to explore possible improvements obtainable by using consistently the same Ansatz also to compute the quantum potential.

Once this improved approximation described above has been defined, we can use the coupled-trajectory mixed quantum classical (CT-MQC) algorithm discussed in Refs.- [53,54,[START_REF] Min | [END_REF] to follow the dynamics of our system. In the next section we describe details of the simulations used to test this approximation.

Simulation details

We select a one-dimensional model to illustrate the performance of the new CT-MQC algorithm. In particular, we (i) investigate how our new quantum potential can provide strong nuclear quantum effects such as the spatial splitting of a nuclear wavepacket induced by tunneling or by electronic nonadiabatic transitions, and (ii) check how well the distribution of trajectories reproduce the real nuclear density of our model.

The chosen system is described by the Hamiltonian in the diabatic basis

Ĥ = P 2 2M + V + (R) d d V -(R) , (37)
where the 2 × 2 matrix is the representation of the "electronic Hamiltonian", Ĥel , in the standard 2-state "electronic" basis (| ↑ , | ↓ ). There, the diagonal elements of the electronic Hamiltonian are

V ± (R) = aR 2 ± bR + c, (38)
and the parameters are a = 1.0 a.u., b = 3.5 a.u., c = 3.0625 a.u., and d = 1.0 a.u. in the first case (set 1), while in the second case (set 2) we change d = 0.25 a.u. . The nuclear mass is chosen to be M = 2000 a.u., as the proton mass in atomic units. The adiabatic PESs (BO PESs) and the adiabatic states are given by

E ± (R) = V + (R) + V -(R) 2 ± 1 2 (V + (R) -V -(R)) 2 + 4d 2 (39)
and

|±; R = N (R) |↑ + E ± (R) -V + (R) d |↓ , (40)
with

N (R) = d/ d 2 + (E ± (R) -V + (R))
2 a normalization term. We report in Fig. 1 the two BO PESs and the nonadiabatic couplings ϕ

(+) R |∂ R ϕ (-) R
for both sets of parameters. As it is clear from the figure, changing the value of the off-diagonal elements in Eq. ( 37) allows to tune the coupling between the two electronic states, and thus to change the nonadiabatic character of the problem. In both cases, the initial condition is chosen such that

|Ψ (R, t = 0) = 4 1 πσ 2 e -(R-R 0 ) 2 2σ 2 e i P0(R-R0) |-; R , (41)
namely a Gaussian-shaped nuclear wavepacket centered at R 0 = -1.75 a.u. with variance σ = 0.2 a.u. is placed "on" the ground-state (lower) surface and, so to say, directed toward the nonadiabatic coupling region, with average linear momentum at time t = 0 is P 0 = 90 a.u. or P 0 = 150 a.u., for the first and for the second study, respectively. The ini- tial nuclear density is also shown in Fig. 1. The two sets of parameters have been chosen to observe different physical processes, namely tunneling through the electronic ground-state barrier (set 1) and nonadiabatic populating the excited state from the ground state (set 2).

Exact propagation is performed in the diabatic basis with time-step dt = 0.1 a.u. using the split-operatortechnique [90] with

e -i Ĥdt e -i Ĥel dt 2 e -i P 2 2M dt e -i Ĥel dt 2 . ( 42
)
The "electronic" part has been solved as shown in Appendix C. As for our approach, we generated N tr = 1000 trajectories by solving the Lagrange equivalent of Eqs. ( 27) and ( 28) using the velocity-Verlet algorithm. The electronic evolution is computed as shown in Refs. [53,54]. The time step dt = 0.1 a.u. is used for the overall integration.

Due to the limited sampling of trajectories, some features of the nuclear density and, thus, of the quantum potential might be not properly reproduced within the CT-MQC procedure. To circumvent this problem, we apply a Gaussian smoothing to the density and the quantum potential computed according to Eqs. ( 33) and ( 12), respectively. The Gaussian smoothing consists in convoluting the function of interest with a suitable Gaussian, thus replacing the initial function with a smoothed representation of the original one. In our simulations, we have used a value of the Gaussian variance between 0.4 and 0.6 a.u., checking the stability of the results.

Illustrative examples

In this section, at first we focus on the process of tunneling in the adiabatic regime. The initial kinetic energy of the travelling wavepacket is not sufficient for it to classically overcome the potential barrier of the ground-state PES. Additionally, once again due to the low momentum of the nuclear wavepacket, and in consequence of the weak nonadiabatic coupling, no population transfer from the electronic ground state to the excited state is observed during the dynamics. Then, we have chosen a stronger coupling between the electronic states so that the initial momentum of the nuclear wavepacket is large enough for the wavepacket to split passing through the avoided crossing, while undergoing a nonadiabatic transition to the excited state.

Tunneling effects in the adiabatic regime

We report in Fig. 2 the transmission probability, i.e., the probability of the nucleus to be in the region R > 0 a.u. The exact calculation shows in Fig. 2 that for t ≥ 70 a.u. almost 30% of the density is transferred across the barrier. In the trajectory-based calculation (CT-MQC algorithm or our present approach) we try to estimate the transmission probability as the fraction of trajectories which have R > 0 a.u. at time t. We see in Fig. 2 that by CT-MQC the transmission probability is less than 2%, while in our approach it is correctly reproduced. A convergence test has been performed on the transmission probability for 100, 1000, and 10000 trajectories. The (systematic) error on the final value of the probability decreases from 16%, to 5% and 3%, respectively. We directly compare the nuclear density Γ (R, t) with our and the CT-MQC approximation given in Eq. (33). Once again, in Fig. 3 we show exact results (black lines), the results based on CT-MQC (green lines) and our results (red lines, undistinguishable from exact). Without the quantum potential, as expected from previous observations, the approximate density is qualitatively different from reference results. The inclusion of the force term containing the effect of the quantum potential clearly improves the behavior of density.

Nonadiabatic regime

The parameters of the potentials and the initial conditions for the dynamics have been modified as described in Section 3 to increase the nonadiabatic character of the simulated process. We show in Fig. 4 the agreement between the population of the ground and excited states computed exactly and according to the two approximations under test. As expected, the ground-state population decreases in time, while the population of the excited state increases with time, since during the passage of the nuclear wavepacket through the region of coupling, a nonadiabatic transfer from the ground to the excited state takes place. To analyze the effect of the quantum potential in this case, we compare in Fig. 5 the exact nuclear probability density with our approximation and with CT-MQC. It is evident that the effect of the quantum potential is weaker than in the previous model study, since CT-MQC is already capable of reproducing the spatial splitting of the nuclear density. However, our approach allows to slightly improve on CT-MQC approximation the trajectory distribution. Specially in the last snapshot shown in Fig. 5 it is clear that our approach produces (i) a good shifting of the peak of the distribution, and (ii) a more convincing zero in the region between the two peaks. Both features improve the agreement with exact results.

Conclusions

In this work we have revisited the CT-MQC approach to quantum nonadiabatic dynamics improving that approximation by including an estimate of the quantum potential which in CT-MQC is put to zero. In doing that we have also intended to clarify the derivation of the approach with special emphasis on the solution of the nuclear dynamics by characteristics. We hope that this clear-cut reconstruction of the meaning of the solution by characteristics together with the choice of an initial nuclear density physically meaningful (a wavepacket) will give a stil intuitive but reasonable justification for the way in which we reconstruct in this approach the nuclear density, i.e., also the quantum potential, at time t. By a suitable illustration using a spin-boson model we have been able to discuss tunneling effects and nonadiabatic dynamics. From these two effects we have shown that the new approach is an objective improvement on CT-MQC.

we have to consider only the infinitesimal variations of the variables x, u, p which leave the value of F unchanged, that is

dF = F p • dp + F u du + F x • dx = 0. ( 64
)
Geometrically this condition means that the vector (dx, du, dp) is orthogonal to the gradient (F x , F u , F p ) of F . The equations of any curve (characteristics) which satisfy this condition, i.e., each curve which stays in the tangent bundle of the manifold generated by F ,

ẋ = dx ds u = du ds ṗ = dp ds (65)
define variables x, u, p which lie on the isosurface F (p, u, x) = 0. From that we get the solution of the PDE for every point x reached by the curve. Normally the number of curves needed to solve completely the PDE is defined by the boundary condition, i.e., by all points x 0 ∈ ∂Ω, u 0 , p 0 , that lie in the boundary of the domain in which we need to solve the PDE. The next step is to derive the explicit form for the characteristic equations (65). To this end, let us look at the function (Eq. ( 60)) F (∇u(x), u(x), x) = F (x) = 0 for any x ∈ Ω as a function of the variable x. That implies

∂F ∂x i = n j=1 F pj (p(s), u(s), x(s))u xj xi (x(s)) + F u (p(s), u(s), x(s))p i (s) + F xi (p(s), u(s), x(s)) = 0, (66)
with i = 1, . . . , n. In Eq. ( 66) we used the definition

p i (x) = ∂u(x) ∂x i = u xi (x) . (67)
If we assume that ẋj (s) = F pj (p(s), u(s), x(s)),

and consider that ṗi

(s) = d ds ∂u(x(s)) ∂x i = n j=1 ∂ ∂x j ∂u(x(s)) ∂x i dx j (s) ds = n j=1 u xj xi (x(s)) ẋj (s), (69)
then we see that the first term on the right-hand side of Eq. ( 66) can be identified with Eq. ( 69), i.e.,

Replacing this expression in Eq. ( 66) we obtain ṗi (s) = -F u (p(s), u(s), x(s))p i (s) -F xi (p(s), u(s), x(s)). where we have used Eqs. ( 67) and (68).

We shall now prove by replacing Eqs. ( 68), [START_REF] Vleck | Proc. Nat. Ac. Sci[END_REF], and (72) in Eq. ( 64) that infinitesimal variations of x, u, p generated moving along the characteristics, we leave F unchanged. In fact, ∀ s

C Decomposition of the electronic Hamiltonian

In this appendix, having separated by Trotter the evolution operator into nuclear (kinetic energy) and "electronic" parts, we write down the explicit expression of the "electronic" evolution operator.

The model Hamiltonian introduced in Section 3 is written as a linear combination of the Pauli matrices and the identity matrix, that is

Ĥel (R) = V + (R) d d V -(R) = α(R)σ x + β(R)σ z + γ(R) Î. ( 74
)
The elements of Ĥel (R) are all real, thus the contribution of σy is identically zero. It is easy to prove that α(R) = d and that, using γ(R) + β(R) = V + (R) and γ(R) -β(R) = V -(R), the explicit values of β(R) and γ(R) are easily computed, but not needed. Using this decomposition, the expression

e -i dt 2 Ĥel = e -i dt 2 [α(R)σx+β(R)σz+γ(R) Î] , (75)
where we have used the property σ2n i = Î and σ2n+1

which gives the action of the Ĥel evolution operator.

Fig. 1 .

 1 Fig. 1. Adiabatic (BO) PESs computed by diagonalizing the 2 × 2 electronic Hamiltonian in Eq. (37), shown as continuous dark-green (parameter set 1) and dashed green (parameter set 2) lines. The nonadiabatic coupling are shown as continuous blue (parameter set 1) and dashed cyan (parameter set 2) lines. The initial nuclear density is also plotted as a red line.

Fig. 2 .

 2 Fig. 2. Transmission probability across the barrier of the ground-state PES for the parameter set 1 of Fig. 1. The black curve shows the quantum-mechanical result, computed from the exactly propagated nuclear density. Our results are shown in red (1000 trajectories), dark-pink (100 trajectories) and pink (10000 trajectories), while the green line shows the CT-MQC result.

Fig. 3 .

 3 Fig.3. Snapshots at 50, 100, 140 a.u. (from left to right) of the exact nuclear density (black lines), and of the nuclear density constructed according to Eq.(33). Red lines, our results (superimposed to the exact); green lines, CT-MQC.

Fig. 4 .

 4 Fig. 4. Population of the electronic states in the adiabatic basis as functions of time. The reference result (black lines) is computed by solving exactly the quantum dynamics as described in Section 3. Our results are shown as red lies, while CT-MQC results are shown as green lines.

Fig. 5 .

 5 Fig. 5. Snapshots at 20, 40, 60 a.u. (from left to right) of the exact nuclear density, of the one from the present work and from CT-MQC. The color code is the same as in Fig. 3.

u

  xj xi (x(s))F pj (p(s), u(s), x(s)) = ṗi (s).

( 71 )

 71 Finally, to determine the last equation for the characteristics, let us differentiate u(s) with respect to s. s)F pj (p(s), u(s), x(s)),[START_REF] Curchod | [END_REF]

F

  pj (p(s), u(s), x(s)) -F u (p(s), u(s), x(s))p j (s)-F xj (p(s), u(s), x(s)) + F u (p(s), u(s), x(s))p j (s)F pj (p(s), u(s), x(s)) + F xj (p(s), u(s), x(s))F pj (p(s), u(s), x(s)) = 0. (73)

  introduced in Eq.(42), can be written as, applying again the Trotter decomposition, ase -i dt 2 [α(R)σx+β(R)σz+γ(R) Î] = e -i dt 2 [α(R)σx+β(R)σz] e -i dt 2 γ(R) Î e -i dt 2 α(R)σx e -i dt 2 β(R)σz e -i dt 2 γ(R) Î .(76)Then, expanding each exponential in a Taylor series, we get e -i dt 2 α(R)σx = cos α

  n a positive integer and i = x, y, z, ande -i dt 2 β(R)σz = cos β2 [γ(R)+β(R)] cos α(R)dt 2 -ie -i dt 2 [γ(R)-β(R)] sin α(R)dt 2 -ie -i dt 2 [γ(R)+β(R)] sin α(R)dt 2 e -i dt 2 [γ(R)-β(R)] cos α(R)dt 2   ,

Acknowledgements

We are delighted to dedicate this small piece of work to our friend Hardy for his 65th birthday. Fond as we are of his highly collaborative and hospitable character we wanted to give a worth contribution to these qualities. Therefore, we made a special effort to contribute to his recent ideas by clarifying once and for all the concept of solution by characteristics. We really hope Hardy will enjoy as much as we have.

Authors contributions

All authors were involved in the preparation of the manuscript. All authors have read and approved the final manuscript.

A Coupled electron-nuclear equations

In this appendix we derive the coupled electronic and nuclear equations, Eqs. ( 2) and (3), starting with the exactfactorization Ansatz of Eq. ( 1) and the general expression of the TDSE for the molecular wavefunction Ψ (r, R, t).

The factored form of Ψ is inserted into the TDSE, that is

When multiplying by Φ * R (r, t) and integrating over electronic coordinates, we get for the term on the left-hand side

and for the term on the right-hand side

In Eqs. ( 44) and ( 45), • r indicates an integration over electronic coordinates and the partial normalization condition Φ R (t)|Φ R (t) r = 1 ∀R, t has been used. The righthand side of Eq. ( 44) can be rewritten using for the second term

and identifying in the third term the time-dependent vector potential,

In sight of Eq. ( 43), and equating ( 45) and ( 47), the timederivative of the nuclear wavefunction is

Adding and subtracting the term A 2 ν /(2M ν ), and defining the TDPES as

Eq. ( 48) becomes

which is a standard TDSE. Now we move to the derivation of the electronic evolution equation (2). Before proceeding, we rearrange Eq. ( 50), wherever χ = 0, as

In Eq. ( 43), we consider the right-hand side

and the left-hand side

Equating the right-hand sides of Eqs. ( 53) and ( 54), we obtain

where in going from the first to the second equality, the first term on the right-hand side has been replaced by Eq. ( 52). To write Eq. ( 56) in the form given in the main text, Eq. ( 2), we need few further steps. We start selecting in the first and second lines of Eq. ( 56) the three terms indicated below, and see that they can be written as

then we gather the terms that in Eq. ( 56) are multiplied by (-i ∇ ν χ)/χ, and we end summing them together to find the quantity

ûν appears in the definition of Û coup en [Φ R , χ], the electronnuclear coupling operator, in Eq. ( 4). The evolution equation for Φ R finally reads

that is Eq. ( 2).

B Characteristics of a partial differential equation

The method of characteristics to solve a partial differential equation (PDE) aims at translating the PDE into infinitely many first-order ordinary differential equations (ODEs) that compute curves (the characteristics) providing a pointwise solution of the PDE. The procedure followed in Section 2.1 to determine the characteristics of the Hamilton-Jacobi equation associated to the nuclear TDSE is an application of the method discussed below [88].

Let us consider the general form of a non-linear PDE

where the unknown function u depends on the coordinates x ∈ Ω ⊆ R n , and we have indicated the derivatives of u with respect to the coordinates as p = ∇u. Let us introduce a family of curves described parametrically x = x(s) = [x 1 (s), . . . , x n (s)] (61) u = f (s) (62) p = g(s) (63) with i = 1, . . . , n and s ∈ I ⊆ R. These curves are the characteristics of the PDE (60) if, for each value of the parameter s, u = u(x) is a solution of the PDE with p(x) = ∇u(x). Therefore, we need to identify "evolution equations" for x(s), u(s), p(s) that guarantee the condition F = 0, i.e., Eq. ( 60), for any value of s. To this end,


Guillaume Guimbretière

Benjamin Pillot

Christophe Révillion

Corrie Mathiak

Romain Authier

Territorial solar baking potential: A socio-technical analysis for planning the sustainable energy transition of the bakery world

Keywords: Artisan bakery, solar energy, spatio-temporal modeling, social-technical transition, multi-scale organization, small island

In the context of a global energy transition, heating remains a highly energy consuming process. One very virtuous and low economic cost solution for cooking is direct solar concentration oven. In this article, in order to estimate the solar baking potential of a territory, a socio-technical solar bakery's workflow model and two socio-economic indicators were developed. From real incoming solar energy as input, the model takes into account the solar oven properties and baker's time management in order to calculate bakery's bread productivity. Then, this model was applied and indicators computed for a small tropical island (Réunion, France). Results show that substantial solar bread production areas match the island population density, and that the level of solar bread production can be viable, though requiring an adaptation of the bakery sector's organization to the solar production variability.

Introduction

Recently, the sharp rise of energy and electricity prices in France has jeopardized many small-scale commercial activities such as artisan bakeries [START_REF] Monde | Crise énergétique : des centaines de boulangers manifestent à paris pour demander davantage d'aides au gouvernement[END_REF]. It reveals that the energy transition framework is still in its infancy, and how the pursuit of a sustainable path must incorporate more than just the energy supply dimension [START_REF] Köhler | An agenda for sustainability transitions research: State of the art and future directions[END_REF][START_REF] Kern | Restructuring energy systems for sustainability? energy transition policy in the netherlands[END_REF]. Political actions still look like small patches on big issues that fail to address the big picture. A typical example are variable renewable energies (VRE) being integrated into centralized electrical networks (e.g. solar & wind farms). In the short and mid-term, this is likely a good option for mitigating CO 2 emissions from the energy sector [START_REF] Clarke | Energy Systems[END_REF] and thus for combating climate change. However, in the longterm, the sustainability of such VRE systems should be questioned : They 1) are much more adapted to small-scale remote applications [START_REF] Pillot | Historical trends in global energy policy and renewable power system issues in Sub-Saharan Africa: The case of solar PV[END_REF], and 2) strongly depend on mining and rare earth elements (REE) [START_REF] Dutta | Global demand for rare earth resources and strategies for green mining[END_REF].

Due to its abundance and worldwide distribution, solar energy has long been identified as the main energy source from which a sustainable society could be built [START_REF] Georgescu-Roegen | Myth about Energy and Matter[END_REF]. The current issue is that nowadays energy sources in general and solar energy in particular are mostly regarded as power sources, i.e. electricity generation sources [START_REF] Kannan | Solar energy for future world: -A review[END_REF]. Overall, the energy-centered view of the energy transition framework is too constraining as it only focuses on the energy part of the technical system. It is key to understand that energy is not an end in and of itself; rather, it is a means towards fulfilling a need. According to the energy services definition by Fell et al. [START_REF] Fell | Energy services: A conceptual review[END_REF], the desired end services or states are obtained or facilitated by functions performed using energy. These functions can be in part unmaterialized but in the end, it is only a hardware tool (here called technology) that can materially act on the real world and fulfill the need (i.e. obtain the desired service). Therefore, by realizing that energy has to be combined with a technology in order to obtain a final end service, it becomes clear that only by addressing both technologies and energy can a sustainable transition of our civilization be achieved. It is important to note that "technologies" adress the whole technological regime, not only the energy "production" technologies. Furthermore, when considering the place and role of the technologies in the energy transition, one has to consider that technologies are not neutral tools isolated from the social framework [START_REF] Ellul | Le système technicien[END_REF][START_REF] Illich | La convivialité[END_REF]. A radical energy transition towards sustainability should only be possible once the socialtechnical system has been considered as a whole [START_REF] Kemp | Technology and the transition to environmental sustainability: the problem of technological regime shifts[END_REF][START_REF] Geels | Technological transitions as evolutionary reconfiguration processes: a multi-level perspective and a casestudy[END_REF]. Accordingly, a broader understanding of the systemic synergies and complexities implied by energy transition incites to apprehend the whole chain of energy supply, [energy] transport and [energy] consumption. It also implies social and cultural changes [START_REF] Coignard | Solar panels-is more always better?: Assessing the carbon footprint of communities[END_REF]: we must transform our manner of living; that is to think about why we need specific types of energy, and how we may find corresponding sustainable alternatives.

What does this shifted perspective of the energy transition imply on the baker's socio-technical transition towards sustainability? In fact, more than 50% of the energy consumed by a bakery is dedicated to baking [START_REF] De | metiers et de l'artisanat Region Rhones-Alpes, Comprendre ses consommations d'energie pour mieux les gererboulangers, patissiers[END_REF]. If our approach is articulated around energy production technologies, it appears that the bakery oven primarily requires heat energy; nowadays this heat is mostly converted directly from combustion of wood [START_REF] Manhiça | Wood consumption and analysis of the bread baking process in wood-fired bakery ovens[END_REF] or fossil fuels, or from the electricity provided by the national grid using an electric oven [START_REF] Therdthai | Recent advances in the studies of bread baking process and their impacts on the bread baking technology[END_REF]. Energy transition is traditionally regarded as an energy-centered transition [START_REF] Slameršak | Energy requirements and carbon emissions for a low-carbon energy transition[END_REF][START_REF] Leach | The energy transition[END_REF]. The energy system being multi-scalar, it is possible to make energy supply more sustainable according to different scales [START_REF] Clarke | Energy Systems[END_REF]. As a result, the current approach consists in replacing utility-scale power plants based on finite resources by renewable energy (RE) systems. Another way is to substitute the electricity supply locally, using a stand-alone RE system such as solar PV panels coupled with storage [START_REF] Meyer | Bake bread while the sun shines: Solar bakery for off-grid rural community development[END_REF]. In both cases, the energy transformation chain is based on two conversions, from primary renewable energy (light, motion, etc.) to electricity and then from electricity to heat. Meanwhile, the related energy technology consumes limited materials, such as REE [START_REF] Dutta | Global demand for rare earth resources and strategies for green mining[END_REF], subject to electronic failures, and complex to the extent that most of the technical maintenance cannot be performed by end users [START_REF] Pillot | Historical trends in global energy policy and renewable power system issues in Sub-Saharan Africa: The case of solar PV[END_REF]. One can immediately see the limits of such an energy-centered view of the energy transition where the baking technology (the oven) is left aside (this is still an electric oven). However, bread is eventually baked in the oven, not in the power plant. A sustainable baking tech should therefore 1) be based on a one-conversion chain, i.e. from primary to heat energy, 2) use the least possible material and REE, 3) be very reliable in the long-term, and 4) be easily usable and maintainable by the bakers. The direct solar-thermal oven does fit this definition: even with the lack of a complete life cycle analysis, with no electrical part, no emissions of CO 2 , GHG or small particles during operation, one can legitimately expect a high degree of sustainability for this kind of technology.

Several solar baking ovens have already been set up throughout the world [START_REF] Müller | Solar community bakeries on the argentinean altiplano[END_REF][START_REF] Ayub | Solar Thermal Application for Decentralized Food Baking Using Scheffler Reflector Technology[END_REF][START_REF] Ayub | Energy-and exergy-based thermal analyses of a solar bakery unit[END_REF][START_REF] Forsythe | Experimental performance assessment of a solar powered baking oven[END_REF][START_REF] Mishra | Solar thermal powered bakery oven[END_REF]; in Normandie (France), a artisan bakery has been using a solar oven for commercial purposes for the past few years [START_REF] Neoloco | Artisanat solaire -Graines NeoLoco[END_REF]. Solar cooking technology has been known for a very long time [START_REF] Kundapur | 250 years saga of solar box cookers[END_REF]. However, in spite of its multiple benefits, its use is not widespread among household users in developing countries due to various reasons identified by P.P. Otte [START_REF] Otte | Solar cookers in developing countries-what is their key to success?[END_REF]. Typically, criteria that influence the adoption of solar cooking are economic, social, cultural, environmental, political and technical. Unfortunately, in the design of cookstoves, engineers mostly tend to pass over usability and the needs of the targeted users in their context of use, while focusing on technical attributes only [START_REF] Moses | Development of a practical evaluation for cookstove usability[END_REF]. Those constraints will probably be quite different for craftsmen against household users. Since they are professionals, their motives will differ: promoters of the direct solar ovens may be more concerned by global warming while bakers may be more concerned by economic viability of the solar technology (even if they are heavily concerned by environmental impact of the baker activity). Besides, as professionals, if the business activity is viable they should also be able to better adapt their practices and time management than household users. In the case of Neoloco [START_REF] Neoloco | Artisanat solaire -Graines NeoLoco[END_REF], the bakery's products and the organization of its activity have been adapted to the use of the direct solar-thermal oven.

With the sustainable transition of part of the bakery sector in mind, the direct solar baking practice needs to be upscaled to the regional scale. To do so, various obstacles have to be overcome. First, only few studies concern collective solar ovens instead of household's ones. Second, a craftsman needs to plan his activity in advance. For this, because of the variability of the solar resource, it is essential to anticipate the potential of solar bread production in the context of the bakery's location site. This solar production depends on the type of solar oven and the solar irradiation level of the site, but also on the organization of the baker's work. In this article, in order to provide some answers to the planning issues of the baker, we propose a sociotechnical model for assessing the production potential of a solar baker. Then, at the territorial scale of the bakery sector, we discuss two socio-economic indicators giving information on the overlap between the consumer market and the solar producers' sector. The insights gained from the socio-technical model and the socio-economic indicators are discussed in the context of Réunion, a French department located in the south-west of the Indian Ocean, where solar energy is abundant and, according to the 2015 Energy Transition Act, France's energy policy is to reach full energy self-sufficiency by 2030 [START_REF]Clean energy for eu islands -la réunion[END_REF].

Material and methods

Réunion island

The Réunion is an island located in the Indian Ocean, about 800 km east of Madagascar (Fig. 1). With an area of 2,512 km 2 , this territory has an ovoid shape with a diameter between 50 and 70 km. It culminates at 3 069 m on the Piton des Neiges. The terrain of Réunion is characteristic of volcanic islands with steep slopes and alluvial plains of limited size [START_REF] Cazes-Duvat | Les littoraux des Mascareignes entre nature et aménagement[END_REF].

With nearly 900,000 inhabitants, the island has a population density close to 350 inhabitants/km 2 . The majority of these inhabitants are located along the coast (Fig. 1) [START_REF]Base de population[END_REF]. The climate is tropical, with a hot and humid season from December to March, and a cooler and drier season from June to September. In addition, the insularity, the circulation of trade winds and the high relief of Reunion Island generate significant climatic contrasts: The windward coast (east and north) can be distinguished from the leeward coast (west and south). The windward coast receives the trade winds from the front, bringing masses of humid air which are blocked on the relief of the center of the island by pouring a great part of their rains there. In fact, the leeward coast has a much drier climate. This geographical contrast also affects the sunshine rate, the coasts and the west of the island being sunnier than the east and the center. It is therefore relevant to study in more detail the spatial distribution of solar radiation over this small area. Solar thermal bakery is mainly fueled by direct solar radiation, which is concentrated through parabolic or flat mirrors into the solar cooker's focus [START_REF] Ayub | Solar Thermal Application for Decentralized Food Baking Using Scheffler Reflector Technology[END_REF][START_REF] Wentzel | The development impact of solar cookers: a review of solar cooking impact research in South Africa[END_REF]. In first instance, knowledge of both 1) spatial and 2) temporal distribution of the solar direct irradiance is essential to assess how solar baking may match the current and future behavior of potential bread consumers. Satellite-based models stand for a good trade-off between accuracy and the need for apprehending the geographical distribution of solar radiation [START_REF] Salazar | Solar irradiance time series derived from high-quality measurements, satellite-based models, and reanalyses at a near-equatorial site in Brazil[END_REF]. However, all models are not equal in terms of temporal and spatial resolution as well as final precision or data access [START_REF] Salazar | Solar irradiance time series derived from high-quality measurements, satellite-based models, and reanalyses at a near-equatorial site in Brazil[END_REF][START_REF] Sengupta | The National Solar Radiation Data Base (NSRDB)[END_REF]. The satellite-based physical model (PSM-3) developed for more than 20 years by the NREL to produce the National Solar Radiation Database (NSRDB) provides good balance between all these requirements for Réunion. Data access is also kept as straightforward as possible in order to get a full dataset of hourly Direct Normal Irradiance (DNI) maps in W/m 2 , and accuracy is close to the one from other satellite-based products such as CAMS and SARAH2 [START_REF] Salazar | Solar irradiance time series derived from high-quality measurements, satellite-based models, and reanalyses at a near-equatorial site in Brazil[END_REF][START_REF] Kenny | Which gridded irradiance data is best for modelling photovoltaic power production in Germany?[END_REF], but with a better spatial resolution than those two over Réunion because of the satellite's location : NSRDB is based on the Meteosat IODC geostationary satellite against Meteosat 11 for CAMS and SARAH2 over Réunion [37,[START_REF] Gracia-Amillo | Adapting PVGIS to trends in climate, technology and user needs[END_REF][START_REF]NSRDB | International Data[END_REF].

Topography shadow correction of the DNI

Another significant information for solar thermal bakery is the topography of the territory, that directly influences available direct solar radiation through terrain shadowing. In order to integrate the shading effect into the DNI maps, one can use the SRTM data released by the CGIAR-CSI, the horizon model instantiated in [START_REF] Pillot | Development and validation of a new efficient SRTM DEMbased horizon model combined with optimization and error prediction methods[END_REF] and the solar position algorithm developed at NREL [START_REF] Reda | Solar position algorithm for solar radiation applications[END_REF]. Accordingly, original NSRDB-based hourly DNI maps have been downscaled to the resolution of the 90m Digital Elevation Model (DEM) of Réunion. Finally, the result is a set of hourly shaded DNI maps over a 3 years range from 2017 to 2019; to give an idea to the reader, the daily mean direct normal irradiation (in kWh/m 2 ) for the 2017-2019 period is depicted in Fig. 2.

A socio-technical solar bakery workflow model

The main objective of our study was to find out how a network of artisan bakers equipped with concentrating solar ovens could fulfill population needs for baked bread using the solar resource. Hence a socio-technical workflow model was built. It is depicted in Fig. 3 with bottom-up flow. On the left of the workflow is a downward arrow representing the evolution of need from the original need for baked bread as food: the need for baking bread leads to a need for energy which must be mobilized by a cooking technology.

The bakery's process workflow model is composed of a series of numeric filters and functions pictured in full red or green rectangles on the Fig. 3. The output data of each step is the input data of the next step and are pictured in rounded rectangles, lightly shaded red, green or blue. For instance, starting at the model's beginning: The first input in the numeric model is the solar energy time series data (yellow rounded rectangle at the bottom). After numeric filtering by the DNI threshold filter, only time steps with DNI value above the DNI threshold (> 700 W/m 2 in the present case study) are kept in the Output Rated Energy time series data. This later data set is next numerically filtered though the Number of continuous hours filter, etc . . . Now, regarding the meaning of the model in details : For solar concentration, aggregated and averaged DNI gives spatial information about what we might call the gross solar income of a given territory. But, for activities such as bread baking, the core of a successful implementation is also linked to heat supply stability. For this to happen, DNI must remain above a given threshold for a time period corresponding to the duration required for preheating the oven plus baking one batch of bread. This threshold will differ from oven to oven, depending on their specific energy focus capacity (based on mirror reflection rate and total surface, oven insulation, etc.). Therefore, the DNI is two times filtered because of the response of the solar oven to solar radiation: a first high-pass power numeric filter removes values of solar power that are under a critical threshold for heating the solar oven (the DNI threshold filter). A second numeric filter is related to the minimum continuous time period above the threshold of the first numeric filter so that the oven maintains the necessary temperature to bake a full batch of bread (the number of continuous hours numeric filter). At the output of these two filters, one obtains the number of periods in the year during which the solar oven can, if one needs it, bake a batch of bread: the number of working periods. Then, depending on the baker's time management, a number of periods used for baking can be set. The oven allows for a specific quantity of bread to be be baked per batch (in kg) and from it, one can estimate the daily productivity (in kg/day). Productivity is mainly based on the daily craftsman's time management and from it, a production (in kg/year) can be calculated, depending on the mid (weekly) to long (yearly) time organization of the baker. This organization results from the integration of social constraints like rate consumption, food practices and the environmental constraint of solar energy availability. Different types of solar ovens can be analyzed based on this process workflow model, through parameterization of the filters: DNI threshold, number of continuous hours and quantity of bread to be baked per batch. Also, the model can be applied to other territories (real or model) by changing the DNI data input. In fact, the DNI and the type of solar oven define what can be called the solar baking physical potential of the bakery's site. Effectively, the available baking time consists of time when the baking chamber remains in a given operating temperature range. The targeted temperature range is fixed by physicochemical processes of bread baking [START_REF] Therdthai | Recent advances in the studies of bread baking process and their impacts on the bread baking technology[END_REF][START_REF] Mondal | Bread baking-a review[END_REF], and the heat comes from a physical process of solar flux conversion in the oven. In the model, the physical potential of solar baking is expressed in number of available working periods. In the following, a day is called bakeable when multiple continuous hours of high DNI allow for baking at least one batch of bread. Next, in order to analyse the physical potential of a Lytefire 11 oven (see subsec. 2.4) in Réunion, two new parameters were built from the number of working periods output: the Maximum number of bakeable days per year and the Average production shortage per year (see Fig. 3, at the right). These parameters are thresholds that correspond to the baker's time periods for baking bread throughout the year, and emphasize both the frequency and intensity issues of the activity. In other words, frequency is given by the number of effective bakeable days as well as the amount of time corresponding to bread shortage (unbakeable days). Intensity corresponds to the number of bread batches that can be baked on a given bakeable day. Shortage period is retrieved by collecting the different spans of consecutive unbakeable days along the year. Finally, the solar baking physical potential managed by the baker through a work strategy results in a bread daily productivity and yearly production.

The craftsman's solar bakery oven

For the present work, the Neoloco [START_REF] Neoloco | Artisanat solaire -Graines NeoLoco[END_REF] solar baking setup is considered: it uses a Lytefire 11 [43], a solar oven model allowing up to 38 kg of bread to be baked per batch. The Lytefire 11 has a 11 m 2 surface of solar energy collection made up of 23 mirrors. Collected light is then concentrated in the oven consisting of two floors: a lower floor into which the light enters through a window, and which contains absorbent elements performing the conversion of light into heat and limiting losses by outgoing reflections. The second floor is the baking chamber to which the craftsman accesses through a door. The two floors are open to each other and the heat circulation is done by natural convection. The solar collector/concentrator and oven assembly form an integral system mounted on wheels to follow the course of the sun. All operations are carried out manually and no electricity is required to operate the oven. Following Neoloco's practice feedback, the DNI threshold high-pass numeric filter is set to DNI > 700 W/m 2 , and the number of continuous hours numeric filter to 3 hours: 2 hours of preheating and 1 hour of baking. The maximum load of the Lytefire 11 is 38kg of bread per batch (the quantity/batch filter).

The baker's work strategy

Once the solar baking physical potential is estimated through the identification of the periods allowing bread baking (the "working periods"), the baker must organize all the stages of his production so that the time of baking matches the solar oven "working period". Thus, the baker has to develop a work strategy for the short (daily) and long (yearly) term. For the daily time management, regarding our hypotheses, it is assumed that the maximum number of batches over the year is estimated by maximizing the daily number of batches up to 3 for every bakeable day. For the yearly strategy, two strategies are considered: the "max production strategy" and the "weekly operating strategy". In the "max production strategy" we consider that the baker produces the maximum bread quantity allowed by the solar baking physical potential and within the limit of the adopted daily time management (max of 3 batches per day). In the "weekly operating strategy", the baker bakes only from one to two days a week. The chosen bakeable days are days where the maximum number of batches is in the 1 to 3 range.

From an artisan bakery to the Réunion's bakery sector

The bakery's process workflow in Fig. 3 allows for the estimation of a baker's production according to the solar oven model, the management of his working time and the solar resource of his baking site. Because of the long term work objective of gaining new insight into the sustainable development transition of the territory, one wants to have an island-wide vision of the potential of a solar-powered bakery sector. Classically, the matching between production and consumption is one of the key-points for the viability of an economic activity. Regarding consumption, in Réunion, the population is heterogeneously distributed (Fig. 1). Concerning the production, one main difference between the solar bakery compared to the classical electrified one is that the energy source is the solar radiation that is also heterogeneously distributed across the island (Fig. 2). For these reasons, two indicators were built in order to analyze the solar baking production/consumption potential match. First, following a systemic view of the territory sustainability, we argue that sole identification of the most efficient areas for solar bread production is not enough: energy and CO 2 emissions savings that may result from the solar oven utilization cannot be transferred to transporting solar bread to the consumers or the consumers to the bakery. That is why a built indicator called Eaten Where Solar Baked gives insight into whether bread consumption matches solar bread production in the immediate vicinity. This indicator is useful for identifying the region where even if the solar bread potential is high, the population of bread consumers is too low and the solar bread production from these regions requires transport to other areas. Based on this indicator, locations that are unsuitable for installing bakeries can be excluded. In order to assess where bakeries should be established, a second indicator called Solar Bread Economic Intensity was built. This later indicator estimates whether production matches consumption and assesses the most relevant regions for a potential solar bakery economy.

Solar bakery sector's indicators

Here, dealing with the proximity between bakery sites for solar bread production and the consumers, a honeycomb mesh of the territory is defined and then, both indicators in every hexagonal cell are computed. The mesh resolution is based on a specific distance threshold from bakery to bread consumer. This threshold may depend on the mobility scenarios to be implemented, including sustainable ones where mobility-based CO 2 emissions would not balance CO 2 emissions spared through the use of a solar bakery. Small enough cells might for instance ensure soft mobility scenarios, where customer travels are realized on foot or by bike.

Overall, bread consumption in cell i is defined by both frequency and intensity: does an individual consume bread and how much quantity if he does? Probability p c that consumed bread quantity will be higher than the bread produced by the solar bakery is given by:

p c = P q c ⩾ X p η c (1)
where q c is the final quantity consumed, η c is the number of bread consumers and X p is the supplied bread quantity. The number of bread consumers in cell i is given by the number of inhabitants η h in that cell and the probability that they actually consume bread. One may retrieve the probability P h,c of a human h to be a bread consumer c in cell i from the distribution of the ratio between the number of bread consumers and the number of inhabitants among all cells. Using the expectation of that ratio E (η c /η h ), the number of bread consumers η i c in cell i is therefore given by:

η i c = η i h • P h,c = η i h • E (η c /η h ) (2)
where η i h is the number of individuals in cell i. Expectation of the ratio between number of consumers and number of inhabitants among geographic cells is given by:

E (η c /η h ) = N i p i η i c η i h ( 3
)
where p i is the probability of the outcome η i c /η i h and N i p i = 1. By definition, considering the expected value as the weighted average of all cell ratios, p i can be expressed as the population weight of cell i, such as :

p i = η i h N h (4)
with N h the total number of inhabitants. One therefore gets:

E (η c /η h ) = 1 N h N i η i c = N c N h (5)
with N c the density of bread consumers and N h the population density in the island. Therefore, from the probability defined in Eq. ( 1), we can define the "Eaten Where Solar Baked" index EWSB i in cell i:

EWSB i = P        q i c ⩾ X i p N h η i h N c        (6)
where X i p is the produced bread quantity in cell i, and q i c the final quantity consumed in cell i. The EWSB index stands for the probability that the local production of solar bread in cell i is actually consumed in cell i. Accordingly, from a sustainability perspective, the bread is eventually baked and consumed in a sustainable manner, i.e. solar bread is supplied to the consumers through sustainable trips. Next, the "Solar Bread Economic Intensity" index SBEI i in cell i is defined as the product between the expected rate of return of the bread production and the EWSB index:

SBEI i = X i p • α • EWS B i ( 7
)
where α is the rate of return ( €) of 1 tonne of bread. The SBEI index expressed in € per cell underlines among the best EWSB sites which one are the best for a sustainable economic bakery's activity; that is the sites where a significant solar bread production can actually be consumed in the surroundings.

Bread consumers in Reunion

To calculate both solar bread indicators described in the next section, a specific population location data is used. This set is calculated from population data on administrative entities, the population being distributed over all buildings of more than 20 m 2 located in these entities, proportionally to the size of the buildings [START_REF]Base de population[END_REF]. Accurate information is retrieved from this data and allows for calculating the number of inhabitants η h and number of bread consumers η c among the cells in which a potential solar bakery could be implemented. Hypotheses about bread consumers profile have been simplified, with yearly individual consumption of bread taken as a normal distribution N (µ, σ) with µ = 50 kg [START_REF]Des habitudes de consommation en évolution[END_REF] and σ = 20 kg.

Reference bakeries

A survey to 63 french bakeries (2008-2010) [START_REF] De | metiers et de l'artisanat Region Rhones-Alpes, Comprendre ses consommations d'energie pour mieux les gererboulangers, patissiers[END_REF] shows that on average, 53% of the energy consumed by the bakeries is for baking. The average energy consumption for baking depends of the fuel type of the oven : ≈ 26000 kWh/y for electric ovens, ≈ 63000 kWh/y for gas ovens, and ≈ 98000 kWh/y for fuel oil ovens. In 2022, we conducted a series of interviews with bakers. Elements are summarized in Table 1

Results and discussion

This work targets an estimation of the solar baking potential of a territory depending on its solar energy and population context. Since solar baking is a solar technology, first a gross estimation can be made from the cumulative data of incoming solar energy (Fig. 2): in the first subsection, such gross estimation is discussed. To go further, in the second and third subsections, the results of the "maximum production strategy" and "weekly operating strategy" described in subsection 2.5, are respectively discussed. For both strategies, first the solar baking physical potential from a purely energy and technology point of view without taking into account the current population distribution, is discussed. Then, the solar baking socio-technical potential by mapping the Eaten Where Solar Baked and Solar Bread Economic Intensity indicators, that take population density into consideration, is discussed.

Solar baking potential from DNI data

Regarding solar variability, the "Daily mean shaded DNI map" (Fig. 2) shows an heterogeneous spatial distribution of the solar energy in the range 1.9 -6.5 kWh/m 2 . As a gross estimation of the solar bakery potential, the most irradiated areas are expected to have the highest potential: the highest DNI values are located in the northeast and west coasts, as well as in the southwest. The slopes of the summits receive little solar radiation due to shadowing effect, while the Piton de la Fournaise volcano in the south-east gets the lowest DNI potential because of micro climate effects. Ultimately, the high summits of the island and the cirque of Mafate are the regions with the highest solar potential. It is important to underline the large size of the satellite pixels projected on a territory with very sharped and rugged field. Future work will engage into enhancement of these satellite-based data, as well as comparison with data from ground stations.

The maximum production strategy

Remember that in the "maximum production strategy", the baker produces the maximum bread quantity allowed by the solar baking physical potential and in the limit of the adopted daily time management : 3 batches per day maximum.

Solar baking physical potential

Once DNI was identified for the period, the related thresholds (subsection 2.3) for solar bread production were derived. "Max number of bakeable days", "production shortage" and "Max number of batches" are respectively depicted in Fig. 4a, Fig. 4b and Fig. 4c. "Maximum solar bread production" is then derived from the range of batches (1-3) per bakeable day combined with the maximum quantity baked per batch. In Lytefire 11 user guide, the batch quantity range is given between 30 and 38 kg, and the result when the batch is set to 38 kg is described in Fig. 4d. It is interesting to note that while the DNI map of Fig. 2 gives information about spatial distribution of the solar resource, the set of figures in Fig. 4 actually gives more insight into the temporal variability of the solar radiation, and thus into the bread production we may derive from it. Color scales are bicolor gradients with tipping values chosen at significant thresholds: for the "max number of bakeable days" map, the tipping value is 175 days which is close to half a year. For the "average production shortage" map, the tipping value is 5 days which is the shelf life of a sourdough bread. For the "max number of batches" map, the tipping point is 280 batches which is close to half the maximum value on the island. 280 batches is also 1 4 of the theoretical maximum number of batches (i.e. 1095 = 3 * 365) and 75% of the theoretical maximum of bakeable days per year (i.e. 365).

As explained in the introduction, it is important to note that the need for energy is not a need per se but a need directly connected to a social service through a technological tool using energy. In the case of the solar bakery, the bakery's process workflow transforms energy and raw materials into baked bread using a direct concentration solar oven. Therefore, this process transfers the issue of solar variability to the issue of temporal and spatial variability of bread production. The maps in Fig. 4 represent several parameters that allows for the consequences of this variability on the production potential of solar bread and its regularity across Réunion to be taken into account. Regarding the various parameters of the solar bread production potential (Fig. 4): essentially, the thresholds maps rather fit the DNI maps (Fig. 2); the sunniest areas are those where there is the largest "max number of bakeable days per year", the lower "average production shortage per year", the "maximum number of batches per year" and therefore ultimately the greater "bread production potential". However, the quantitative outcome is significantly different with a solar bread production range being two times the DNI range. In more details:

• The max number of bakeable days per year map (Fig. 4a) shows that in the areas exhibiting the highest number of bakeable days (in green), baking can be performed more than half of the days a year.

• The average production shortage per year map (Fig. 4b) shows that except in the Piton de la Fournaise volcano area at the south-east (in pink), in the rest of the island there is less than 5 consecutive unbakeable days in average. In a nutshell, it means that baking can be performed at least once a week. It is compatible with some sourdough bakery's business model like the Neoloco one [START_REF] Neoloco | Artisanat solaire -Graines NeoLoco[END_REF].

• The max number of batches per year map (Fig. 4c) shows that in the area where the number of bakeable is greater than 175 (more than half of the days a year), the number of batches is well above 175: two or three batches per day should be a common situation.

• The solar bread production potential map (Fig. 4d) shows information close to the max number of batches per year map: in order to obtain the solar bread quantity, the number of batches are just multiplied by the oven capacity of : 38kg per batch. Nevertheless, on one hand this map gives the bread quantity in tons that can be directly compared to a typical bread production or consumption. On the other hand, the bread quantity scale's levels are chosen in order to identify sites in the sunniest areas where the solar bread production potential is the best: the summits and the Mafate mountain cirque at the island's center are the most relevant areas. On the coast, Sainte-Marie, Sainte-Suzanne (north-east) and Saint-Pierre (south) are the best populated areas for solar baking.

Solar baking socio-technical potential

Here, the study is positioned in the context of the 2023 population density and soft mobility system (bicycles and pedestrians). Obviously, it is not the scenario of the whole mobility system shift to soft mobility because the current distribution of the population on the island is highly dependent on the dominant mode of transportation, which is the private car. It is likely that a whole transition to a soft mobility system would lead to a change in the distribution of the population on the island. Soft mobility is considered only for solar bread collection, because it would be counterproductive to reduce electrical energy consumption and polluting emissions with solar baking while transferring this energy consumption to polluting trips required for collecting the bread. Moreover, the soft mobility framework is deeply sustainable and robust, even when facing with a context of energy crises which affect not only energy required for baking bread, but also mobility.

The spatial distribution of both the EWSB and SBEI index over Réunion is depicted in Fig. 5, for hexagonal cells of 500 m radius and a fixed rate of return α = 8 €/kg. The EWSB index (Eq. 6, Fig. 5a) gives information on the ability of the population to consume the entire solar bread production in an area very close to the solar bakery site. It is an accessibility index showing how accessible a batch of bread would be to the surrounding consumers when it is actually produced. The index range is from 0 (the bread is not consumed) to 1 (all bread is consumed). High accessibility does not mean that the solar bakery would eventually be viable, but only that the bread which is effectively produced is highly accessible to the local consumers. Globally, without surprise, the EWSB indicator map shows high index values in the area where the population density is the highest: in the cities of the north (Saint-Denis), the north-west (Le Port), the east (Saint-Andre and Saint-Benoit), and the south (Saint-Louis, Saint-Pierre and Le Tampon). The high potential solar bread production areas of the center of the island (see Fig. 4d) where population density is very low exhibit very low EWSB values.

The SBEI index (Eq. 7, Fig. 5b) shows the areas where the solar bread production is the highest AND can be consumed in the immediate vicinity of the solar bakery site. It gives combined information about solar bread production and accessibility, and hence shows how viable a solar bakery might be. Viability is defined locally, that is the density of potential consumers is computed from the number of inhabitants and surface area of each pre-defined cells. Finally, it is important to note that obviously the resulting map does not give a full realistic view of how much a solar bakery might be viable locally, because of the simplified hypotheses about individual consumption, the rate of return, etc. However, it gives a very interesting relative information about the most suitable locations for developing solar baking in the island. These most suitable locations are especially in cities of Saint-Pierre, Saint-Andre, Saint-Benoît and the eastern part of Saint-Denis. Finally, Le Tampon, Le Port and the western part of Saint-Denis where the EWSB index is high are not suitable places for solar baking because of the too low solar bread production potential.

A detailed discussion of the map is beyond the scope of this paper, but one can note here the sharp split-up between west Saint-Denis and east Saint-Denis. This pattern is probably due to the large size of the DNI pixels (see Fig. 2). No higher resolution DNI map is currently available, but as said previously (see section 3.1), we will study the pixel size effect in future works.

The weekly operating strategy

Remember that the "weekly operating strategy" consists in baking at most two days a week. The chosen bakeable days are days with the maximum number of batches in the range 1-3 (see Fig. 3 at right). This solar baker organization comes from the consideration of breads' types that the artisan can bake: for instance, Neoloco bakes sourdough breads one day a week. The sourdough breads are long shelf life breads (3 to 5 days) compare to only yeast breads (1 to 2 days), and they suit the solar baking variability.

Solar baking physical potential

Fig. 6 shows the maps of the solar bread production potential for the weekly strategy. According to this strategy, the baker manages his time with at least 5 days a week dedicated to activities besides baking. Following the solar energy availability, the baker could possibly bake only one day a week or consider some weeks with no production. One can compare this map with the map corresponding to the maximum production strategy (Fig. 4d). The color scales are the same. The map of Fig. 6 shows that globally, the best area for solar baking remains the same with the weekly strategy as for max production strategy: the coast, especially the north-east, and in a lesser extent the south.

Figure 6: Yearly solar bread production for a weekly operating strategy in Réunion Fig. 7 plots the yearly solar bread production of the weekly strategy versus the production of the max production strategy (black dots). Lines are empirical references that correspond to a percentage of the max production strategy minus one ton: the blue, yellow, red and green lines respectively correspond to a percentage of 100, 70, 60 and 50 %. The trend of the production switch between the maximum and the weekly strategies is linear with a 1 tonne systematic decrease (about 26 batches). First, one can see that all the locations are between the 100 % (blue dashed) and the 50 % (green dashed) lines. Next, one immediately distinguishes a first group of locations where the production is very high according to the max production strategy (between 16 and 23 t/y), and falls down between 8 and 11 t/y based on the weekly strategy. The production of these locations is very impacted by the implementation of the weekly strategy: points are close to the 60 % line (dotted red), while the scatter plot of all other locations is centered on the 70 % line (dotted yellow). However, the production of the locations following the weekly strategy are close to the maximum theoretical production of 11,856 kg (38 * 3 * 2 * 52). The maximum theoretical production with the maximum production strategy is 41,610 kg (38 * 3 * 365) and almost every site under study is below half this value. Finally, for a weekly strategy where more than 70 % of the time (5 days out of 7) is dedicated to other activity than baking, in the majority of places considered for solar baking, the average decline of the maximum production is only of 30 % minus 1 t.

In order to compare the calculated solar productions to those of existing bakeries, the yearly productions of reference bakeries (subsection 2.6.3) are reported on the plot by mean of bakeries' names located on the axis at their respective range production. The solar production matches well the bakeries' production : Neoloco and LPN level of production (their work strategy is based on one baking day per week) can be reached with the weekly strategy. The level of production of ACDL (3 baking days per week) and Cassava (5 baking days per week) can be reached with the maximum production strategy.

Solar baking socio-technical potential

The soft mobility context is the same as for the "max production strategy" and the spatial distribution of both the EWSB and SBEI indexes over Réunion is depicted in Fig. 8. Fig. 8a shows the mapping of the EWSB index for the weekly strategy. With the weekly strategy, a lower amount of bread is produced (see the the solar bread production potential, Fig. 6) and an increase of the EWSB is mechanically expected. This is the case: areas where the EWSB index was the highest for the max production strategy still stand for higher values according to the weekly strategy, and the low value areas of the coast in the "max production strategy" show higher values in the "weekly operating strategy". Globally, with the "weekly strategy", all solar bread baked along the inhabited coast of Réunion should eventually be consumed.

Fig. 8b shows the SBEI index's map for the weekly strategy. Compared to the maximum production strategy map (Fig. 5b), the decline of the SBEI index between both strategies in areas where it was formerly the most significant is only 30 % in average (from 100-115 k € in Fig. 5b, to 64-81 k € in Fig. 8b). Globally, the most suitable area lie in the same color level and are the same as for the "max production strategy": the eastern part of Saint-Denis, cities of Saint-Pierre, Saint-Andre and Saint-Benoit. But even more interesting, north-east areas that had intermediate suitability in the "maximum production strategy" (in orange in Fig. 5b) are now the best areas (in red) in the weekly strategy, following a decline of only 20 % in average (80-100 k€ range to 64-81 k€ range) or less on some points. Fig. 9 plots the EWSB index versus the solar bread production for all the hexagonal cells and for both the max and weekly production strategies (respectively blue and red dots). For the max production strategy, the sites that produce more than 15 t or less than 5 t exhibit EWSB index values close to zero: they are irrelevant for solar baking. For the weekly strategy, irrelevant sites for solar baking are those corresponding to a production higher than 10 t or lower than 2 t. Hence we assume that the most robust organization for the bakery sector should consider a work strategy with site production between 5 and 10 t per year (with the Lytefire 11 oven). This range of production includes the production of Neoloco and LPN. The production level of ACDL can be reached with the max production strategy in some areas showing a high EWSB index. The production level of Cassava cannot be reached in Réunion when produced bread is consumed in close vicinity (EWSB index close to zero), for the sustainable mobility scenario.

Based on the average energy consumption for baking introduced in subsection 2.6.3 (≈ 26000 kWh/year for electric ovens) and the average ratio of direct emissions per electrical kWh consumed in 2019 in Réunion (719 gCO 2 /kWh) [START_REF] Gilboire | Bilan énergétique de la réunion 2019[END_REF], a rough estimation of of the saved gCO 2 emissions with a solar weekly strategy can be performed: considering that 26000 kWh/year corresponds to 6 electrical baking days per week, the emission of a classical "electrical" french bakery in Réunion with a 2 days weekly strategy is 719 × 26000/3 ≈ 6.2 tCO 2 .

In conclusion, the results of a simple weekly strategy show that both local optimization of the solar baker's time management and global optimization of the solar bakery sector can be performed in order to exploit the maximum of the solar bakery potential. For this, technology and workflow have to adapt to the solar resource variability in time and space, instead of adapting the energy production system to the technology and workflow already in place. In the weekly strategy, 5 days a week are released for other activities and the whole coastal area of the island (except the south-east volcano) becomes suitable

Conclusion

This article reports a methodology based on a socio-technical approach to a direct solar economy. First, a solar bakery workflow model was developed for assessing solar bread production based on 1) the solar energy resource, 2) the solar oven technology and 3) the baker's time management. This model allows the solar bread production potential of a territory to be estimated, and was applied to Réunion island. Second, we defined two socio-economical indexes that identify the best locations for solar bakeries based on the radical sustainability perspective of slow mobility. This approach is based on a shift from the traditional view of energy transition as an energy-centered concept, to a focus on the energy -technology -society system. Based on this socio-technical view of the energy transition, the corresponding solar bakery workflow model shows how it is not energy which is eventually produced, but baking time for bread production. In this way, even if the direct solar baking is a very sustainable option for baking and energy consumption management, it can hardly be taken into account in the energy mix calculation. Similarly, it is not rigorously possible to calculate energy savings as it was shown how an appropriate implementation of the direct solar bakery technology would impact both the baker's time management and the organization of the whole bakery sector: the shift from a part of the electrical bakery sector to direct solar baking cannot be reduced to the sole replacement of the electrical oven by direct solar ovens.

That is why the EWSB and SBEI indexes were calculated for a slow mobility scenario. However, even in this case, one is still looking at energy transition as a substitution of fuel tap, while the results strongly suggest that consuming solar bread would necessarily be different from the current consumption of electrical bread. If the heat energy fuel tap is only substituted, one implicitly considers that a solar bakery can deliver the exact same service as a current classical bakery. This is very unlikely, precisely because of the inherent volatile nature of direct solar radiation against the dispatchable nature of electrical supply. In order to look beyond sole substitution, one must understand how changes in bread consumer's behavior could meet the new paradigm of the solar baker's activity.

As an example, in Réunion, by considering one or two full bakeable days a week, i.e. with at least 6 continuous hours for bread production, a solar bakery could significantly enhance its viability. But it would require to bake a different type of bread on those specific days, a bread that could be stored and consumed during a greater period of time. In some way, this solar bread would serve as a bio-chemical storage of the solar radiation. It could address a current need for long-term bread storage, as well as a new manner of living that would have been co-constructed between the baker and its community of consumers. Different groups of consumers may be considered even within the same community, whose behavior would match the different bakery's ways and days of operating. Short-term or long-term bread consumption, and conservation, emphasize how solar radiation variability might be handled by the solar baker, without even considering another dispatchable backup resource. Apprehending this social time (see Fig. 3) is fundamental for ensuring a win-win between both production time and consumption time, and eventually promoting solar bakeries over electrical ones. This aspect sheds light on the social-cultural dimension of energy transition and on how it ultimately affects energy service sustainability.

In conclusion, the direct solar economy is radically sustainable compared to any other energy systems, and even if it cannot be implemented for all crafts or industries, it is a very interesting socio-technical alternative for the bakery sector's transition towards sustainability. Geographical aggregation is common strategy to mitigate variability from VRE power systems across territories [START_REF] Jewell | Limits on cloud-induced fluctuation in photovoltaic generation[END_REF]. One may apply the same idea regarding the bakery sector when solar energy is considered the main supply. It would ask for an economic reorganization where artisan bakeries are economically interconnected, e.g. belong to the same corporation. In Reunion, it makes sense as meteorological conditions are extremely variable over short distances. For small areas, and associated with soft mobility, it would then be theoretically possible for one solar bakery to run and compensate the lack of solar energy of the other ones. Another potential solution is to use mobile solar ovens pulled on trailers, so the production site can follow the sun. To keep the solutions sustainable, it must obviously be associated with short distances and soft mobility. Also, as with seasonal fruits, bread could be considered to be a product that is produced and consumed only at specific times during the year, though weekly or biweekly strategies ensure bread's availability almost over the entire year.

Figure 1 :

 1 Figure 1: Population density of the Réunion Island

2. 2 .

 2 The solar resource 2.2.1. Direct normal solar radiation

Figure 2 :

 2 Figure 2: Daily mean shaded direct normal irradiation (in kWh/m²) for the 2017-2019 period in Réunion

Figure 3 :

 3 Figure 3: Direct solar bakery's process workflow. Numeric filters' values in the purple boxes are given based on the Lytefire 11 solar oven model.

Figure 4 :

 4 Figure 4: Solar bread production thresholds in Réunion following the maximum production strategy. (a) Max number of bakeable days per year, (b) Average production shortage per year, (c) Max number of batches per year, and (d) Max yearly solar bread production.

Figure 5 :

 5 Figure 5: (a) Eaten Where Solar Baked and (b) Solar Bread Economic Intensity indicators in Réunion

  production for Max production strategy (t) Mass production for Weekly strategy (t) Maximum theoretical production

Figure 7 :

 7 Figure 7: Mass production for the weekly strategy vs mass production for the max production strategy (black dots). Colored dotted lines are empirical references (see text) and the grey dotted line shows the maximum theoretical production level with the weekly strategy.

Figure 8 :

 8 Figure 8: (a) Eaten Where Solar Baked and (b) Solar Bread Economic Intensity index for a weekly operating strategy in Réunion

Figure 9 :

 9 Figure 9: EWSB index vs production for both strategies

:

  Neoloco is a solar baker located in Europe (Normandie, France). Le pain nourricier (LPN), Au coin du levain (ACDL) and Cassava are bakeries located in Réunion. All of them bake sourdough breads.

	Name	Location	Fuel	Sales model	Work strategy	Yearly production 1 Price range
	Neoloco	Normandie	Sun or Wood Pre-sale	(baking days/week) 1	(tons/year) 7-8	( €/kg) 6-10
		(France)		& bread deposit
	Le pain nourricier Entre-Deux	Wood	Market place	2	6	6-12
	Au coin du levain	Sainte-Marie Electricity	Shop	3	12-15	8-15
	Cassava	Saint-Paul	Electricity	Shop	5	16-20	-

1

Yearly production is calculated from daily or weekly with a year equal to 50 weeks

Table 1 :

 1 Reference bakeries data

Acknowledgements

The authors would like to thank the CNRS and the IRD for their financial support through the MITI interdisciplinary programs.