Y Lee

M Bescond
email: marc.bescond@im2np.fr

N Cavassilas

D Logoteta

L Raymond

M Lannoo

M Luisier

Quantum treatment of phonon scattering for three-dimensional atomistic transport modeling

Based on the non-equilibrium Green's function (NEGF) formalism, we show a numerically efficient method to treat inelastic scattering in multi-dimensional atomistic codes. Using a simple rescaling approach, we detail the calculations of the lowest order approximation (LOA) 1 series to the usual, computationally intensive, self-consistent Born approximation (SCBA). This, combined with the analytic continuation technique of Padé approximants, is applied to an atomistic code based on a tight-binding sp 3 d 5 s * model for electrons and holes, and a modified valence-force-field method for phonons. Currents in Si and Ge gate-all-around nanowire transistors are then computed considering the main crystallographic transport directions ( 100 , 110 , 111 ) for both n-type and p-type devices. Our results show that in most configurations, 3 rd order LOA currents are enough to achieve a high agreement with SCBA results, while reducing the calculation time by about one order. In addition, we propose a criterion to determine the validity of such expansion techniques.

The recent introduction of three-dimensional (3D) Fin structured tri-gate transistors in the semiconductor industry has opened the door for a steady continuation of Moore's scaling law in the sub-20 nm technology node 2 . However, despite the success of the device architecture for the mass production, there is a remaining question: Will this device structure be still valid for transistors with sub-10 nm gate lengths? In this perspective, the ongoing researches mainly focus on i) improving the electrostatic control of the device architecture like gate-allaround (GAA) Si nanowire [START_REF] Mertens | Symp. VLSI Technol[END_REF] , or ii) searching novel materials exhibiting high mobilities and/or high band-gap [START_REF] Li | [END_REF][5][6][7][8][9] , or iii) the combination of both aspects, like Wrap-Gated (WG) InGaAs III-V nanowire transistors [10][11][12] .

In all the aforementioned devices both theoretical [13][14][15][16] and experimental 17 works underline the impact of the position of each individual atom on the obtained performance. Indeed, physics of nano-scaled devices is mainly governed by quantum effects (confinement 18 and tunneling 19,20 ), surface roughness scattering 21 and electron-phonon interactions 22 . Among those phenomena, it has been theoretically predicted 9,23 that inelastic interactions between electrons and phonons play a major role even in ultra scaled devices with dimensions of the order of the electron mean free path. Thus, in order to analyze the physical properties of the devices, the development of atomistic quantum simulation tools including inelastic electron-phonon scattering is urgently needed. During recent decades, several quantum methods have been proposed to treat this issue precisely. Among them the non-equilibrium Green's function (NEGF) formalism [START_REF] Datta | Electronic Transport in Mesoscopic Systems[END_REF][START_REF] Ferry | Transport in Nanostructures[END_REF][START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF] is one of the most advanced quantum methods that has attracted intensive interests due to its suitability for addressing inelastic carrier transport in nanostructures.

Unfortunately, atomistic NEGF quantum transport simulations of realistic devices require important numerical resources. The situation becomes even more complicated when inelastic interactions are included. In-deed phonon scattering is usually treated with the selfconsistent Born approximation (SCBA) that induces an additional self-consistent scheme to the conventional "Poisson-Schrödinger" loop. The SCBA was first applied to simple effective mass or k • p Hamiltonians 27,28 . It was then implemented in 3D full-band atomistic NEGF codes based on tight-binding Hamiltonian 22,29 . Among those advanced codes, the simulator OMEN developed by Luisier and co-workers is one of the most sophisticated that couples inelastic transport of phonons and electrons with an atomistic precision 30,31 . However, this simulator is only manageable with more than 1000 CPUs due to unavoidable numerical burden.

Recently we proposed a technique that allows to highly reduce the computational complexity for taking into account the interaction between electrons and phonons within the NEGF framework. The approach 1,32,33 is based on a lowest order approximation (LOA) of the interactions combined with an analytic continuation technique of Padé approximants [START_REF] Hardy | Divergent Series[END_REF][START_REF] Baker | Padé Approximants[END_REF] . Its main advantage is to avoid the numerous iterations required by SCBA. However, it implies the inversions and multiplications of potentially large matrices which may be impractical for 3D systems. We here present a simple rescaling technique able to calculate a series of LOA using the first iterations of the conventional SCBA algorithm. We apply this efficient technique to the atomistic code OMEN 30 to model 3D nanowire (NW) transistors for which electronphonon scattering is important. In particular, we show that the physical quantities (i.e. the currents) based on the 3 rd order LOA combined with Padé analytic continuation technique are enough to obtain the results of the full SCBA iterations. Our investigations focus on both ntype and p-type Si GAA NW considering the three main crystallographic orientations ( 100 , 110 , 111 ). Application to Ge NW in 110 transport direction is also discussed.

The device structure investigated in this paper is schematized in Fig. 1. The chosen gate length is 5 nm since it is considered as an ultimate limit for conventional logic transistors. The channel is surrounded by a 3 nm thick HfO 2 dielectric layer (with R = 20), while source and drain regions are covered by low κ dielectric layer (with R = 5). A supply voltage

| V DS | = | V DD | = 0.
6 V is applied between the source and drain contacts, and a gate-to-source bias V GS then controls the currents flowing inside the channel. Since electron-phonon scattering is found to be stronger in the saturation regime 22 , we focus here on this operation mode. We consider steady-state electron transport described within the NEGF framework employing a full band tightbinding sp 3 d 5 s * model without spin-orbit coupling [START_REF] Slater | [END_REF]37 . The lateral surface of the semiconductor NW in which we solve the Schrödinger's equation is passivated by increasing the energy of dangling bonds 38 , and no surface roughness scattering is considered. Phonon relation dispersion in the confined NW is obtained via a modified valenceforce-field method including four interaction terms 39,40 . Here, as described in Ref. 30, electrons are coupled to phonons through the following diagonal lesser or greater self-energy

Σ ≶σ1σ2 nn (E) = i l ij σ3σ4 ∞ -∞ d( w) 2π ∇ i H σ1σ3 nl × G ≶σ3σ4 ll (E ∓ ω)∇ j H σ4σ2 ln D ≶ij ln (ω) -D ≶ij ll (ω) -D ≶ij nn (ω) + D ≶ij nl (ω) , (1)
where is the reduced Planck's constant, E the electron energy, ω the phonon frequency, H nl the Hamiltonian (tight-binding matrix element), G(E) the electron Green's function, and D(ω) the phonon Green's function expressed in the real-space (the indices i and j run over x, y, and z directions). The other indices n and l denote all atomic positions, and σ represent the orbitals. ∇ i H nl represents the first derivative of the tight-binding Hamiltonian matrix H nl along the i direction (x, y, or z). It defines the electron-phonon coupling due to atom oscillations around their equilibrium position (please see Ref. [30] for more details). We assume a bath of phonons at thermal equilibrium by considering no temperature gradients inside the device and by neglecting the phonon self-energies which can contribute to the phonon Green's function D(ω) (see Ref. 30). We should note that the diagonal approximation of the self-energy, which results from computational limitations in atomistic codes, may lead to an underestimation of the strength of electronphonon coupling 41,42 . However, Ref. 43 showed that the diagonal approximation was still valid in ultra-scaled Si nanowires for electrons, except for transport along the 110 crystal orientation. More generally, the LOA combined with Padé approximant can be applied to nondiagonal self-energies. The self-energy of Eq.( 1) combined with Dyson's equation [START_REF] Mahan | Many-Particle Physics[END_REF] allows to calculate the electron interacting Green's function G(E) starting from the non-interacting one g(E).

G = g + gΣ[G]G, (2)
where the simplified matrix notation

[1 = (r 1 , t 1 ), gΣ[G]G = d2 d2 g(1;
2)Σ(2; 2 )G(2 ; 1 )] is used. Dyson's equation ( 2) is typically solved using the socalled iterative SCBA scheme with a specified convergence tolerance such as

G N = [g -1 -Σ[G N -1 ]] -1 (3)
where G N is the Green's function at the N th SCBA iteration and G 0 = g. Green's functions resulting from first few iterations (i.e. N small in Eq.( 3)) of the SCBA scheme are generally not conserving for strong electron-phonon interactions 1,33 . The SCBA might then require more than several dozens of iterations to provide accurate conserving physical properties. On the other hand, LOA Green's functions are conserving at any order. It has been shown that the first few orders of LOA Green's functions coupled with Padé approximants can highly accelerate the calculations of physical quantities (both current and charge density) 1 . However, the LOA Green's function algorithm presented in Ref. 1 implies numerous inversions and multiplications of Green's function matrices. Its implementation is then numerically very demanding when applied to 3D nanostructures. Here, we show that the derivation of LOA Green's functions from the usual SCBA algorithm can be easily reached by introducing a scaling parameter to Eq.(1).

As an example, let us rewrite Eq.(3) in the Taylor series expansion as

G N = g + gΣ[G N -1 ]g + gΣ[G N -1 ]gΣ[G N -1 ]g + • • • . (4)
The first SCBA Green's function G 1 is then defined from Eq.( 4) as

G 1 = g + gΣ[g]g + gΣ[g]gΣ[g]g + • • • . ( 5
) (a) G λ1 1 = + 1 λ 1 + 1 λ 1 2 Conserving terms + 1 λ 1 3 +• • • Suppressed non-conserving higher order terms of G λ1 1 (b) G λ2 2 = + 1 λ 2 + 1 λ 2 2 + Conserving terms + 1 λ 2 3 + + +• • •
Suppressed non-conserving higher order terms of G G 1 is not conserving due to second and higher order terms in interaction 33 . Its first order LOA counterpart G LOA1 is conserving and expressed as

G LOA1 = g + gΣ[g]g. (6)
G LOA1 can then be deduced from G 1 by introducing a scaling parameter 1/λ 1 in the self-energy of Eq.( 1). G 1 rescaled by λ 1 is now defined by:

G λ1 1 = g + 1 λ 1 gΣ[g]g + 1 λ 1 2 gΣ[g]gΣ[g]g + • • • . (7)
Taking a relevant value of λ 1 , the factors (1/λ 1 ) n (n > 1) can eliminate the infinite sum of non-conserving terms while preserving the first order conserving one. In NEGF, since any relevant expectation value (O) is linear in the one-electron Green's function, the expectation value of the first order term to G LOA1 can be constructed as

∆O 1 = O(gΣ[g]g) = λ 1 [O(G λ1 1 ) -O(g)].
The same technique can be applied to the second SCBA iteration of Eq.( 4) considering a different scaling factor λ 2 to obtain the rescaled G 2 as

G λ2 2 = g + 1 λ 2 gΣ[G λ2 1 ]g + 1 λ 2 2 gΣ[G λ2 1 ]gΣ[G λ2 1 ]g + • • • , (8)
where

G λ2 1 is G 1 rescaled by λ 2 .
Since the SCBA selfenergy is linear with respect to the Green's function 33 , the rescaled G 2 can be expressed as

G λ2 2 = g + 1 λ 2 gΣ[g]g + 1 λ 2 2 (gΣ[g]gΣ[g]g + gΣ [gΣ [g] g] g) + 1 λ 2 3 (gΣ[g]gΣ[g]gΣ[g]g + • • • ) + • • • . ( 9
)
where λ 2 should be chosen to remove higher order terms associated to (1/λ 2 ) n (n > 2), but to maintain the first and second order conserving terms. With this scaling parameter, the expectation value of the second order terms to second order LOA Green's function G LOA2 can be obtained as Once the rescaled conserving Green's functions at each order are obtained, the expectation values of the desired observable can be reconstructed using the linearity property of the corresponding operator. In particular, N th order LOA current can be expressed as

∆O 2 = O(gΣ[g]gΣ[g]g + gΣ [gΣ [g] g] g) = λ 2 2 [O(G λ2 2 ) -O(g)] -λ 2 ∆O 1 . Figures 2 (a
I N th LOA = I 0 + N n=1 I n , (10)
with

I N = (λ N ) N I N th SCBA (G λ N N ) - N -1 n=0 (λ N ) N -n I n , I 0 = I 0 . (11)
This procedure demonstrates that conserving Green's function at any order N can be in principle obtained from conventional SCBA scheme. At each order N , λ N has to be large enough to remove all non-conserving terms (whose order in interaction is larger than N ) but small enough to maintain conserving terms (whose order in interaction is ≤ N ). This approach relies on the relevant choice of λ N which can be rapidly verified through the current conservation. Explicitly we determine the λ N 's based on a semi-empirical method noting that they are associated to a Taylor development of the Green's function. We first find the scaling parameter λ 1 , which ensures that the current is conserved after the first SCBA iteration. We then deduce an approximated value of λ 2 as 2 × (λ 1 ) 1/2 . Based on the same argument, we take λ 3 equal to 5 × (λ 1 ) 1/3 . The multiplication factors 2 and 5 for λ 2 and λ 3 , respectively represent the weighting according to the number of diagrams at each order. We are aware of the empirical character of this approach, but the assessment of scaling parameters at any order in interactions and to any systems is a very complex issue beyond the scope of the present work. Here, we calculate up to 3 rd order currents and apply Padé approximants 0/1, 1/1, and 1/2. Calculation of higher order LOA currents is also possible, but it usually requires additional SCBA iterations that could reduce the relevance of the presented technique. Moreover, the approach could have been also applied to the calculation of the charge carrier density by using the same scaling parameters as those for the current. For the sake of clarity, we decided to only focus on the current considerations. The comparisons between SCBA and LOA/Padé approximant are performed with the same electrostatic potential profile as the one resulting from the self-consistent loop between the Poisson equation and the NEGF formalism including the interactions within the SCBA scheme.

Table I shows comparisons of the first three LOA orders and corresponding Padé currents with SCBA values for n-Si (n-type silicon) and n-Ge (n-type germanium) devices along several principal crystallographic orientations. We first see that all LOA current series (LOA1, LOA2, and LOA3) are diverging from SCBA currents, meaning that the considered devices operate in a relatively strong electron-phonon interaction regime 1 . Padé 0/1 applied to the first order LOA currents shows an underestimation of the current degradation compared with the SCBA currents, except for n-Si 111 -oriented device. In addition, Padé 1/1 approximants based on up to 2 nd order LOA currents are in most cases not reliable as pointed out in previous work 1 . However, the application of Padé 1/2 on LOA current series shows significant improvements, with high accuracies compared with the SCBA values. For n-Si and n-Ge 110 -oriented devices, Padé 1/2 provides almost the same values as the SCBA currents with a relative difference (ε) of less than 2% on average. It is also worth noting that Padé 1/2 can provide a reasonably good result (< 9% relative difference) for 111 transport direction where SCBA needs more than 40 iterations to converge. Therefore, currents resulting from Padé 1/2 can be seen as a good estima-tion to the SCBA values with a relative difference less than 10% for all the considered n-type devices. These results confirm the previous work 1 where it has been shown that Padé N/N + 1 can give better results even in strong electron-phonon interacting system.

Table II shows similar results as Table I for p-type devices. We see that all LOA currents (LOA1, LOA2, and LOA3) are divergent series like in the case of n-type devices. Table II also reports very favorable results for the 110 direction for which even the simplest Padé 0/1 guarantees a high agreement of the current (ε < 6%) with respect to the SCBA approach. Padé 1/1 also produces reliable results for the 110 direction while remaining unstable for the other configurations. Padé 1/2 still gives current value very close to those of SCBA with ε less than 6%. The result is also remarkable for the 100 direction where Padé 1/2 provides similar accuracy while SCBA requires more than 80 iterations. However, the comparison is less successful for the p-Si 111 configuration for which all LOA and Padé currents fail to reproduce the SCBA result.

The transport along the 111 direction is quite singular for both n-type and p-type devices. For that orientation first order LOA currents are larger than the ballistic values (see Tables I andII). Concerning the n-type device the converged SCBA current is also larger than the ballistic one and the Padé 1/2 finally succeeds to provide an accurate value. In n-type 111 -oriented NW the counterintuitive increase of the current in the presence of phonon scattering with respect to the ballistic value results from the interplay between two competing effects: first, backscattering caused by phonon emission or absorption tends to decrease the current magnitude. Secondly, inter-subband scattering induced by electronphonon interactions connects bands with a small energy bandwidth that would otherwise not carry any current in the ballistic limit of transport [START_REF] Chang | [END_REF][START_REF] Szabó | Dissipative quantum transport simulations in two-dimensional semiconductor devices from first principles[END_REF] . In ultra-short devices with bandstructure exhibiting several sub-bands with a low energy bandwidth, as 111 -oriented Si nanowires, the second effect dominates, thus leading to an increase of the current. However, in longer devices, backscattering plays the major role and the situation goes back to nor- mal, i.e. the current magnitude decreases when electronphonon scattering is turned on. The case of p-type device is more pathological since the converged SCBA becomes smaller than the ballistic value. Such a configuration seems to involve more complex physical phenomena that require to consider higher order LOA currents. It clearly illustrates the limit of our technique to reproduce SCBA currents. We also mention that it has been recently shown [START_REF] Lee | [END_REF]48 that the SCBA with Luttinger-Ward Functional self-energy may produce unphysical convergence for strong-coupling system. 111 p-type direction could be one of these configurations. We can also note that a similar approach using Hypergeometric resummation has been recently proposed and could be tested in such case 49 . However, we suggest to take as a warning signal of the accuracy of the present approach the situation where first order LOA current is larger than the ballistic one.

In summary we have proposed a numerically efficient technique to reproduce the SCBA currents. It relies on the LOA coupled to Padé approximants and only requires the first few iterations of the SCBA procedure. The method has been applied to NW transistors using an atomistic 3D quantum transport code. We focused on the saturation regime current of the devices at high gate bias and compared the LOA and Padé currents with the SCBA values. In particular, we showed that the Padé 1/2 with only first three order LOA currents can successfully reproduce the SCBA currents with a small relative difference (< 10%) in most of considered devices except for p-Si 111 device where the electronphonon interaction exhibits singular effects. We also suggested as a warning signal to SCBA solution the case where first order LOA current is larger than ballistic one. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) and funded by the NOODLES contract (ANR-13-NANO-0009) of the ANR-French National Research Agency. The author Y. Lee would like to thank the AMIDEX for financial support.

FIG. 1 .

 1 FIG. 1. (Color online) Schematic representation of the silicon GAA NW transistors considered in this work (with a 100 crystallographic orientation in the present case). The length of the gate is LG = 5 nm, the source and drain extensions L S/D = 15 nm, and the diameter of NW dNW = 3 nm. Source and drain regions are doped with a doping concentration of donors or acceptors N D/A = 1 × 10 20 cm -3 . Si 110 , Si 111 , and Ge 110 are also considered (not shown).

λ2 2 FIG. 2 .

 22 FIG. 2. (Color online) Rescaled Green's functions, (a) λ 1 1 of first SCBA iteration and (b) G λ 2 2 of second SCBA iteration, described by Feynman diagrams. Conserving and non-conserving terms are arranged in ascending order in interaction. Thin lines with an arrow indicate electron non-interacting Green's function. Dashed lines are for free phonon propagators.

  ) and 2 (b) illustrate this scaling technique in terms of Feynman diagrams for the first and second SCBA iterations, respectively.

TABLE I .

 I Comparisons of 3 first LOA currents and corresponding Padé approximants (0/1, 1/1, and 1/2) with the ballistic and SCBA currents in n-Si 100 , n-Si 110 , n-Si 111 , and n-Ge 110 devices. Relative differences (ε) of the LOA and Padé currents with respect to the SCBA values are also reported (ε = 100 × |ISCBA -I|/ISCBA) where I is ballistic, LOA, or Padé current.

		Ballistic	LOA1	LOA2	LOA3	Padé 0/1	Padé 1/1	Padé 1/2	SCBA	Device
	Current [A]	8.12e-6	6.40e-6	-2.37e-5	2.47e-4	6.70e-6	8.22e-6	5.05e-6	5.26e-6	n-Si 100
	ε [%]	54.4	21.7	550.6	4595.8	27.4	56.3	4.0	0.0
	Current [A]	6.55e-6	4.58e-6	-3.92e-5	1.61e-3	5.04e-6	6.64e-6	4.42e-6	4.45e-6	n-Si 110
	ε [%]	47.2	2.9	980.9	36079.8	13.3	49.2	0.7	0.0
	Current [A]	4.34e-6	1.72e-5	-4.38e-5	2.03e-4	-2.20e-6	6.59e-6	5.48e-6	6.01e-6	n-Si 111
	ε [%]	27.8	186.2	828.8	3277.7	136.6	9.6	8.8	0.0
	Current [A]	5.27e-6	4.76e-6	2.09e-6	7.11e-5	4.80e-6	5.39e-6	4.72e-6	4.62e-6	n-Ge 110
	ε [%]	14.1	3.0	54.8	1439.0	3.9	16.7	2.2	0.0

TABLE II .

 II Comparisons of 3 first LOA currents and corresponding Padé approximants (0/1, 1/1, and 1/2) with the ballistic and SCBA currents in p-Si 100 , p-Si 110 , p-Si 111 , and p-Ge 110 devices. Relative differences (ε) of the LOA and Padé currents with respect to the SCBA values are also reported (ε = 100 × |ISCBA -I|/ISCBA) where I is ballistic, LOA, or Padé current.

		Ballistic	LOA1	LOA2	LOA3	Padé 0/1	Padé 1/1	Padé 1/2	SCBA	Device
	Current [A]	6.10e-6	1.96e-6	3.53e-8	-3.36e-5	3.63e-6	-1.65e-6	3.89e-6	4.31e-6	p-Si 100
	ε [%]	41.5	54.5	99.2	879.6	15.8	138.3	9.7	0.0
	Current [A]	6.26e-6	3.77e-6	5.07e-6	1.81e-5	4.48e-6	4.66e-6	4.48e-6	4.74e-6	p-Si 110
	ε [%]	32.1	20.5	7.0	281.9	5.5	1.7	5.5	0.0
	Current [A]	5.29e-6	8.18e-6	-4.12e-5	3.70e-3	1.16e-5	5.45e-6	9.08e-6	4.85e-6	p-Si 111
	ε [%]	9.1	68.7	949.5	76188.7	139.2	12.4	87.2	0.0
	Current [A]	5.80e-6	4.11e-6	4.70e-6	5.03e-7	4.49e-6	4.55e-6	4.49e-6	4.74e-6	p-Ge 110
	ε [%]	22.4	13.3	0.8	89.4	5.3	4.0	5.3	0.0


F Michelini
email: fabienne.michelini@im2np.fr

K Beltako

M Bescond

N Cavassilas

L Raymond

Time-resolved quantum transport for optoelectronics

We investigate time-resolved energy currents in a molecular optoelectronic junction made of two donors and an acceptor sandwiched between two electrodes and excited by a Gaussian femtosecond laser pulse. Features of direct energy currents are thus correlated to the intra-molecular structure.

I. INTRODUCTION

Recent experimental trend toward femtosecond multidynamical excitations has triggered an interest for investigating time-dependent regimes in nanostructured semiconductor, molecular or hybrid electronic systems [START_REF] Bakulin | Ultrafast Spectroscopy with photocurrent detection: Watching Excitonic Optoelectronic Systems at work[END_REF]. Molecular electronic devices form indeed a promising alternative to standard electronic devices due to their fast response on ultra-short time scale. The quest for understanding the time-resolved response of nanometer scale junctions interacting with external electromagnetic fields is thus fueled by both the fundamental questions it raises and the premise for long term applications in optoelectronics. In this delicate context, the problem of energy transfer is emerging as one of the most crucial issues at several levels, from the impasses meat by quantum thermodynamics [START_REF] Esposito | The nature of heat[END_REF] to biologically inspired high efficiency concepts [START_REF] Creatore | Efficient biologically insired photocell enhanced by delocalized states[END_REF]. In this paper, we derive and examine time-resolved energy currents induced by femtosecond laser pulses in a molecular optoelectronic junction sandwiched between metallic electrodes.

II. THEORETICAL ORIENTATIONS

By the end of the last century, nanodevices reached the nanometer scale, and lived a tremendous change of paradigm as they become quantum systems. Such a breakthrough has revealed powerful methods of quantum statistics for simulating nanowire transistor characteristics, in particular the method of non-equilibrium Green's functions (NEGFs) [START_REF] Stefanucci | Nonequilibrium Many-Body Theory of Quantum Systems: a Modern Introduction[END_REF]. Today, this formalism is widely used for several device functionalities, but challenges remain, especially for timedependent regimes and quantum thermoelectricity.

Few years ago, a technique was proposed by Gaury and co-workers in order to manage time-dependence in mesoscopic systems [START_REF] Gaury | Numerical simulations of time-resolved quantum electronics[END_REF]. Elementary objects of the technique are wave functions (WFs) instead of NEGFs, vectors instead of matrices. In fact, WFs are components of the retarded Green's function projected on propagation modes that are defined from reservoir self-energies. Despite being applied only in non-interacting cases so far, the technique seems also promising to tackle interactions. The technique aims to manage time-dependence by means of simpler simulations while accounting for the complete energy dependence of the self-energies, i.e. beyond the wide band limit.

On the other side, quantum thermoelectricity has recently emerged as a part of the solution to the problem of sustainable emission-free energy production. However, these advances have raised some fundamental questions about energy and heat [START_REF] Michelini | Entropy production in photovoltaic-thermoelectric nanodevices from the non-equilibrium Green's function formalism[END_REF], which are crucial for thermodynamics and energy conversion at the nanoscale especially in time-dependent regimes [START_REF] Esposito | The nature of heat[END_REF]. Morevor, few teams around the world develop specific time-resolved techniques to directly address the photo-carrier dynamics [START_REF] Bakulin | Ultrafast Spectroscopy with photocurrent detection: Watching Excitonic Optoelectronic Systems at work[END_REF], and even less the photo-carrier thermodynamics at the nanometer scale.

We present a practical methodology based on lattice Hamiltonian to investigate energy exchanges on ultimate time and size scales. The approach is developed for the study of a molecular optoelectronic junction that interacts with femtosecond laser pulses.

III. ENERGY CURRENT IN THE MOLECULAR JUNCTION

We consider an optoelectronic junction made of three subsystems: two donors (D 1 and D 2 ) interact with a Gaussian femtosecond laser pulse, and they are both connected to an acceptor (A). The whole complex is in contact with two thermal reservoirs. The junction is shown Fig. 1. Each donor is described by a two-level model of HOMO and LUMO states. For the acceptor, we use a two-site model with a single level model to study energy transport. All parameters are shown in the molecular complex level structure of Fig. 2.

For the loop circuit designed by the molecular complex, see Fig. 3, we derive an energy current operator from the continuity equation. In Ref. [START_REF] Wu | Energy lux operator, current conservation and the formal Fourrier's law[END_REF], the authors developed an energy current operator for a one-dimensional lattice Hamiltonian. We here address the loop case. The difficulty underlying the definition of an energy current operator resides in the fact that in contrast with the charge current operator for which: N = n s s ∑ , the total Hamiltonian of the system doesn't have the form of a sum of subsystem Hamiltonians.

Indeed H = h s 0 + V s, r ( ) r(s) ∑ s ∑ : all sites s have a local
Hamiltonian, h s 0 , but they are also interconnected with coupling terms V s, r ( ) = t sr a s + a r + h.c. , where t sr is the coupling parameter, t sr = 0 if sites s and r are not connected.

We have to redefine local Hamiltonians that include intersite coupling terms. We adopt the same proposal as [START_REF] Wu | Energy lux operator, current conservation and the formal Fourrier's law[END_REF] with

h s = h s 0 + 1 2 V s, r ( ) r(s) ∑ .
For this peculiar circuit, we identify four contributions to energy current operator by comparison with the onedimensional discrete case treated by Wu and Segal [START_REF] Wu | Energy lux operator, current conservation and the formal Fourrier's law[END_REF]. This viewpoint makes possible to distinguish between the two single line terms, two-site J E (2) and four-site J E (4) terms, the additional branch term, J E (B) , and finally two loop contribution in J E (L ) . All these contributions are schematically represented Fig. 3. Time-resolved mean values of these contributions and total energy current were calculated using WF technique [START_REF] Gaury | Numerical simulations of time-resolved quantum electronics[END_REF], still in the wide band limit.

IV. FINDINGS

In Fig. 4, the time-resolved energy current exhibits oscillations whose frequency varies as a function of the D 1 -D 2 coupling. These oscillations manifest the proper frequencies of the internal structure of the molecular complex. Indeed, the donor-donor coupling modifies the internal structure of the complex: frequencies change. As long as electrons have enough time to explore all the hybridized states in the system, they oscillate between these states before tunneling to reservoirs. To study the effects of couplings and excitation characteristics, we now rely on the direct energy current, that is the energy current integrated time.

The direct energy current together with its different components are plotted as a function of the D 1 -D 2 coupling in Fig. 5. The J E (2) term represents the current which is directly proportional to charge current. This contribution is the major one. The J E (4) term represents the four-site contribution along the D1-A line, while J E (B) involves the D 2 -A branch. However, term J E (B) is several orders of magnitude lower than J E (4) because carriers cannot reach the left reservoir via this branch (tunneling parameter is zero for left LUMO states). Finally, the loop term, J E (L ) , provides a negative contribution to the total energy current. The D 1 -D 2 coupling drives energy in the loop formed by the two donors and the acceptor. Increasing simultaneously both D-A and Γ couplings naturally emphasize both J E (4) and J E (L ) contributions (see Fig. 3 for schematic definitions).

We finally explore impacts on energy current of a dissymmetry between the two donors.

The first difference is introduced via excitation. We suppose that the two laser pulses are not synchronized on a femtosecond timescale, a delay occurs between the two lightdonor interactions. For this study, two molecule configurations are considered: the case where the donors are not connected, and the loop in its highly symmetric configuration, where the two donors are connected with a coupling parameter equal to the D-A coupling value, both cases are shown Fig. 6. While a monotonous behavior of both currents is observed in the first case, in the loop case the two currents flowing from D 1 and D 2 towards A are oscillating with an opposite phase and a frequency about D-A coupling.

The second difference between the two donors is implemented in the D-A coupling parameter. The D 2 -A coupling value increases from zero to twice the D 1 -A one. In the D 1 -A branch, energy current decreases while it increases starting from zero in the D 2 -A branch, as plotted Fig. 7. Interestingly, the sum of the two currents exhibits a maximum close to the symmetric configuration also shown Fig. 7.

At last, the roles of D-D and D-A couplings are exchanged: the D 2 -A coupling is zero while the D 1 -D 2 coupling increases from zero to twice the D 1 -A coupling value. We observe two remarkable behaviors on Fig. 8. On the one hand, J D1-A E increases when D 1 is connected to D 2 , and on the other hand, J D2-A E is non-zero while the two sites D 2 and A are not connected. This intriguing last result is originating from one of the J E (L ) components. Indeed, from the definition we have chosen, J E (L ) includes a term which is not proportional to the coupling between D 2 and A (see Fig. 8). What is the meaning of this contribution? We believe that it is likely related to transport channels activated by the time-dependent electromagnetic field. Otherwise, the energy current operator shall be accordingly adapted.

V. CONCLUSION

We studied the dynamics of energy transfer at donor-acceptor interfaces in a molecular optoelectronic junction using a mean-field approach inside the wave function technique. This exploratory work still ignores electron-electron interaction and vibrations, but nonetheless provides basic insights towards the design of new energy-conversion schemes for time-dependent optoelectronics.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Schematics of the molecular complex junction. The Gaussian pulse has a central frequency equal to the donor gap (1.4 eV) , a width of 30 fs and an amplitude of 1 eV.

Fig. 4 .

 4 Fig. 4. Time-resolved energy current as a function of the donordonor coupling.

Fig. 5 .

 5 Fig. 5. Direct energy current as a function of the donor-donor coupling for two values of D-A coupling. All components of energy current are shown.

Fig. 6 .

 6 Fig. 6. Direct energy current as a function of the delay between the two femtosecond laser pulses for two D-D coupling values.

Fig. 7 .

 7 Fig. 7. Direct energy currents, , and their sum, as a function of the D 2 -A coupling.

Fig. 8 .

 8 Fig. 8. Direct energy current as a function of the D-D coupling when sites D 2 and A are disconnected. Schematics of inter-site couplings involved in the component of .

ACKNOWLEDGMENT

The authors acknowledge financial support from the ANR-French National Research Agency [NOODLES, contract No. ANR-13-NANO-0009].


Benoit Galvani

Fabienne Michelini

Marc Bescond

Masakazu Sugiyama

Jean-François Guillemoles

Nicolas Cavassilas

Minibands modeling in strained balanced InGaAs/GaAs/GaAsP cells

Keywords: multi quantum wells, quantum modeling, minibands, carrier collection

. This approach is efficient thanks to the strain-balanced materials which, at a well under compressive strain, associates a barrier under tensile strain. This permits to consider a large number of wells while preventing the formation of dislocations during crystal growth.

. Finally the recombination rate is large and impacts both open-circuit voltage and shortcircuit current.

. The wave functions of carriers in minibands are Bloch waves, meaning that propagation is efficient. Our theoretical study, based on quantum simulation (Green functions formalism) in InGaAs/GaAs/GaAsP cells, sheds light on minibands in which the average velocity of carriers is around 10 7 cm s -1 . However, we also show that, without an adapted design, such minibands are inefficient since they connect only a few wells. We will present some improvements related to the distance between barriers and the positioning of the MQW inside the cell.

INTRODUCTION

Over the last decades, carriers transport phenomena in semiconductor materials have been investigated to increase the performances of optoelectronic devices. The study was focused on properties of materials used in such devices. Then, Leo Esaki and Raphael Tsu have proposed in 1970 a theoretical study [START_REF] Esaki | Superlattice and Negative Differential Conductivity in semiconductors[END_REF] of structures with periodic potential called superlattices (SL). They showed that SLs could exhibit remarkable optical and electronical properties, such as negative differential conductance and minibands. Even if the physical concepts predicted by Esaki and Tsu were controversial, the study of these quantum mechanical effects has opened a brand new research field in semiconductors. The development of epitaxial growth techniques has permit to fabricate SLs-based devices made of layers of a few nanometers, and thus to validate the theoretical models experimentally. In 1986, Capasso et al. observed for the first time resonant tunneling in a AlInAs/GaInAs SL [START_REF] Capasso | Sequential resonant tunneling through a multiquantum well superlattice[END_REF]. Miniband conduction was also investigated by Schneider et al. in a GaAs/AlAs structure [START_REF] Schneider | Tunneling resonances and miniband conduction in superlattices[END_REF]. Mendez et al. worked on Stark localization in GaAs/GaAlAs SLs to highlight the influence of an applied electric field on the miniband conduction [START_REF] Mendez | Stark Localization in GaAs-GaAlAs Superlattices under an Electric Field[END_REF]. During the early 1990's, a lot more experiments [START_REF] Sibille | Limits of semiclassical transport in narrow miniband GaAs/AlAs superlattices[END_REF][START_REF] Mendez | Wannier-Stark Ladders and Bloch Oscillations in Superlattices[END_REF] were carried out to observe optical and electronical properties of SLs, even at room temperature [START_REF] Dekorsy | Bloch oscillations at room temperature[END_REF]. Possible applications for these interesting properties were discussed [START_REF] Esaki | Superlattice and Negative Differential Conductivity in semiconductors[END_REF][START_REF] Mendez | Wannier-Stark Ladders and Bloch Oscillations in Superlattices[END_REF]. Capasso et al. observed sequential resonant tunneling in SLs that could be used for lasers and wavelength selective detectors [START_REF] Capasso | Resonant Tunneling Through Double Barriers, Perpendicular Quantum Transport Phenomena in Superlattices, and Their Device Applications[END_REF]. Rudolf F. Kazarinov and Robert A. Suris suggest the possibility of using intersubbands transitions in a semiconductor device made of SLs [START_REF] Kazarinov | Possibility of the amplification of electromagnetic waves in a semiconductor with a superlattice[END_REF] to achieve laser emission, which will lead to the creation of quantum cascade lasers. Most of SLs devices are light emitting devices such as lasers, where SLs dimensions and materials can be tuned to modify the band structure and thus to control the wavelength of emitted photons. Frequency modulators exploit the negative differential conductance property of semiconductor SLs to operate at high frequency. Solar cells can also benefit from the use of MQW structures to tailor the optical absorption [START_REF] Browne | Triple-junction quantum-well solar cells in commercial production[END_REF][START_REF] Adams | Recent results for single-junction and tandem quantum well solar cells[END_REF].

In this paper, we discuss the possibility of using minibands in SLs to improve the carrier transport in solar cells and thus to increase their efficiency. Our numerical study was carried out on a InGaAs/GaAs/GaAsP MQW nanostructure used in p-i-n single junction solar cells. We have implemented a numerical code using Green Function formalism [START_REF] Cavassilas | Modeling of nanoscale solar cells: The Green's function formalism[END_REF] in order to calculate the band diagram of the whole structure and the corresponding local density of states (LDOS). We have studied the influence of SL layer widths on both minibands formation and carrier transport. We figured out that some structural optimizations can allow to extend both valence and conduction minibands across the whole MQW region, and thus to improve the transport and the collection of carriers.

THEORY OF MINIBANDS FORMATION

Before showing our results, we will briefly explain the physical processes behind minibands formation in SLs. Let's consider a periodic SL made of N periods with two different semi-conductor materials A and B. The energy gap difference between A and B will form a barrier/well potential configuration. The Hamiltonian of this system can be written using an N-by-N matrix representation shown below :

where E 1 represents the eigenvalue of the fundamental energy level for an isolated well and V represents the tunnel coupling terms between wells (hopping). We consider in this case nearest-neighbor interactions only. It have been demonstrated that the E(k) dispersion relation of electrons in this system is a sinusoidal function in the form:

In this equation, k is the wave vector and a is the SL period. The dispersion relation of this system defines allowed energy states called miniband. This phenomenon will modify the carrier transport properties of the SL since electrons (and holes) are delocalized in the whole region crossed by the miniband. From Eq. ( 2), we can calculate the effective mass of electrons in the miniband :

(1)

(2)

The effective mass depends on the SL period a and on the hopping term V. The average velocity can also be obtained from :

Minibands in SL could potentially increase the current generated by the cell. Indeed, if a carrier reaches one of the energy states composing the miniband, it will be completely delocalized along the whole structure, meaning that such carriers propagate with an average velocity which can be high.

STRUCTURE AND PARAMETERS

The MQW region of the studied solar cell is an undoped InGaAs/GaAs/GaAsP SL (20 periods) made of alternation of 22 Å GaAs 0.58 P 0.42 barrier, 38 Å In 0.21 Ga 0.79 As well and two 31 Å GaAs interlayers between each barrier and well. The SL is sandwiched between two undoped GaAs layers, with an adapted thickness to get a total undoped region (GaAs layers + SL) of 1µm. GaAs interlayers helps to reduce lattice relaxation during epitaxial growth [START_REF] Fujii | 100-period[END_REF]. The whole calculated band diagram is shown in Fig. 1. Self-consistent calculations were made by resolving Schrödinger-Poisson equations using Green function formalism. This way permits to assume quantum properties of electronic transport such as tunneling, interference, and confinement.

The energy band gaps of each material with such compositions are 1.42 eV for GaAs, 1.99 eV for GaAsP and 1.13 eV for InGaAs. The effective masses of electrons are equal to 0.067 m 0 for GaAs, 0.1102 m 0 for GaAsP and 0.0585 m 0 for InGaAs, where m 0 is the free electron mass [START_REF] Vurgaftman | Band parameters for III-V compound semiconductors and their alloys[END_REF]. However, in strained materials, the band structure will be modified, which will lead to some differences in the values of energy gaps and effective masses. In order to use the correct values as entry parameters, we have calculated the band structure of each strained material (InGaAs and GaAsP in our system) using a 8-band k.p model. Calculated energy gaps in strained materials are equal to 1.87 eV and 1.15 eV for GaAsP and InGaAs respectively. Effective masses of electrons were also calculated and are equal to 0.097 m 0 and 0.068 m 0 for GaAsP and InGaAs respectively. (3) (4)

MINIBAND WITHOUT ELECTRIC FIELD

For a first calculation, we assume no electric field in order to extract the minibands parameters. Our numerical LDOS calculations for a 20-period SL as defined above are shown in Fig. 2 as a function of the energy and the position in the transport direction. LDOS shows minibands in the GaAsP barriers in both conduction band (CB) and valence band (VB). Such minibands are extended across the whole MQW region. We note that the miniband in CB is made of second-order energy states. The first-order energy states are located deeper inside the well, but they are not visible here. Indeed, in this calculation, we do not assume inelastic behavior such as electronphonon or electron-photon scatterings. In this case, states in the QWs are no coupled to the contacts, and the corresponding LDOS is not visible [17,[START_REF] Cavassilas | Theoretical comparison of multiple quantum wells and thicklayer designs in InGaN/GaN solar cells[END_REF].

For the minibands, which are directly coupled to the contacts, the effective masses for CB and VB calculated with Eq. ( 3) are equal to 0.029 m 0 and 0.173 m 0 respectively. These results show that both electrons and holes have effective masses lower than those in a bulk GaAs material, which are equal to 0.063 m 0 for electrons and 0.51 m 0 for holes. Average velocities calculated from Eq. ( 4) are equal to 3.22 10 7 cm/s for electrons and 5.52 10 6 cm/s for holes. These velocities are much higher than velocities in a SL without minibands, estimated to be 2.1 10 4 cm/s for both electrons and holes [START_REF] Toprasertpong | Comparison of Electron and Hole Mobilities in Multiple-Quantum-Well Solar Cells Using a Timeof-Flight Technique[END_REF], since electronic transport is a succession of thermal escape, assisted tunnel escape and direct tunneling across barriers. For this SL, ballistic escape time for electrons and holes are estimated to be 7.58 10 -13 s and 4.42 10 -12 s respectively. Since miniband are narrower than the optical phonon energy (ħω=35 meV) the ballistic behavior is expected to be significant.

Our numerical study was made considering a perfectly periodic system with fixed layer widths. However, in a realistic system, layer widths are controlled during the epitaxial growth in order to balance strains in order to avoid dislocations. This leads to induce some small variations in each layer width depending on the strains within the structure and the minimum width that can be grown. In order to take such variations into account, we modeled the same system as shown above, but we added a +0.2, 0 or -0.2nm on each layer width. Variations are randomly chosen for each layer. Fig. 3, which represents the corresponding LDOS, shows that minibands in both CB and VB are not perfectly formed. This is due to the fact that the SL periodicity is altered. In order to demonstrate the importance of the periodicity, we also assumed the randomly generated structure, but with GaAs interlayers widths adjusted in order to keep the SL periodicity. Fig. 4 shows the LDOS calculations for this structure. This system exhibits extended minibands in both CB and VB, which prove that periodicity is a key factor to have extended minibands.

CALCULATIONS IN THE SOLAR CELL

In a SL without any external electric field, this miniband extends across the whole MQW region. However, in a solar cell electric field is required and the LDOS within the MQW region change as well. Fig. 5 shows the LDOS calculations for the 20-period GaAs 0.58 P 0.42 /GaAs/In 0.21 Ga 0.79 As solar cell under a 0.75 V applied bias. The LDOS in CB for this solar cell exhibits a miniband that cross around 13 wells of the SL. Indeed, the electric field tends to shift the energy states in each quantum well, and so, the coupling between wells is not strong enough to allow tunneling process across the whole quantum well region. We note that the second-order miniband in VB is almost not visible.

In order to expand miniband across the whole MQW region, we have modeled a 20-period SL with the GaAs interlayer widths increasing from the n-side to the p-side following a quadratic gradient. The quadratic evolution of the interlayer width permits to align correctly the energy level since the energy states in a quantum well are function of the squared value of the well width. Fig. 6 shows the LDOS calculations under a 0.75 V external bias with interlayer width increasing from 48 Å to 84 Å. The initial value of the interlayer width is higher than the value in the original solar cell (31 Å) to see the evolution of the miniband energy.

In this configuration, the second-order miniband is lower in energy due to the larger interlayer widths. Moreover, the system exhibits a third-order miniband that crosses 18 periods. In this case, more periods are coupled to the contacts. However, there is no miniband in VB. Indeed, if the interlayer width gradient is adapted to expand the miniband in the CB, it is the opposite in VB. This simple gradient is thus adapted either to the CB or the VB. However it can be interested for the collection of carrier generated in the GaAs. If, as shown in Fig. 1 the MQW are located near to the p-doped side of the solar cell, it could be interesting to have an extended miniband in VB, in order to make the holes transport more efficient. In order to expand both CB and VB minibands, we have assumed a double quadratic evolution of the interlayer widths. The interlayer width increases from the p-side to the middle of the SL, and then it decreases from the middle of the SL to the n-side of the solar cell. Fig. 7 shows the LDOS calculated for the double gradient with, from p-side to the middle of the SL to the n-side, 22/53/22 Å interlayers. The solar cell exhibits CB and VB minibands, expanded across 13 and 9 periods respectively.

CONCLUSION

In this paper, we have shown that miniband engineering could help to increase the efficiency of solar cells. Carriers in minibands have lower effective masses than their bulk counterparts. The average velocity of carriers in the miniband is higher than in bulk and strongly higher than a SL without miniband.

We have made LDOS calculations using Green function formalism to show the evidence of miniband formation in a 20-period GaAs 0.58 P 0.42 /GaAs/In 0.21 Ga 0.79 As SL without any electric field in the active region. We showed the importance of the system periodicity on the miniband formation phenomenon.

LDOS calculations for the same solar cell under a 0.75 V applied bias showed that the valence miniband is not expanded across the MQW region, and thus because the internal electric field tends to shift the energy states, reducing the coupling between periods.

Finally, we have proposed a solution to expand both conduction and valence miniband across the SL. By using an interlayer width gradient, we can improve the coupling between more periods of the system and so permit to expand the miniband. However, in order to achieve both conduction and valence miniband improvements, an interlayer width double gradient is needed.

Fig. 1 :

 1 Fig. 1: Band diagram of the modeled solar cell. The SL region is made of alternation of 22 Å GaAs 0.58 P 0.42 barrier, 38 Å In 0.21 Ga 0.79 As well and 31 Å GaAs interlayers between each barrier and well. Red lines represent the CB and blue lines represent the VB.

Fig. 2 :

 2 Fig. 2: LDOS of the MQW region of an undoped InGaAs/GaAs/GaAsP 20-period superlattice made of alternation of 22 Å GaAs 0.58 P 0.42 barrier, 38 Å In 0.21 Ga 0.79 As well and 31 Å GaAs interlayers. The LDOS (color scale on the right) is calculated in arbitrary units, as a function of the energy and the position. a) LDOS calculated in CB. b) LDOS calculated in VB.

Fig. 3 :Fig. 4 :

 34 Fig. 3: LDOS of the MQW region of an undoped InGaAs/GaAs/GaAsP 20-period superlattice with randomly generated errors on each layer width. Errors can be equal to +0.2, +0 or -0.2nm for each layer. The LDOS (color scale on the right) is calculated in arbitrary units, as a function of the energy and the position. a) LDOS calculated in CB. b) LDOS calculated in VB.

Fig. 5 :

 5 Fig. 5: LDOS of the MQW region of an undoped InGaAs/GaAs/GaAsP 20-period superlattice under a 0.75V external bias. The LDOS (color scale on the right) is calculated in arbitrary units, as a function of the energy and the position. a) LDOS calculated in CB. b) LDOS calculated in VB.

Fig. 6 :

 6 Fig. 6: LDOS of the MQW region of an undoped InGaAs/GaAs/GaAsP 20-period superlattice under a 0.75V external bias with an interlayer width gradient from 84 Å to 48 Å. The LDOS (color scale on the right) is calculated in arbitrary units, as a function of the energy and the position. a) LDOS calculated in CB. b) LDOS calculated in VB.

Fig. 7 :

 7 Fig. 7: LDOS of the MQW region of an undoped InGaAs/GaAs/GaAsP 20-period superlattice under a 0.75V external bias with an interlayer width (22/53/22Å) double gradient. The LDOS (color scale on the right) is calculated in arbitrary units, as a function of the energy and the position. a) LDOS calculated in CB. b) LDOS calculated in VB.


G Zheng

D Boutat

J P Barbot

SINGLE OUTPUT DEPENDENT QUADRATIC OBSERVABILITY NORMAL FORM

In this paper, two quadratic observability normal forms dependent on output respectively correspondent to drift dominant term and forced dominant term are firstly discussed. Moreover characteristic numbers are studied in order to simplify the calculation of these two normal forms.

INTRODUCTION

Normal form is a powerful tool to analyze the properties of dynamical systems, such as stability [START_REF] Poincaré | Approximate linearization by state feedback and coordinate change[END_REF], controllability [START_REF] Kang | Extended quadratic controller normal form and dynamic state feedback linearization of non linear systems[END_REF]. In 2001, this concept was firstly introduced in [START_REF] Boutat-Baddas | Quadratic Observability normal form[END_REF] in order to analyze the observability of dynamical system, in which the quadratic observability normal form was studied and the interests of this technique was highlighted. In 2005, a new observability normal form dependent on its output was studied in [START_REF] Zheng | Output Dependent Observability Normal Form[END_REF], and the necessary and sufficient conditions were proposed. Moreover, an extension to multi-outputs case has also been analyzed in [START_REF] Boutat | Observer Error Linearization Multi-Output Depending[END_REF]. Since the observability normal form dependent on its output (or outputs) were exactly linear, another reasonable extension based on the work of [START_REF] Boutat-Baddas | Quadratic Observability normal form[END_REF] and [START_REF] Chabraoui | Observability quadratic characteristic numbers[END_REF] is the study of quadratic observability normal form dependent on its output (or outputs). Therefore, in this paper, we consider the following single input single output system:

ζ = f (ζ) + g(ζ)u y = h(ζ) (1)
where ζ ∈ D ⊂ IR n , u ∈ IR , f : IR n → IR n , g : IR n → IR n and h : IR n → IR are analytic functions, and assume that for all ζ ∈ D, we have rank dh, dL f h, ..., dL n-1

f h T = n.
As an extension of our work in [START_REF] Zheng | Output Dependent Observability Normal Form[END_REF], if the theorem in [START_REF] Zheng | Output Dependent Observability Normal Form[END_REF] is not verified, we assume that system (1) can be transformed into the following form 1 :

   η = β(y) + A(y)η + γ [2] y (η) +B(y)u + ϑ [1] y (η)u + O [3] y (η, u) y = Cη (2) where η ∈ IR n , η = [η 1 , ...η n-1 ] T , γ [2] y (η) = [γ [2] y 1 (η), ..., γ [2] yn (η)] T , ϑ [1] y (η) = [ϑ [1] y 1 (η), ..., ϑ [1] yn (η)] T ,
and for 1 ≤ i ≤ n, γ [2] y i (η) and ϑ [1] y i (η) are function of η with order 2 and 1 respectively parameterized by y, β(y) = [β 1 (y), ..., βn(y)] T , and

A(y) =     0 ... 0 0 α 1 (y) ... 0 0 . . . . . . ... . . . 0 0 α n-1 (y) 0     ,
and B(y) is defined as follows:

B(y) = b(y), 0, • • • , 0 T (3)
As it was shown in [START_REF] Boutat-Baddas | Quadratic Observability normal form[END_REF] that the equivalence modulo an output injection is justified by the fact that the output injection can be canceled in the observation error dynamics. Therefore, by an output injection β(y), system (2) is equivalent modulo an output injection to the following system

   ẋ = A(y)x + γ [2] y (x) +B(y)u + ϑ [1] y (x)u + O [3] y (x, u) y = xn = Cx (4) where x = [x 1 , . . . , x n-1 ] T .
Our problem is how to characterize the fact that all the quadratic terms can be canceled by a diffeomorphism? If this kind of diffeomorphism does not exist, then what's its normal form and its resonant terms?

2. QUADRATIC EQUIVALENCE MODULO AN OUTPUT INJECTION Definition 1. System (4) is quadratically equiva- lent to system    ż = A(y)z + γ[2] y (z) +B(y)u + θ[1] y (z)u + O [3] y (z, u) y = Cz (5)
if there exists a diffeomorphism of the following form:

z = x + φ [2] y (x) (6)
which transforms the quadratic term γ [2] y (x) into another quadratic term γ [2] y (z), where φ

[2] y (x) = [φ [2] y,1 (x), ..., φ [2]
y,n (x)] T and φ [2] y ,i (x) are the homogenous polynomials with order 2 in Remark 1. i) In order to keep the output unchanged, we choose the output equal to xn, which means the diffeomorphism z = x + φ [2] y (x) should verify φ [2] y,n (x) = 0.

ii) It should be noted that this choice is not obligatory. In fact, we can also choose φ [2] y,n (x) = φ

[0] y,n (y), i.e., a function of the output. Proposition 1. System (4) is quadratically equivalent to system (5) modulo an output injection, if and only if the following homologic equations are satisfied:

   γ [2] y (x) + Γ [2] y (x) = A(y)φ [2] y (x) + γ[2] y (x) ϑ [1] y (x) + ∂φ [2] y (x) ∂x 1 b(y) = θ[1] y (x) (7)
where

Γ [2] y (x) = ∂φ [2] y (x) ∂x 2 α 1 (y), • • • , ∂φ [2] y (x) ∂x n-1 α n-2 (y), 0, 0 x. Assume that z = x + φ [2] y (x), hence ż = ẋ + ∂φ [2] y (x) ∂x ẋ.
According to equation ( 4) and ( 5), we have

ż = 1 + ∂φ [2] y (x) ∂x A(y)x + γ [2] y (x) + B(y)u +ϑ [1] y (x)u + O [3] y (x, u) = A(y) x + φ [2] y (x) + γ[2] y (x + φ [2] y (x)) +B(y)u + θ[1] y (x)u + O [3] y (x, u)
So we obtain

     γ [2] y (x) + ∂φ [2] y (x) ∂x A(y)x = A(y)φ [2] y (x) + γ[2] y (x) ϑ [1] y (x) + ∂φ [2] y (x) ∂x B(y) = θ[1] y (x) (8)
where

∂φ [2] y (x) ∂x A(y)x = ∂φ [2] y (x) ∂x 2 α 1 (y), • • • , ∂φ [2] y (x) ∂y α n-1 (y), 0, 0 x = Γ [2] y (x) + O [3] y (x)
and finally equation ( 8) becomes:

   γ [2] y (x) + Γ [2] y (x) = A(y)φ [2] y (x) + γ[2] y (x) ϑ [1] y (x) + ∂φ [2] y (x) ∂x 1 b(y) = θ[1] y (x).

QUADRATIC OBSERVABILITY NORMAL FORM DEPENDENT ON ITS OUTPUT

Since there exist two quadratic terms in system (5): γ[2] y (z) and θ [1] y (z)u, we will study two normal form which correspond respectively to the drift term and forced term in this section.

Normal form with drift dominant term

In this subsection, we study the normal form correspondent to drift dominant term by simplifying the quadratic term γ [2] y (z).

Theorem 1. Normal form correspondent to drift dominant term system (5) by quadratically equivalent modulo an output injection is in the following normal form:

                                 ξ =        n-1 j≥i=1 h ij (y)ξ i ξ j 0 . . . 0        +               n-1 i=1 c1 i (y)ξ i . . . n-1 i=1 cn-1 i (y)ξ i n-1 i=1 cn i (y)ξ i               u +A(y)ξ + B(y)u + O [3] y ( ξ, u) y = ξn (9)
Since the objective of this normal form is to make

γ[2]
y (x) = 0, so the first homologic equation in ( 7) becomes

γ [2] y (x) + Γ [2] y (x) = A(y)φ [2] y (x) Define      φ [2] y (x) = [φ [2] y, 1 (x), ..., φ [2] y,n (x)] T , γ [2] y (x) = [γ [2] y, 1 (x), ..., γ [2] y,n (x)] T , ϑ [1] y (x) = ϑ [1] y, 1 (x), ..., ϑ [1] y,n (x) T , where φ [2] yn (x) = 0.
And we obtain

                         α 1 (y)φ [2] y, 1 (x) = γ [2] y, 2 (x) + n-2 i=1 ∂φ [2] y, 2 (x) ∂x i+1 α i (y)x i . . . α n-2 (y)φ [2] y, n-2 (x) = γ [2] y, n-1 (x) + n-2 i=1 ∂φ [2] y, n-1 (x) ∂x i+1 α i (y)x i α n-1 (y)φ y, n-1 (x) = γ [2] y,n (x) (10)
and the first line of equation ( 7) gives

γ [2] y, 1 (x) + n-2 i=1 ∂φ [2] y, 1 (x) ∂x i+1 x i α i (y) = 0 (11)
Equation ( 10) can be used to deduce φ [2] y, i (x) in order to cancel the quadratic terms from γ

[2] y, 2 (x) to γ [2] y,n (x) respectively. Moreover, if γ [2] y, i (x) and φ [2] y,i (x)
verify also equation ( 11), then this system can be quadratically linearizable. Otherwise, it gives the following resonant terms:

γ [2] y, 1 (x) + n-2 i=1 α i (y) ∂φ [2] y ,1 (x) ∂x i+1 x i . With ϑ [1] y (x)+ ∂φ [2] y (x) ∂x 1 b(y) = θ[1] y (x), we have θ[1] y,n (x) = ϑ [1] y,n (x). Assume θ[1] y,n (x) = n-1 i=1 cn i (y)x i and ϑ [1] y,n (x) = n-1 i=1 c n i (y)x i , then cn i (y) = c n i (y).
And it is not possible to cancel other components. Hence we obtain normal form (9). order to highlight the proposed method, we consider the following example.

Example 1. Consider the following system:

         ẋ1 = x 2 1 -2x 2 3 x 1 x 2 + x 2 3 x 2 2 + x 3 -4x 2 3 x 1 -2x 4 3 x 2 u ẋ2 = x 3 x 1 + x 2 3 x 2 2 -x 3 x 2 + 2x 2 3 x 1 u ẋ3 = x 3 x 2 + x 2 3 x 1 u y = x 3 (12)
where y = 0. According to equation (10), we have

     φ [2] y ,3 (x) = 0 α 2 (y)φ [2] y, 2 (x) = γ [2] y ,3 (x) α 1 (y)φ [2] y, 1 (x) = γ [2] y ,2 (x) + α 1 (y) ∂φ [2] y, 2 (x) ∂x 2
x 1

and we obtain

z 1 = x 1 + x 3 x 2 2 z 2 = x 2 z 3 = x 3
. Hence the quadratic resonant terms are:

γ [2] y, 1 (x) + α 1 (y) ∂φ [2] y, 1 (x) ∂x 2 x 1 = x 2 1 + x 2 3 x 2 2
Finally we have the following normal form:

         ż1 = z 2 1 + z 2 3 z 2 2 + z 3 -4z 2 3 z 1 -2z 4 3 z 2 u +O [3] y (z 1 , z 2 , u) ż2 = z 3 z 1 -z 3 z 2 + 2z 2 3 z 1 u + O [3] y (z 1 , z 2 , u) ż3 = z 3 z 2 + z 2 3 z 1 u + O [3] y (z 1 , z 2 , u) y = z 3 (13)

Normal form with forced dominant term

This subsection is devoted to study another normal form by simplifying quadratic terms: θ [1] y (z)u, in system (5).

Theorem 2. Normal form forced dominant term of system (5) by quadratically equivalent modulo an output injection is in the following form:

                                 ξ =               n-1 j≥i=1 d 1 i (y)ξ i ξ j ξ 1 n-1 i=1 d 2 i (y)ξ i . . . ξ 1 n-1 i=1 d n i (y)ξ i               + u        0 . . . 0 n-1 i=1 cn i (y)ξ i        +A(y)ξ + B(y)u + O [3] y ( ξ, u) y = ξn (14)
Assuming φ [2] y,n (x) = 0, we have

ϑ [1] y,n (x) = θ[1] y,n (x)
and if we set

ϑ [1] y (x) + ∂φ [2] y (x) ∂x 1 b(y) = 0 (15) then we obtain θ[1] y, 1 (x) = ... = θ[1] y, n-1 (x) = 0.
According to the first homologic equation in (7),

γ [2] y (x) + Γ [2] y (x) = A(y)φ [2] y (x) + γ[2] y (x) (16) we obtain                                        γ[2] y, 1 (x) = γ [2] y, 1 (x) + n-2 i=1 ∂φ [2] y, 1 (x) ∂x i+1 x i α i (y) γ[2] y ,2 (x) + α 1 (y)φ [2] y, 1 (x) = γ [2] y, 2 (x) + n-2 i=1 ∂φ [2] y, 2 (x) ∂x i+1 α i (y)x i . . . γ[2] y, n-1 (x) + α n-2 (y)φ [2] y, n-2 (x) = γ [2] y, n-1 (x) + n-2 i=1 ∂φ [2] y, n-1 (x) ∂x i+1 α i (y)x i γ[2] y,n (x) + α n-1 (y)φ [2] y, n-1 (x) = γ [2] y,n (x)
From equation ( 15), the term φ [2] y (x) can be used to cancel all the quadratic terms form the second line to the last one, except for the terms x 1 n-1 i=1

d j i x i , j ∈ [1, n],
hence we obtain:

γ[2] y (x) = n-1 j≥i=1 d 1 i (y)x i x j , x 1 n-1 i=1 d 2 i (y)x i , • • • , x 1 n-1 i=1 d n i (y)x i T
Finally we have the normal form (14).

The following example is to illustrate the proposed normal form.

Example 2. (Example 1 continue) According to the above method, by simple calculation, we have

z 1 = x 1 + 2x 3 x 2 1 + 2x 3 3 x 1 x 2 + x 3 x 2 2 z 2 = x 2 + x 1 x 2 + x 3 x 2 1 z 3 = x 3 (17)
with this diffeomorphism we have the following normal form:

                 ż1 = 1 + 2z 4 3 z 2 1 + z 2 3 z 2 2 z 3 u +O [3] y (z 1 , z 2 , u) ż2 = z 3 z 1 + z 3 -2z 2 3 z 2 1 -2z 4 3 z 1 z 2 +O [3] y (z 1 , z 2 , u) ż3 = z 3 z 2 -z 3 z 1 z 2 -z 2 3 z 2 1 + z 2 3 z 1 u +O [3] y (z 1 , z 2 , u) y = z 3 (18)

CHARACTERISTIC NUMBERS

In order to simplify the calculation of diffeomorphism proposed before, a new method will be proposed in this section which permit us to determine the diffeomorphism (6) in a easer way, with which we need not to solve the homologic equation ( 7).

Characteristic numbers for normal form with drift dominant term

In order to simplify equation ( 7), we assume:

   φ [2] y, i (x) = x T φy, i x, γ[2] y, i (x) = x T γy, i x, γ [2] y, i (x) = x T γy, i x,
and

                                       φy, i :=     φ i 11 (y) • • • φ i 1,n-1 (y) 0 . . . . . . . . . . . . φ i 1,n-1 (y) • • • φ i n-1,n-1 (y) 0 0 • • • 0 0     , γy, i :=     γi 11 (y) • • • γi 1,n-1 (y) 0 . . . . . . . . . . . . γi 1,n-1 (y) • • • γi n-1,n-1 (y) 0 0 • • • 0 0     , γy, i :=     γ i 11 (y) • • • γ i 1,n-1 (y) 0 . . . . . . . . . . . . γ i 1,n-1 (y) • • • γ i n-1,n-1 (y) 0 0 • • • 0 0     .
We obtain

Γ [2] y (x) = [Γ [2] y, 1 (x), ..., Γ [2] y,n (x)] T where Γ [2]
y, i (x) = x T Γy, i x and Γy, i = ĀT (y)φy, i + φy, i Ā(y),

where

Ā(y) :=      0 ... 0 0 0 α 1 (y) ... 0 0 0 . . . . . . . . . . . . . . . 0 ... α n-2 (y) 0 0 0 ... 0 0 0      .
Setting A(y)φ [2] y (x) := x T φyx, for all x, homologic equation ( 7) can be written as follows:

x T γyx + x T Γyx = x T φyx + x T γyx
And we have γy + Γy = φy + γy ( 19)

if γy = 0, we can cancel all the quadratic terms in system (4).

Because

A(y)φ [2] y (x) =     0 α 1 (y)φ [2]
y, 1 (x) . . .

α n-1 (y)φ [2] y, n-1 (x)     if γ[2] y = 0, then     γy, 1 γy, 2 . . . γy, n     + ĀT (y)     φy, 1 φy, 2 . . . φy, n     +     φy, 1 φy, 2 . . . φy, n     Ā(y) = 0, α 1 (y)φy, 1 , • • • , α n-1 (y)φ y,n-1 T we have γy, 1 + ĀT (y)φy, 1 + φy, 1 Ā(y) = 0 α i (y)φy, i = γy, i+1 + ĀT (y)φy, i+1 + φy, i+1 Ā(y), for 1 ≤ i ≤ n -1.
Finally, by recurrence, we get

φy, i = n-1-i k=0          k j=0 C k j ĀT (y) k-j γy, i+k+1 Āj (y) k m=0 α i+m          (20) for 1 ≤ i ≤ n -1, where C k j denotes the combinato- rial coefficient.
With this diffeomorphism, according to the following equality:

ϑ [1] y (x) + ∂φ [2] y (x) ∂x 1 b(y) = θ[1] y (x) we have ϑ [1] y,n (x) = θ[1] y,n (x). Setting            ϑ [1] y,n (x) = n-1 j=1 c n j (y)x j , θ[1] y,n (x) = n-1 j=1 cn j (y)x j , then cn j (y) = c n j (y). Because ϑ [1] y (x) + ∂φ [2] y (x) ∂x 1 b(y) = θ[1] y (x), assuming            ϑ [1] y, i (x) = n-1 j=1 c i j (y)x j , θ[1] y, i (x) = n-1 j=1 ci j (y)x j .
then, since φ [2] y, i (x) = x T φy, i x, we obtain:

ci j (y) = c i j (y) + 2b(y)φ i 1,j (y) , cn j (y) = c n j (y) (21)
for

1 ≤ i, j ≤ n -1.
Definition 2. We define the characteristic matrix for system (4) as follows:

My = γy, 1 + ĀT (y)φy, 1 + φy, 1 Ā(y) (22)
and

Cy =   c1 1 (y) • • • c1 n-1 (y) 0 . . . . . . . . . . . . cn 1 (y) • • • cn n-1 (y) 0   ( 23
)
where My and Cy are only function of y.

We are now able to set the following theorem.

Theorem 3. Normal form correspondent to drift dominant term of system ( 5) is as follows:

             ξ = ξ T     My 0 . . . 0     ξ + A(y)ξ + B(y)u +Cyξu + O [3] y ( ξ, u) y = ξn (24)
where My is defined in (22).

Remark 2. i) My(i, j) in equation ( 24) depends on coefficients h ij (y) of equation ( 9) as follows:

My(i, j) = My(j, i) = 1 2 h ij (y), i < j My(i, i) = h ii (y)
ii) there are 3n(n -1)/2 characteristic numbers in the normal form. . Finally we get the same normal form correspondent to drift dominant term in (13).

Characteristic numbers for normal form with forced dominant term

With the same argument, in order to simplify equation ( 7), setting

           ϑ [1] y,n (x) = n-1 j=1 c n j (y)x j , θ[1] y,n (x) = n-1 j=1 cn j (y)x j ,
and

cn j (y) = c n j (y). ( 25
)
Since ϑ [1] y (x) +

∂φ [2] y (x) ∂x 1 b(y) = θ[1]
y (x), and

φ [2] y, i (x) = x T     φ i 11 (y) • • • φ i 1,n-1 (y) 0 . . . . . . . . . . . . φ i 1,n-1 (y) • • • φ i n-1,n-1 (y) 0 0 • • • 0 0     x, if we set φ i 1,j (y) = - c i j (y) 2b(y) 1, we obtain θ[1] y, 1 (x) = ... = θ[1]
y, n-1 (x) = 0. Then we have γy + Γy = φy + γy, which gives γy, 1 + Ā(y)φy, 1 + φy, 1 Ā(y) = γy, 1 γy, i+1 + α i (y)φy, i = γy, i+1 + ĀT (y)φy, i+1 + φy, i+1 Ā(y)

for 1 ≤ i ≤ n -1. Defining Υy, i =     Υ i 11 (y) • • • Υ i 1,n-1 (y) 0 . . . . . . . . . . . . Υ i 1,n-1 (y) • • • Υ i n-1,n-1 (y) 0 0 • • • 0 0     = n-1-i k=0          k j=0 C k j ĀT (y) k-j γy, i+k+1 -γy, i+k+1 Āj (y) k m=0 α i+m          (26)
and

Ῡy, i =     Ῡi 11 (y) • • • Ῡi 1,n-1 (y) 0 . . . . . . . . . . . . Ῡi 1,n-1 (y) • • • Ῡi n-1,n-1 (y) 0 0 • • • 0 0     (27) = n-1-i k=0          k j=0 C k j ĀT (y) k-j γy, i+k+1 Āj (y) k m=0 α i+m          Because φ i 1,j (y) = - c i j (y) 2b(y) , i, j ∈ [1, n -1], if we set φ i l,s (y) = Ῡi l,s (y) , i ∈ [1, n -1], l, s ∈ [2, n -1] we have γi l,s = 0, i ∈ [2, n], l, s ∈ [2, n -1].
Hence we get the following diffeomorphism:

φy, i : φ i 1,j (y) = - c i j (y) 2b(y) , i, j ∈ [1, n -1] φ i l,s (y) = Ῡi l,s (y) , i ∈ [1, n -1], l, s ∈ [2, n -1] Because we have γi l,s = 0, i ∈ [2, n], l, s ∈ [2, n -1],
then γi 1,j , i ∈ [2, n], j ∈ [1, n -1] can be determined by the following equation:

φ i 1,j (y) = - c i j (y) 2b(y) = Υ i 1,j (y)
and we note

My k (i, j) = My k (j, i) = γk 1,j , i = 1 0, i = 1, n , (28)
for 1 ≤ j ≤ n -1 and 2 ≤ k ≤ n.

Finally because γy, 1 + Ā(y)φy, 1 + φy, 1 Ā(y) = γy ,1 , we obtain:

My ,1 (i, j) = My ,1 (j, i) = γ1 i,j , i < j (29)

Then we can give the following theorem.

Theorem 4. Normal form correspondent to forced dominant term of system ( 5) is as follows: where cn i and My, i are defined by equations ( 25), ( 28) and ( 29).

                     ξ = ξ T     My ,1 My ,2 . . . My, n     ξ +        0 . . .
Remark 3. The number of the free coefficients is 3n(n -1)/2.

CONCLUSION

This paper is devoted to study the quadratic observability normal form parameterized by its output. Above all, two homologic equations were given in order to guarantee the equivalence of quadratic transformation. After that, we have studied two normal forms respectively correspondent to the drift dominant term and forced dominant term. Both representations are equivalent. In order to simplify the calculation, we proposed to apply quadratic terms' characteristic matrix for these two normal forms.

Example 3 .

 3 (Example 1 continue) With the guide of the computation process proposed above, we have My = γy, 1 + ĀT (y)φy, 1 + φy, 1 Ā(y

The necessary and sufficient condition for this transformation is still an open question.


Salma Jiddi
email: salma.jiddi@technicolor.com

Philippe Robert
email: technicolorphilippe.robert@technicolor.com

Eric Marchand
email: eric.marchand@irisa.fr

Photometric Registration using Specular Reflections and Application to Augmented Reality

Keywords: Photometry, Reflectance, Lambertian, Specular, Illumination, Augmented Reality, Relighting. ACM Classification I.4.8. Image Processing and Computer Vision: Scene Analysis-Photometry; I.4.1. Image Processing and Computer Vision: Digitization and Image Capture-Reflectance

Photometric registration consists in blending real and virtual scenes in a visually coherent way. To achieve this goal, both reflectance and illumination properties must be estimated. These estimates are then used, within a rendering pipeline, to virtually simulate the real lighting's interaction with the scene. In this paper, we are interested in indoor scenes where light bounces off of objects with different reflective properties (diffuse and/or specular). In these scenarios, existing solutions often assume distant lighting or limit the analysis to a single specular object. We address scenes with various objects captured by a moving RGB-D camera and estimate the 3D position of light sources. Furthermore, using spatio-temporal data, our algorithm recovers dense diffuse and specular reflectance maps. Finally, using our estimates, we demonstrate photo-realistic augmentations of real scenes (virtual shadows, specular occlusions) as well as virtual specular reflections on real world surfaces.

INTRODUCTION

The image brightness of a 3D object is a function of three components: object geometry, surface reflectance and illumination distribution. Provided these three components are accurately estimated, one is able to photo-realistically render a virtual scene that is visually coherent with a camera's capture of the real scene. This output is of interest for many applications such as Augmented Reality (AR) and scene relighting.

Geometry estimation is a well-known problem in the computer vision community. In the last decade, the field of active RGB-D sensors has experienced significant improvements. Hence, sensors such as Microsoft Kinect, Google Tango and Intel RealSense are nowadays commonly used in photometry registration problems. In our method, we use the R200 sensor that provides a satisfactorily accurate model of the scene.

Problem statement. As far as reflectance and illumination estimation is concerned, existing solutions often address scenes with relatively simple reflectance (e.g. Lambertian and/or non-textured surfaces) and use cast shadows [START_REF] Sato | Illumination from shadows[END_REF], light probes [START_REF] Debevec | Rendering synthetic objects into real scenes: Bridging traditional and image-based graphics with global illumination and high dynamic range photography[END_REF][3] or fisheye cameras [START_REF] Rohmer | Interactive near-field illumination for photorealistic augmented reality on mobile devices[END_REF] to estimate illumination.

On the contrary, when only specular reflections are considered, they are either recovered as highlights using a single frame [5][6] or the analysis is limited to a small region (e.g. a single object) [7][8]. Also, lighting is often assumed to be distant. Nonetheless, in indoor environments where the size of the scene is not small with regard to the distance to light, this approach is likely to deliver inaccurate renderings.

Our contributions. In this paper, we focus on specular reflections to estimate diffuse and specular reflectance properties of indoor scenes with arbitrary geometry and texture. Moreover, we aim at recovering the 3D position of existing light sources without using any light probes or external assistance. We only consider as input the RGB video and 3D geometry acquired with an active sensor. Our estimation is achieved via spatio-temporal processing of the input image sequence and 3D scene analysis. Our key contributions are: (a) diffuse and specular reflectance estimation of 3D objects with arbitrary surface properties; (b) 3D position recovery of light sources responsible for specular reflections.

RELATED WORK

Illumination and reflectance estimation requires fitting a Bidirectional Reflectance Distribution Function (BRDF) and lighting distribution to the color intensity of an image. In this task, also called inverse rendering, several reflection models [START_REF] Phong | Illumination for computer generated pictures[END_REF][10] and rendering techniques [START_REF] Debevec | Rendering synthetic objects into real scenes: Bridging traditional and image-based graphics with global illumination and high dynamic range photography[END_REF] [START_REF] Keller | Instant radiosity[END_REF] have been considered. In the following, we give an overview of related work that takes into account specular reflections. Specular surfaces and distant lighting. Nishino et al. [START_REF] Nishino | Determining reflectance parameters and illumination distribution from a sparse set of images for view-dependent image synthesis[END_REF] separate the reflectance components based on pixel intensity variations throughout a registered image sequence. In [START_REF] Jachnik | Real-time surface light-field capture for augmentation of planar specular surfaces[END_REF], Jachnik et al. capture a surface light field using a hand-held camera. This light field is split into its diffuse component and specular component. The latter is used to recover the environment map. Both [START_REF] Nishino | Determining reflectance parameters and illumination distribution from a sparse set of images for view-dependent image synthesis[END_REF] and [START_REF] Jachnik | Real-time surface light-field capture for augmentation of planar specular surfaces[END_REF] consider a single specular object. Using a fisheye camera to capture lighting, Knecht et al. [START_REF] Knecht | Interactive brdf estimation for mixed-reality applications[END_REF] detect specular reflections as highlight in each frame independently. The scene is segmented using color informa-tion and an albedo value is attributed to each segment. Hence, considered surfaces are mainly textureless.

Specular surfaces and non-distant lighting. In [START_REF] Plopski | Reflectance and light source estimation for indoor ar applications[END_REF], Plopski et al. estimate the position of a dominant light source by recovering the intersection point of reflection rays raised from the center of detected highlights. In [START_REF] Neverova | Lighting estimation in indoor environments from low-quality images[END_REF], Neverova et al. separate an image into its specular and diffuse components then use them as inputs of an optimization process to recover the 3D position of light sources. Both methods use [START_REF] Shen | Simple and efficient method for specularity removal in an image[END_REF] to estimate reflectance. The latter can deliver inaccurate diffuse estimates when the size of the specularity is significant.

PHOTOMETRIC REGISTRATION

We choose Phong reflection model [START_REF] Phong | Illumination for computer generated pictures[END_REF] to describe the way a point p on a surface reflects light as a combination of diffuse and specular reflections:

I p = I p d + I p s
(1) where I p , I p d and I p s are respectively the color, diffuse and specular intensities. Diffuse reflectance is the property that defines an ideal "matte" surface, also called Lambertian surface. Its apparent brightness to an observer is the same regardless of his angle of view. On the other hand, specular reflection is the mirror-like reflection of light from a surface, in which light from a single incoming direction is reflected into a single outgoing direction. Using [START_REF] Phong | Illumination for computer generated pictures[END_REF], diffuse and specular intensities can be described as:

I p = l m=1 [k p d i m (N p • L p m ) + k p s i m (R p m • V p ) αp ] (2)
where k p d is the diffuse reflectance parameter of point p, k p s is its specular reflectance parameter, N p is its normal vector, V p is its viewpoint vector, α p is its shininess parameter, i m is the intensity of the light source m, L p m is its light direction vector, R p m is its reflection vector, and l is the number of light sources present in the scene. Our goal is to estimate these parameters for each 3D point in the scene.

Reflectance Estimation

We take advantage of dense 3D data provided by the sensor to spatio-temporally analyze surface color variations throughout camera trajectory. First, we sub-sample the input image sequence with respect to camera poses in order to avoid redundancies. Then we split the resulting sequence into subsets (this allows us to avoid registration errors resulting from distant camera poses). Finally, we register each subset with regard to a keyframe and provide a spatio-temporal color profile of each 3D point (Fig. 1) represented by a pixel in the image.

In this work, the lighting and scene are supposed to be static. Thus, the intensity variations present in the profile can only originate from the specular component as described in [START_REF] Phong | Illumination for computer generated pictures[END_REF]. In [START_REF] Jachnik | Real-time surface light-field capture for augmentation of planar specular surfaces[END_REF], the median of the profile's values is used to compute I p d . Although this value is robust in presence of shadows or registration errors, it often gives an over-estimation of the diffuse component. We have chosen, as in [START_REF] Nishino | Determining reflectance parameters and illumination distribution from a sparse set of images for view-dependent image synthesis[END_REF], to estimate the diffuse intensity as the minimum observed value in the profile as it should be closer to the correct value (Fig. 2). Moreover, we constrain the recovered diffuse intensity to have a minimum number of occurrences throughout the The specular intensity is retrieved, at each frame t, as the difference between I p (t) and I p d . The unknown specular parameters are then (k p s i m ), R p m and α p (Eq. 2). If pixel intensity at the profile peak is not saturated, the product (k p s i m ) is retrieved as the lobe's peak intensity value. In fact, specular reflection occurs when the vectors R p m and V p are aligned (their scalar product is equal to 1). Hence, we set R p m to be roughly equal to the peak's corresponding view vector V p . Finally, estimated (k p s i m ) values are reported in the chosen keyframe images.

Light Sources 3D Position Estimation

In this section, we aim at estimating the 3D position of light sources represented by point lights. For each keyframe, we use surface normals N p and recovered reflection vectors R p m to compute the light direction vectors (

L p m = 2.(R p m • N p ) • N p -R p m ).
Light rays originating from specular reflections that fall within the same 3D region in the scene are clustered using Euclidean distance. Hence, small clusters are processed as outliers since they generally result from registration errors or inaccurate normals. Then, a mean light direction vector is computed for each significant cluster. Finally, the problem of finding the position of light sources is similar to computing the intersection points of a set of 3D lines.

Optimization

Specular reflectance parameters have been previously recovered only for 3D points roughly viewed along the observed reflection direction. Our goal is to estimate dense reflectance maps. To achieve this, we initially assume uniform specular reflectance for each 3D object in the scene. A first step consists then in clustering the 3D mesh of the scene using Euclidean distance between vertices and normals smoothness constraint [START_REF]Point Cloud Library[END_REF]. Then, provided that each cluster contains at least one 3D point with its specular intensity estimate k s i m , we spread its value to the entire cluster points. Finally, considering all keyframes, we compute the specular reflectance of each object as the maximum encountered k s i m value.

We now address the possibility that an object/cluster may not exhibit a unique specular reflectance. First, we render specular reflections using previously recovered k s i m values and light sources position. Rendered specular maps are correlated with observed specular maps. If correlation fails, we proceeed to a color-based segmentation using the k-means algorithm and set final k s i m values to each color segment as follows: observed points in the direction of light sources that do not exhibit specular effects are considered to be Lambertian (k s = 0), points that are left keep the recovered k s value. Finally, we refine k s i m and estimate α using the Levenberg Marquardt algorithm with the following cost function:

F j = i [I i -(I d,i + m k s i m (R i m • V i (t)) α )] 2 (3)
where i and m respectively iterate over pixels that belong to cluster j and over recovered light sources. I i correspond to observed pixel intensities. In (3), the diffuse component I d,i is fixed and only k s i m and α can be varied by the solver.

EXPERIMENTAL RESULTS

We use the Intel RealSense R200 sensor that provides camera positions and a 3D model of the scene. Also, we photometrically calibrate our sensor as in [START_REF] Robertson | Dynamic range improvement through multiple exposures[END_REF] by taking three images of a color checker with patches of known reflectance at three different shutter speeds. Finally, a calibrated camera with fixed aperture, shutter speed and color gain browses the scene.

Reflectance Evaluation

Figure 4 shows our photometric registration results for three different scenes (S1, S2 and S3). Our algorithm recovers dense diffuse (Fig. 4, col-2) and specular (Fig. 4, col-4) maps.

To evaluate the accuracy of these maps, we achieve two main applications: relighting and Augmented Reality (AR). In Fig. 3, we use recovered reflectance and illumination to virtually relight the specular book. One can notice that the position and intensity of the specularity are well estimated as they are visually coherent with the input image. Furthermore, column 5 in Figure 4 shows different augmented scenes where virtual shadows and specular occlusions are correctly rendered. In Fig. 4, row-1/col-5, one can notice a correct shape of virtual shadows as well as consistent specular occlusions by the virtual green and red spheres. In Figure 4. row-3/col-5, the number of light sources is correctly recovered as well.

Illumination Evaluation

Using a fish-eye lens, we capture the environment map and qualitatively compare it to the recovered lighting distribution (Fig. 5) for three scenes (S1, S2 and S3 with respectively one, two and three light sources). Our algorithm estimates correct lighting distribution as it recovers the correct number of light sources responsible for specular reflections. Furthermore, we use salient control points and evaluate the accuracy of recovered positions. Table 1 shows a comparison between measured (Measured D) and recovered distances to light sources (Estimated D). Our algorithm is tested on various indoor scenes under various lighting (e.g single and multiple spot lights and/or led lights) and recovers light sources positions with an average error of 21cm for a mean distance of 1.95m to the light source.

CONCLUSION AND FUTURE WORK

We presented an algorithm that estimates dense diffuse and specular reflectance maps. Furthermore, our proposed approach recovers the 3D position of light sources for indoor scenes with arbitrary geometry and texture. For future work, we are interested in estimating the color and shape of the illumination in order to produce both hard and soft shadows.

Figure 1 :

 1 Figure 1: Left. 3D/2D registration : t-1 and t+1 are frames registered with regard to keyframe K. Right. Recovered color intensity profile of a point p using a frame subset (red dots) registered with respect to a keyframe K (green dot).

Figure 2 :

 2 Figure 2: Left. Input image of a specular book in the scene. Right. Estimated diffuse reflectance map of the book.

Figure 3 :

 3 Figure 3: Comparison between input image (left) and virtually rendered image using estimated reflectance and lighting.

Figure 4 :

 4 Figure 4: Column 1. Camera views of three indoor scenes (S1, S2 and S3). Column 2. Recovered diffuse maps for surfaces with various textures and reflective properties. Column 3. 3D mesh clustering. Column 4. Estimated specular reflectance parameter k js for each cluster j. Blue pixels correspond to 3D points frequently occluded during scene browsing. Brighter k j s values correspond to more specular surfaces. Column 5. Augmented scenes, with different reflectance and illumination conditions. We demonstrate correctly rendered virtual objects as they occlude real specular effects and show realistic shadows.

Figure 5 :

 5 Figure 5: Left. Captured environment maps using a fish-eye lens. Right. Recovered lighting respectively for S1 (row-1), S2 (row-2) and S3 (row-3) with respectively one, two and three main light sources.

Table 1 .

 1 Comparison between measured and estimated distances to light sources for scenes under different lightings.

	Scenes	Measured D (cm)	Estimated D (cm)
	S1	72.2	-	-	86.1	-	-
	S2	94.7 105.6	-	88.9 127.1	-
	S3	314.7 293.6 306 297.1 261.6 283.2